AliRoot Core  a565103 (a565103)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
driftITSTPC.C
Go to the documentation of this file.
1 
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 
31 void Init(){
32 
33  gSystem->Load("libSTAT");
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 }
57 
58 
59 
60 void 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 }
81 
82 //
83 // variable array
84 //
85 npoints=0;
86 TVectorD vside;
87 TVectorD vz;
88 TVectorD vdr;
89 TVectorD vtl;
90 TVectorD vsa;
91 TVectorD vca;
92 TVectorD vdz;
93 TVectorD vdtl;
94 TVectorD vez;
95 TVectorD vetl;
96 
97 TVectorD errors;
98 
99 
100 void FillVar(){
102 
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 }
138 
139 
140 void FitI1(){
142 
143  TString fstringTl1="";
144  fstringTl1+="(side)++";
145  fstringTl1+="(tl)++";
146  fstringTl1+="side*(tl)++";
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  //
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);
170  cut1 = "abs(dz-fitZ1)<2 && abs(dtl-fitTl1)<0.02";
171 }
172 
173 void Fit2I(){
175 
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 
220 void Fit1(){
224 
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 
271 void Fit2(){
275 
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 
347 
348 
349 
350 
351 
352 
353 
354 
355 
356 
357 
358 
359 
360 
361 TString fstringP3="";
362 fstringP3+="tl++";
363 fstringP3+="sa++";
364 fstringP3+="ca++";
365 fstringP3+="side*tl++";
366 fstringP3+="side*sa++";
367 fstringP3+="side*ca++";
368 
369 TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA+cut38009, chi2,npoints,fitParam,covMatrix);
370 strP3->Tokenize("+")->Print();
371 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
372 tree->SetAlias("fitP3",strP3->Data());
373 tree->Draw("dtl-fitP3",cut0+cut38009);
374 
375 
376 
377 TString fstringP3="";
378 //
379 fstringP3+="side++";
380 fstringP3+="tl*(run<40000)++";
381 fstringP3+="tl*(abs(run-45600)<200)++";
382 fstringP3+="tl*(run>46000)++";
383 //
384 fstringP3+="sa++";
385 fstringP3+="ca++";
386 fstringP3+="sa*tl++";
387 fstringP3+="ca*tl++";
388 //
389 fstringP3+="side*tl++";
390 fstringP3+="side*sa++";
391 fstringP3+="side*ca++";
392 fstringP3+="side*sa*tl++";
393 fstringP3+="side*ca*tl++";
394 
395 
396 
397 TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA, chi2,npoints,fitParam,covMatrix);
398 strP3->Tokenize("+")->Print();
399 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
400 tree->SetAlias("fitP3",strP3->Data());
401 tree->Draw("dtl-fitP3",cut0);
402 
403 
404 
405 TString fstringP1="";
406 //
407 fstringP1+="side++";
408 fstringP1+="dr*(run<40000)++";
409 fstringP1+="dr*(abs(run-45600)<200)++";
410 fstringP1+="dr*(run>46000)++";
411 fstringP1+="side*dr*(run<40000)++";
412 fstringP1+="side*dr*(abs(run-45600)<200)++";
413 fstringP1+="side*dr*(run>46000)++";
414 //
415 fstringP1+="tl*(run<40000)++";
416 fstringP1+="tl*(abs(run-45600)<200)++";
417 fstringP1+="tl*(run>46000)++";
418 fstringP1+="side*tl*(run<40000)++";
419 fstringP1+="side*tl*(abs(run-45600)<200)++";
420 fstringP1+="side*tl*(run>46000)++";
421 //
422 fstringP1+="sa++";
423 fstringP1+="ca++";
424 fstringP1+="sa*tl++";
425 fstringP1+="ca*tl++";
426 //
427 fstringP1+="side*sa++";
428 fstringP1+="side*ca++";
429 fstringP1+="side*sa*tl++";
430 fstringP1+="side*ca*tl++";
431 
432 TCut cuttl="abs(dtl-fitP3)<0.02"
433 
434 TString *strP1 = toolkit.FitPlane(tree,"dz",fstringP1->Data(),cuttl+cutA, chi2,npoints,fitParam,covMatrix);
435 strP1->Tokenize("+")->Print();
436 printf("Chi2/npoints = %f\n",TMath::Sqrt(chi2/npoints));
437 tree->SetAlias("fitP1",strP1->Data());
438 tree->Draw("dz-fitP1",cuttl+cut0);
439 
440 */
printf
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
Fit2I
void Fit2I()
Definition: driftITSTPC.C:173
vdtl
TVectorD vdtl
Definition: driftITSTPC.C:93
vetl
TVectorD vetl
Definition: driftITSTPC.C:95
vez
TVectorD vez
Definition: driftITSTPC.C:94
f
TFile f("CalibObjects.root")
fstringP3
TString fstringP3
Definition: driftITSTPC.C:361
vz
TVectorD vz
Definition: driftITSTPC.C:87
tree
TTree * tree
Definition: AnalyzeLaserCE.C:160
npoints
npoints
Definition: driftITSTPC.C:85
vdr
TVectorD vdr
Definition: driftITSTPC.C:88
vdz
TVectorD vdz
Definition: driftITSTPC.C:92
strP3
TString * strP3
Definition: driftITSTPC.C:369
Fit2
void Fit2()
Definition: driftITSTPC.C:271
InitCuts
void InitCuts()
Definition: driftITSTPC.C:60
chi2
Double_t chi2
Definition: AnalyzeLaser.C:7
fstringP1
TString fstringP1
Definition: driftITSTPC.C:405
Init
void Init()
Definition: driftITSTPC.C:31
cutA
TCut cutA
Definition: CalibLaserExBscan.C:27
toolkit
TStatToolkit toolkit
Definition: gainCalib.C:18
cuttl
TCut cuttl
Definition: driftITSTPC.C:432
errors
TVectorD errors
Definition: driftITSTPC.C:97
Fit1
void Fit1()
Definition: driftITSTPC.C:220
vca
TVectorD vca
Definition: driftITSTPC.C:91
FitI1
void FitI1()
Definition: driftITSTPC.C:140
vsa
TVectorD vsa
Definition: driftITSTPC.C:90
FillVar
void FillVar()
Definition: driftITSTPC.C:100
vtl
TVectorD vtl
Definition: driftITSTPC.C:89
vside
TVectorD vside
Definition: driftITSTPC.C:86