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