]> 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
Update variable description file
[u/mrichter/AliRoot.git] / TPC / CalibMacros / driftITSTPC.C
CommitLineData
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 29void 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 58void 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//
83npoints=0;
84TVectorD vside;
85TVectorD vz;
86TVectorD vdr;
87TVectorD vtl;
88TVectorD vsa;
89TVectorD vca;
90TVectorD vdz;
91TVectorD vdtl;
92TVectorD vez;
93TVectorD vetl;
94
95TVectorD errors;
96
97
98void 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 138void 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 172void 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 220void 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
272void 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
363TString fstringP3="";
364fstringP3+="tl++";
365fstringP3+="sa++";
366fstringP3+="ca++";
367fstringP3+="side*tl++";
368fstringP3+="side*sa++";
369fstringP3+="side*ca++";
370
0e943608 371TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA+cut38009, chi2,npoints,fitParam,covMatrix);
34aa7526 372strP3->Tokenize("+")->Print();
373printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
374tree->SetAlias("fitP3",strP3->Data());
375tree->Draw("dtl-fitP3",cut0+cut38009);
376
377
378
379TString fstringP3="";
380//
381fstringP3+="side++";
382fstringP3+="tl*(run<40000)++";
383fstringP3+="tl*(abs(run-45600)<200)++";
384fstringP3+="tl*(run>46000)++";
385//
386fstringP3+="sa++";
387fstringP3+="ca++";
388fstringP3+="sa*tl++";
389fstringP3+="ca*tl++";
390//
391fstringP3+="side*tl++";
392fstringP3+="side*sa++";
393fstringP3+="side*ca++";
394fstringP3+="side*sa*tl++";
395fstringP3+="side*ca*tl++";
396
397
398
399TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA, chi2,npoints,fitParam,covMatrix);
400strP3->Tokenize("+")->Print();
401printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
402tree->SetAlias("fitP3",strP3->Data());
403tree->Draw("dtl-fitP3",cut0);
404
405
406
407TString fstringP1="";
408//
409fstringP1+="side++";
410fstringP1+="dr*(run<40000)++";
411fstringP1+="dr*(abs(run-45600)<200)++";
412fstringP1+="dr*(run>46000)++";
413fstringP1+="side*dr*(run<40000)++";
414fstringP1+="side*dr*(abs(run-45600)<200)++";
415fstringP1+="side*dr*(run>46000)++";
416//
417fstringP1+="tl*(run<40000)++";
418fstringP1+="tl*(abs(run-45600)<200)++";
419fstringP1+="tl*(run>46000)++";
420fstringP1+="side*tl*(run<40000)++";
421fstringP1+="side*tl*(abs(run-45600)<200)++";
422fstringP1+="side*tl*(run>46000)++";
423//
424fstringP1+="sa++";
425fstringP1+="ca++";
426fstringP1+="sa*tl++";
427fstringP1+="ca*tl++";
428//
429fstringP1+="side*sa++";
430fstringP1+="side*ca++";
431fstringP1+="side*sa*tl++";
432fstringP1+="side*ca*tl++";
433
434TCut cuttl="abs(dtl-fitP3)<0.02"
435
436TString *strP1 = toolkit.FitPlane(tree,"dz",fstringP1->Data(),cuttl+cutA, chi2,npoints,fitParam,covMatrix);
437strP1->Tokenize("+")->Print();
438printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
439tree->SetAlias("fitP1",strP1->Data());
440tree->Draw("dz-fitP1",cuttl+cut0);
441
0e943608 442*/