ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/HFmon2012/macros/fitData.C
Revision: 1.1
Committed: Mon Aug 13 14:57:42 2012 UTC (12 years, 8 months ago) by yilmaz
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Error occurred while calculating annotation data.
Log Message:
update analysis

File Contents

# Content
1 #include "Laser.h"
2 #include "TLegend.h"
3
4 double getChi2(TFile* mcFile, TProfile* p, int i = 3, int j = 2, int k = 3, bool draw = 0){
5
6 double xMin = 0.5;
7 double xMax = 0.8;
8
9 // cout<<"i : "<<i<<" j : "<<j<<" k : "<<k<<endl;
10
11 TProfile* pmc = (TProfile*)mcFile->Get(Form("p_%d_%d_%d",i,j,k));
12 pmc->SetMarkerStyle(24);
13 pmc->SetMarkerColor(2);
14 pmc->SetLineColor(2);
15
16 if(draw){
17
18 p->SetTitle(";phase/25ns+arbitrary const;R");
19 p->Draw();
20 pmc->Draw("same");
21 TLegend* leg = new TLegend(0.52,0.73,0.91,0.91);
22 leg->SetFillColor(0);
23 leg->SetBorderSize(1);
24 leg->SetFillStyle(1001);
25 leg->SetTextFont(43);
26 leg->SetTextSize(13);
27 leg->AddEntry(p,"Data","pl");
28 leg->AddEntry(pmc,"Model","pl");
29 leg->Draw();
30 }
31
32 double sum = 0;
33 int ndof = 0;
34
35 for(int ib = 1; ib < p->GetNbinsX(); ++ib){
36
37 double x = p->GetBinCenter(ib);
38 if(x < xMin || x > xMax) continue;
39
40 double y1 = p->GetBinContent(ib);
41 double y2 = pmc->GetBinContent(ib);
42 double e1 = p->GetBinError(ib);
43 double e2 = pmc->GetBinError(ib);
44 double ee = e1*e1 + e2*e2;
45 double dy = (y2-y1);
46 double d = dy*dy/ee;
47
48 sum+=d;
49 ndof++;
50 }
51
52 if(sum > 0 && ndof > 0) sum /= ndof;
53
54 return sum;
55
56 }
57
58
59
60
61 double analyzeChannel(TNtuple* nt, int ieta = 30, int iphi = 21, int idepth = 1, int phase0 = 1066, const char* mcLabel = "12_6_5"){
62
63 double deltaT = 23;
64 nt->SetAlias("R","sts3/sts4");
65 nt->SetAlias("P",Form("(p-%d)/%f",phase0,deltaT));
66
67
68 double r[200] = { 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.97, 0.98, 1. };
69 double w[100] = { 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5 };
70 double n[100] = { 0.5, 1, 2, 5, 10, 20, 100 };
71
72 int Nr = 100;
73 int Nw = 9;
74 int Nn = 7;
75
76 for(int i = 0; i < Nr; ++i){
77 r[i] = 1./Nr*(i+1);
78 }
79
80 r[Nr] = 1.01;
81 w[Nw] = 0.55;
82 n[Nn] = 105;
83
84 Nn = 3;
85 Nw = 9;
86
87 TFile* outf = new TFile("nothing.root","recreate");
88 TH2D* hChi2 = new TH2D("hChi2",";R;W;#Xi^{2}",Nr,r,Nw,w);
89
90 TH2D* h = new TH2D(Form("h"),"",200,-2,6,200,-0.2,1.2);
91 TProfile* p = new TProfile(Form("p"),"",200,-2,6);
92 TH2D* hw = new TH2D(Form("hw"),"",200,-2,6,200,-0.2,1.2);
93 TProfile* pw = new TProfile(Form("pw"),"",200,-2,6);
94
95 TH2D* hrw = new TH2D(Form("hrw"),"",200,-2,6,200,-0.2,1.2);
96 TProfile* prw = new TProfile(Form("prw"),"",200,-2,6);
97
98 TProfile* wdata = new TProfile(Form("wdata"),"",20,0,1);
99 TProfile* wmc = new TProfile(Form("wmc"),"",20,0,1);
100
101
102 TFile* mcFile = new TFile("resultsD.root");
103 outf->cd();
104
105
106 TCut channel(Form("ieta == %d && iphi == %d && idepth == %d",ieta,iphi,idepth));
107 TCut rlimits("R > -0.2 && R < 1.2");
108 TCut wing(" w > 0.95");
109
110 nt->Draw("R:P>>h",channel&&rlimits);
111 nt->Draw("R:P>>p",channel&&rlimits);
112
113 nt->Draw("R:P>>hw",channel&&rlimits&&wing);
114 nt->Draw("R:P>>pw",channel&&rlimits&&wing);
115
116 nt->Draw("w:P>>hrw",channel&&rlimits);
117 nt->Draw("w:P>>prw",channel&&rlimits);
118
119
120 TCanvas* c1 = new TCanvas("c1","c1",500,500);
121
122 p->SetMarkerStyle(20);
123 p->SetMarkerColor(1);
124 p->SetLineColor(1);
125
126 p->Draw();
127
128 wdata->SetMarkerStyle(20);
129 wdata->SetMarkerColor(1);
130 wdata->SetLineColor(1);
131 wmc->SetMarkerStyle(24);
132 wmc->SetMarkerColor(2);
133 wmc->SetLineColor(2);
134
135 int kNoise = 5;
136
137 for(int i = 0; i < Nr; ++i){
138 for(int j = 0; j < Nw; ++j){
139 double chi2 = getChi2(mcFile,p,i,j,kNoise);
140 hChi2->Fill(r[i],w[j],chi2);
141 }
142 }
143
144 hChi2->Draw("lego2");
145
146 TCanvas* c2 = new TCanvas("c2","c2",500,500);
147 int ibest = 0;
148 int jbest = 0;
149 int kbest = kNoise;
150
151 int run = 180633;
152
153 hChi2->GetMinimumBin(ibest,jbest,kbest);
154
155 getChi2(mcFile,p,ibest,jbest,kbest,1);
156
157 TNtuple* ntmc = (TNtuple*)mcFile->Get(Form("nt_%d_%d_%d",ibest,jbest,kbest));
158
159 TLegend* leg = new TLegend(0.52,0.53,0.91,0.75);
160 leg->SetFillColor(0);
161 leg->SetBorderSize(1);
162 leg->SetFillStyle(1001);
163 leg->SetTextFont(43);
164 leg->SetTextSize(13);
165 leg->AddEntry(hChi2,Form("Run %d",run),"");
166 leg->AddEntry(hChi2,Form("i#eta,i#phi,depth = (%d,%d,%d)",ieta,iphi,idepth),"");
167 leg->AddEntry(hChi2,Form("Fit R = %0.2f",r[ibest]),"");
168 leg->AddEntry(hChi2,Form("Fit pulse width (ns) = %0.2f",w[jbest]),"");
169 leg->Draw();
170
171 c2->Print(Form("FitResult_run%d_ieta%d_iphi%d_idepth%d.gif",run,ieta,iphi,idepth));
172
173
174 TCanvas* c3 = new TCanvas("c3","c3",500,500);
175
176 nt->Draw("w:P>>wdata",channel&&rlimits&&" P >0 && P < 1","");
177 ntmc->Draw("w:p>>wmc","r < 1.8 && r > -0.2 && p > 0 && p < 1","same");
178
179 wmc->Draw();
180 wdata->Draw("same");
181
182 return 0;
183 }
184
185
186
187
188
189
190
191 void fitData(int run = 180633){
192
193 TFile *inf = new TFile(Form("res%d.root",run));
194 TNtuple* nt = (TNtuple*)inf->Get("nt");
195
196 int phase0 = 1066;
197 // analyzeChannel(nt,38,21,1,1066);
198
199 // analyzeChannel(nt,30,1,2,1062);
200
201 analyzeChannel(nt,30,37,2,1059);
202
203 // analyzeChannel(nt,30,57,1,1061);
204 // analyzeChannel(nt,32,1,1,1061);
205
206 // analyzeChannel(nt,32,21,2,1061);
207 // analyzeChannel(nt,32,37,1,1061);
208
209 // analyzeChannel(nt,32,57,2,1062);
210
211 // analyzeChannel(nt,34,21,1,1060);
212
213
214
215 }
216