]>
Commit | Line | Data |
---|---|---|
34aa7526 | 1 | /* |
2 | Formulas: | |
3 | ||
4 | z = s* (z0 - vd*(t-t0)) | |
5 | ||
6 | s - side -1 and +1 | |
7 | t0 - time 0 | |
8 | vd - nominal drift velocity | |
9 | zs - miscalibrated position | |
10 | ||
11 | zs = s*(z0 - vd*(1+vr)*(t-(t0+dt)) | |
12 | vr - relative change of the drift velocity | |
13 | dzt - vd*dt | |
14 | dr = zz0-s*z | |
15 | .. | |
16 | ==> | |
17 | zs ~ z - s*vr*(z0-s*z)+s*dzt | |
18 | -------------------------------- | |
19 | 1. Correction function vr constant: | |
20 | ||
21 | ||
22 | dz = zs-z = -s*vr *(z0-s*z)+s*dzt | |
23 | dzs/dl = dz/dl +s*s*vr*dz/dl | |
24 | d(dz/dl) = vr*dz/dl | |
25 | ||
26 | */ | |
27 | ||
28 | ||
0e943608 | 29 | void Init(){ |
30 | ||
31 | gSystem->Load("libSTAT.so"); | |
32 | TStatToolkit toolkit; | |
33 | Double_t chi2; | |
34 | TVectorD fitParam, fitParam1,fitParam2; | |
35 | TVectorD errParam, errParam1,errParam2; | |
36 | TMatrixD covMatrix,covMatrix1,covMatrix2; | |
37 | ||
38 | Int_t npoints; | |
39 | // | |
40 | TFile f("driftitsTPC.root"); | |
41 | TTree * tree = (TTree*)f.Get("Test"); | |
42 | tree->SetAlias("side","(-1+(pTPC.fP[1]>0)*2)"); //side | |
43 | tree->SetAlias("z","(pTPC.fP[1]+0.0)"); //z position | |
44 | tree->SetAlias("dr","(1-abs(pTPC.fP[1])/250.)"); //norm drift length | |
45 | tree->SetAlias("tl","(pTPC.fP[3]+pITS.fP[3])*0.5"); //tan lampbda | |
46 | tree->SetAlias("sa","sin(pTPC.fAlpha+0.)"); //sin alpha | |
47 | tree->SetAlias("ca","cos(pTPC.fAlpha+0.)"); //cos alpha | |
48 | ||
49 | tree->SetAlias("dz","(pTPC.fP[1]-pITS.fP[1])"); //z delta | |
50 | tree->SetAlias("dy","(pTPC.fP[0]-pITS.fP[0])"); //z delta | |
51 | tree->SetAlias("dtl","(pTPC.fP[3]-pITS.fP[3])"); //delta tan lampbda | |
52 | tree->SetAlias("etl","sqrt(pITS.fC[9]+0.)"); //error tan lampbda | |
53 | tree->SetAlias("ez","sqrt(pITS.fC[2]+0.)"); //error z | |
54 | } | |
34aa7526 | 55 | |
56 | ||
57 | ||
0e943608 | 58 | void InitCuts(){ |
59 | TCut cut0("abs(dy)<1&&abs(pTPC.fP[1])>10"); | |
60 | //TCut cutOut="(side<0)" | |
61 | TCut cutA=cut0; | |
62 | TCut cut38009("run==38009"); // 524 | |
63 | TCut cut38010("run==38010"); // 60 | |
64 | TCut cut38286("run==38286"); | |
65 | TCut cut38532("run==38532"); | |
66 | TCut cut38576("run==38576"); | |
67 | TCut cut38586("run==38586"); | |
68 | TCut cut38591("run==38591"); | |
69 | TCut cut45639("run==45639"); | |
70 | TCut cut46993("run==46993"); | |
71 | TCut cut47164("run==47164"); | |
72 | TCut cut47175("run==47175"); | |
73 | TCut cut47180("run==47180"); | |
74 | TCut cut47274("run==47274"); | |
75 | TCut cutA = cut38591+cut0; | |
76 | TCut cut1 = "1"; | |
77 | TCut cut2 = "1"; | |
78 | } | |
34aa7526 | 79 | |
0e943608 | 80 | // |
81 | // variable array | |
82 | // | |
83 | npoints=0; | |
84 | TVectorD vside; | |
85 | TVectorD vz; | |
86 | TVectorD vdr; | |
87 | TVectorD vtl; | |
88 | TVectorD vsa; | |
89 | TVectorD vca; | |
90 | TVectorD vdz; | |
91 | TVectorD vdtl; | |
92 | TVectorD vez; | |
93 | TVectorD vetl; | |
94 | ||
95 | TVectorD errors; | |
96 | ||
97 | ||
98 | void FillVar(){ | |
99 | // | |
100 | // | |
101 | npoints =tree->Draw("side",cutA+cut1+cut2); | |
102 | vside.ResizeTo(npoints); vside.SetElements(tree->GetV1()); | |
103 | // | |
104 | npoints =tree->Draw("z",cutA+cut1+cut2); | |
105 | vz.ResizeTo(npoints); vz.SetElements(tree->GetV1()); | |
106 | // | |
107 | npoints =tree->Draw("dr",cutA+cut1+cut2); | |
108 | vdr.ResizeTo(npoints); vdr.SetElements(tree->GetV1()); | |
109 | // | |
110 | npoints =tree->Draw("tl",cutA+cut1+cut2); | |
111 | vtl.ResizeTo(npoints); vtl.SetElements(tree->GetV1()); | |
112 | // | |
113 | npoints =tree->Draw("sa",cutA+cut1+cut2); | |
114 | vsa.ResizeTo(npoints); vsa.SetElements(tree->GetV1()); | |
115 | // | |
116 | npoints =tree->Draw("ca",cutA+cut1+cut2); | |
117 | vca.ResizeTo(npoints); vca.SetElements(tree->GetV1()); | |
118 | // | |
119 | npoints =tree->Draw("dz",cutA+cut1+cut2); | |
120 | vdz.ResizeTo(npoints); vdz.SetElements(tree->GetV1()); | |
121 | // | |
122 | npoints =tree->Draw("dtl",cutA+cut1+cut2); | |
123 | vdtl.ResizeTo(npoints); vdtl.SetElements(tree->GetV1()); | |
124 | // | |
125 | npoints =tree->Draw("ez",cutA+cut1+cut2); | |
126 | vez.ResizeTo(npoints); vez.SetElements(tree->GetV1()); | |
127 | // | |
128 | npoints =tree->Draw("etl*0.5",cutA+cut1+cut2); | |
129 | vetl.ResizeTo(npoints); vetl.SetElements(tree->GetV1()); | |
130 | // | |
131 | npoints = tree->Draw("run","run>0"); | |
132 | Int_t *np=new Int_t[npoints]; for (Int_t i=0;i<npoints;i++) np[i]=tree->GetV1()[i]; | |
133 | Int_t *nsort=new Int_t[npoints]; | |
134 | Int_t nruns = toolkit.Freq(npoints,np,nsort,kTRUE); | |
135 | } | |
34aa7526 | 136 | |
34aa7526 | 137 | |
0e943608 | 138 | void FitI1(){ |
139 | // | |
140 | // Independent fit 1 | |
141 | // | |
34aa7526 | 142 | TString fstringTl1=""; |
143 | fstringTl1+="(side)++"; | |
144 | fstringTl1+="(tl)++"; | |
0e943608 | 145 | fstringTl1+="side*(tl)++"; |
34aa7526 | 146 | // |
147 | // | |
148 | TString fstringZ1=""; | |
149 | fstringZ1+="(side)++"; | |
150 | fstringZ1+="(-side*(250.0-side*z))++"; | |
151 | fstringZ1+="((250.-side*z))++"; | |
152 | // systematic shift | |
153 | fstringZ1+="sa++"; | |
154 | fstringZ1+="ca++"; | |
155 | fstringZ1+="side*sa++"; | |
156 | fstringZ1+="side*ca++"; | |
157 | // | |
34aa7526 | 158 | TString *strTl1 = toolkit.FitPlane(tree,"dtl:etl",fstringTl1->Data(), cutA+cut1, chi2,npoints,fitParam,covMatrix); |
159 | strTl1->Tokenize("+")->Print(); | |
160 | printf("Tl1: Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); | |
161 | tree->SetAlias("fitTl1",strTl1->Data()); | |
162 | tree->Draw("dtl-fitTl1",cutA+cut1); | |
163 | // | |
164 | TString *strZ1 = toolkit.FitPlane(tree,"dz:ez",fstringZ1->Data(), cutA+cut1, chi2,npoints,fitParam,covMatrix); | |
165 | strZ1->Tokenize("+")->Print(); | |
166 | printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); | |
167 | tree->SetAlias("fitZ1",strZ1->Data()); | |
168 | tree->Draw("dz-fitZ1",cutA+cut1); | |
0e943608 | 169 | cut1 = "abs(dz-fitZ1)<2 && abs(dtl-fitTl1)<0.02"; |
34aa7526 | 170 | } |
171 | ||
0e943608 | 172 | void Fit2I(){ |
173 | // | |
174 | // Independent fit | |
175 | // | |
34aa7526 | 176 | TString fstringTl2=""; |
177 | fstringTl2+="(side)++"; | |
178 | fstringTl2+="(tl)++"; | |
179 | fstringTl2+="(side*tl)++"; | |
180 | // | |
181 | fstringTl2+="(sa)++"; | |
182 | fstringTl2+="(side*sa)++"; | |
183 | fstringTl2+="(ca)++"; | |
184 | fstringTl2+="(side*ca)++"; | |
185 | // | |
186 | fstringTl2+="(sa*tl)++"; | |
187 | fstringTl2+="(side*sa*tl)++"; | |
188 | fstringTl2+="(ca*tl)++"; | |
189 | fstringTl2+="(side*ca*tl)++"; | |
190 | // | |
191 | // | |
192 | // | |
193 | // | |
194 | TString fstringZ2=""; | |
195 | fstringZ2+="(side)++"; | |
196 | fstringZ2+="(-side*(250.0-side*z))++"; | |
197 | fstringZ2+="((250.-side*z))++"; | |
198 | // systematic shift | |
199 | fstringZ2+="sa++"; | |
200 | fstringZ2+="ca++"; | |
201 | fstringZ2+="side*sa++"; | |
202 | fstringZ2+="side*ca++"; | |
203 | // | |
204 | ||
205 | TString *strTl2 = toolkit.FitPlane(tree,"dtl:etl",fstringTl2->Data(), cutA+cut1+cut2, chi2,npoints,fitParam,covMatrix); | |
206 | strTl2->Tokenize("+")->Print(); | |
207 | printf("Tl2: Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); | |
208 | tree->SetAlias("fitTl2",strTl2->Data()); | |
209 | tree->Draw("dtl-fitTl2",cutA+cut1); | |
210 | // | |
211 | TString *strZ2 = toolkit.FitPlane(tree,"dz:ez",fstringZ2->Data(), cutA+cut1+cut2, chi2,npoints,fitParam,covMatrix); | |
212 | strZ2->Tokenize("+")->Print(); | |
213 | printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); | |
214 | tree->SetAlias("fitZ2",strZ2->Data()); | |
215 | tree->Draw("dz-fitZ2",cutA+cut1); | |
216 | cut2 = "abs(dz-fitZ2)<2 && abs(dtl-fitTl2)<0.02"; | |
217 | } | |
218 | ||
219 | ||
0e943608 | 220 | void Fit1(){ |
221 | // | |
222 | // | |
223 | // dz = zs-z = -s*vr *(z0-s*z)+s*dzt | |
224 | // dzs/dl = dz/dl +s*s*vr*dz/dl | |
225 | // d(dz/dl) = vr*dz/dl | |
226 | TLinearFitter fitter1(5, "hyp4"); | |
227 | // parameters | |
228 | // 0 - P3 offset | |
229 | // 1 - P1 offset | |
230 | // 2 - vdr factor | |
231 | // | |
232 | // 3 - side offset P3 | |
233 | // 4 - side offset P1 | |
234 | ||
235 | // dz = zs-z = -s*vr *(z0-s*z)+s*dzt = [1] + [2]*(-s*(z0-s*z))+ [4]*s | |
236 | // d(dz/dl) = vr*dz/dl = [0] + [2]*dz/dl + [3]*s | |
237 | ||
238 | Double_t xxx[100]; | |
239 | for (Int_t i=0;i<npoints;i++){ | |
240 | for (Int_t jc =0;jc<10;jc++) xxx[jc]=0; | |
241 | // P1 | |
242 | xxx[0] = 1; | |
243 | xxx[1] = -vside[i]*(250.0-vside[i]*vz[i]); | |
244 | xxx[3] = vside[i]; | |
245 | fitter1.AddPoint(xxx,vdz[i],vez[i]); | |
246 | for (Int_t jc =0;jc<10;jc++) xxx[jc]=0; | |
247 | xxx[1] = vtl[i]; | |
248 | xxx[2] = vside[i]; | |
249 | fitter1.AddPoint(xxx,vdtl[i],vetl[i]); | |
250 | } | |
251 | fitter1.Eval(); | |
252 | fitter1.GetParameters(fitParam1); | |
253 | fitter1.GetErrors(errParam1); | |
254 | // dz = zs-z = -s*vr *(z0-s*z)+s*dzt = [1] + [2]*(-s*(z0-s*z))+ [4]*s | |
255 | // d(dz/dl) = vr*dz/dl = [0] + [2]*dz/dl + [3]*s | |
256 | ||
257 | TString fstrP1f1 =""; | |
258 | TString fstrP3f1 =""; | |
259 | fstrP1f1+=fitParam1[1];fstrP1f1+="-side*(250.0-side*z)*(";fstrP1f1+=fitParam1[2]; | |
260 | fstrP1f1+=")+side*("; fstrP1f1+=fitParam1[4]; fstrP1f1+=")"; | |
261 | // | |
262 | fstrP3f1+=fitParam1[0];fstrP3f1+="+tl*(";fstrP3f1+=fitParam1[2]; | |
263 | fstrP3f1+=")+side*("; fstrP3f1+=fitParam1[3]; fstrP3f1+=")"; | |
264 | // | |
265 | tree->SetAlias("fP1f1",fstrP1f1->Data()); | |
266 | tree->SetAlias("fP3f1",fstrP3f1->Data()); | |
267 | // | |
268 | } | |
269 | ||
270 | ||
271 | ||
272 | void Fit2(){ | |
273 | // | |
274 | // | |
275 | // dz = zs-z = -s*vr *(z0-s*z)+s*dzt | |
276 | // dzs/dl = dz/dl +s*s*vr*dz/dl | |
277 | // d(dz/dl) = vr*dz/dl | |
278 | TLinearFitter fitter2(7, "hyp6"); | |
279 | // parameters | |
280 | // 0 - P3 offset | |
281 | // 1 - P1 offset | |
282 | // 2 - vdr factor | |
283 | // | |
284 | // 3 - side offset P3 | |
285 | // 4 - side offset P1 | |
286 | // | |
287 | // 5 - vdr difference -A side - C side | |
288 | // 6 - vdr gradient -sin(alpha)*tl dependence | |
289 | ||
290 | ||
291 | // dz = zs-z = -s*vr *(z0-s*z)+s*dzt = [1] + [2]*(-s*(z0-s*z))+ [4]*s +[5]* (-s*s*(z0-s*z)) | |
292 | // dz+= [6]*(-s*(z0-s*z))*sa | |
293 | // | |
294 | // d(dz/dl) = vr*dz/dl = [0] + [2]*dz/dl + [3]*s +[5]* s*dz/dl | |
295 | // d(dz/dl)+= [6]*dz/dl*sa | |
296 | Double_t xxx[100]; | |
297 | for (Int_t i=0;i<npoints;i++){ | |
298 | for (Int_t jc =0;jc<10;jc++) xxx[jc]=0; | |
299 | // P1 | |
300 | xxx[0] = 1; | |
301 | xxx[1] = -vside[i]*(250.0-vside[i]*vz[i]); | |
302 | xxx[3] = vside[i]; | |
303 | xxx[4] = -vside[i]*vside[i]*(250.0-vside[i]*vz[i]); | |
304 | xxx[5] = -vside[i]*(250.0-vside[i]*vz[i])*vsa[i]; | |
305 | fitter2.AddPoint(xxx,vdz[i],vez[i]); | |
306 | //P3 | |
307 | for (Int_t jc =0;jc<10;jc++) xxx[jc]=0; | |
308 | xxx[1] = vtl[i]; | |
309 | xxx[2] = vside[i]; | |
310 | xxx[4] = vside[i]*vtl[i]; | |
311 | xxx[5] = vtl[i]*vsa[i]; | |
312 | fitter2.AddPoint(xxx,vdtl[i],vetl[i]); | |
313 | } | |
314 | fitter2.Eval(); | |
315 | fitter2.GetParameters(fitParam2); | |
316 | fitter2.GetErrors(errParam2); | |
317 | // dz = zs-z = -s*vr *(z0-s*z)+s*dzt = [1] + [2]*(-s*(z0-s*z))+ [4]*s +[5]* (-s*s*(z0-s*z)) | |
318 | // dz+= [6]*(-s*(z0-s*z))*sa | |
319 | // | |
320 | // d(dz/dl) = vr*dz/dl = [0] + [2]*dz/dl + [3]*s +[5]* s*dz/dl | |
321 | // d(dz/dl)+= [6]*dz/dl*sa | |
322 | ||
323 | ||
324 | TString fstrP1f2 =""; | |
325 | TString fstrP3f2 =""; | |
326 | fstrP1f2+=fitParam2[1];fstrP1f2+="-side*(250.0-side*z)*(";fstrP1f2+=fitParam2[2]; | |
327 | fstrP1f2+=")+side*("; fstrP1f2+=fitParam2[4]; | |
328 | fstrP1f2+=")-side*side*(250.0-side*z)*("; fstrP1f2+=fitParam2[5]; | |
329 | fstrP1f2+=")-side*(250.0-side*z)*sa*tl*(";fstrP1f2+=fitParam2[6];fstrP1f2+=")"; | |
330 | // | |
331 | fstrP3f2+=fitParam2[0];fstrP3f2+="+tl*(";fstrP3f2+=fitParam2[2]; | |
332 | fstrP3f2+=")+side*("; fstrP3f2+=fitParam2[3]; | |
333 | fstrP3f2+=")+side*tl*("; fstrP3f2+=fitParam2[5]; | |
334 | fstrP3f2+=")*tl*sa*("; fstrP3f2+=fitParam2[6];fstrP3f2+=")"; | |
335 | ||
336 | // | |
337 | tree->SetAlias("fP1f2",fstrP1f2->Data()); | |
338 | tree->SetAlias("fP3f2",fstrP3f2->Data()); | |
339 | // | |
340 | } | |
341 | ||
342 | ||
343 | ||
344 | ||
345 | ||
346 | ||
347 | ||
348 | ||
34aa7526 | 349 | |
350 | ||
351 | ||
352 | ||
353 | ||
354 | ||
355 | ||
356 | ||
357 | ||
358 | ||
359 | ||
360 | ||
361 | ||
362 | ||
363 | TString fstringP3=""; | |
364 | fstringP3+="tl++"; | |
365 | fstringP3+="sa++"; | |
366 | fstringP3+="ca++"; | |
367 | fstringP3+="side*tl++"; | |
368 | fstringP3+="side*sa++"; | |
369 | fstringP3+="side*ca++"; | |
370 | ||
0e943608 | 371 | TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA+cut38009, chi2,npoints,fitParam,covMatrix); |
34aa7526 | 372 | strP3->Tokenize("+")->Print(); |
373 | printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); | |
374 | tree->SetAlias("fitP3",strP3->Data()); | |
375 | tree->Draw("dtl-fitP3",cut0+cut38009); | |
376 | ||
377 | ||
378 | ||
379 | TString fstringP3=""; | |
380 | // | |
381 | fstringP3+="side++"; | |
382 | fstringP3+="tl*(run<40000)++"; | |
383 | fstringP3+="tl*(abs(run-45600)<200)++"; | |
384 | fstringP3+="tl*(run>46000)++"; | |
385 | // | |
386 | fstringP3+="sa++"; | |
387 | fstringP3+="ca++"; | |
388 | fstringP3+="sa*tl++"; | |
389 | fstringP3+="ca*tl++"; | |
390 | // | |
391 | fstringP3+="side*tl++"; | |
392 | fstringP3+="side*sa++"; | |
393 | fstringP3+="side*ca++"; | |
394 | fstringP3+="side*sa*tl++"; | |
395 | fstringP3+="side*ca*tl++"; | |
396 | ||
397 | ||
398 | ||
399 | TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA, chi2,npoints,fitParam,covMatrix); | |
400 | strP3->Tokenize("+")->Print(); | |
401 | printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); | |
402 | tree->SetAlias("fitP3",strP3->Data()); | |
403 | tree->Draw("dtl-fitP3",cut0); | |
404 | ||
405 | ||
406 | ||
407 | TString fstringP1=""; | |
408 | // | |
409 | fstringP1+="side++"; | |
410 | fstringP1+="dr*(run<40000)++"; | |
411 | fstringP1+="dr*(abs(run-45600)<200)++"; | |
412 | fstringP1+="dr*(run>46000)++"; | |
413 | fstringP1+="side*dr*(run<40000)++"; | |
414 | fstringP1+="side*dr*(abs(run-45600)<200)++"; | |
415 | fstringP1+="side*dr*(run>46000)++"; | |
416 | // | |
417 | fstringP1+="tl*(run<40000)++"; | |
418 | fstringP1+="tl*(abs(run-45600)<200)++"; | |
419 | fstringP1+="tl*(run>46000)++"; | |
420 | fstringP1+="side*tl*(run<40000)++"; | |
421 | fstringP1+="side*tl*(abs(run-45600)<200)++"; | |
422 | fstringP1+="side*tl*(run>46000)++"; | |
423 | // | |
424 | fstringP1+="sa++"; | |
425 | fstringP1+="ca++"; | |
426 | fstringP1+="sa*tl++"; | |
427 | fstringP1+="ca*tl++"; | |
428 | // | |
429 | fstringP1+="side*sa++"; | |
430 | fstringP1+="side*ca++"; | |
431 | fstringP1+="side*sa*tl++"; | |
432 | fstringP1+="side*ca*tl++"; | |
433 | ||
434 | TCut cuttl="abs(dtl-fitP3)<0.02" | |
435 | ||
436 | TString *strP1 = toolkit.FitPlane(tree,"dz",fstringP1->Data(),cuttl+cutA, chi2,npoints,fitParam,covMatrix); | |
437 | strP1->Tokenize("+")->Print(); | |
438 | printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints)); | |
439 | tree->SetAlias("fitP1",strP1->Data()); | |
440 | tree->Draw("dz-fitP1",cuttl+cut0); | |
441 | ||
0e943608 | 442 | */ |