4 z = s* (z0 - vd*(t-t0))
8 vd - nominal drift velocity
9 zs - miscalibrated position
11 zs = s*(z0 - vd*(1+vr)*(t-(t0+dt))
12 vr - relative change of the drift velocity
17 zs ~ z - s*vr*(z0-s*z)+s*dzt
18 --------------------------------
19 1. Correction function vr constant:
22 dz = zs-z = -s*vr *(z0-s*z)+s*dzt
23 dzs/dl = dz/dl +s*s*vr*dz/dl
31 gSystem->Load("libSTAT.so");
34 TVectorD fitParam, fitParam1,fitParam2;
35 TVectorD errParam, errParam1,errParam2;
36 TMatrixD covMatrix,covMatrix1,covMatrix2;
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
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
59 TCut cut0("abs(dy)<1&&abs(pTPC.fP[1])>10");
60 //TCut cutOut="(side<0)"
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;
101 npoints =tree->Draw("side",cutA+cut1+cut2);
102 vside.ResizeTo(npoints); vside.SetElements(tree->GetV1());
104 npoints =tree->Draw("z",cutA+cut1+cut2);
105 vz.ResizeTo(npoints); vz.SetElements(tree->GetV1());
107 npoints =tree->Draw("dr",cutA+cut1+cut2);
108 vdr.ResizeTo(npoints); vdr.SetElements(tree->GetV1());
110 npoints =tree->Draw("tl",cutA+cut1+cut2);
111 vtl.ResizeTo(npoints); vtl.SetElements(tree->GetV1());
113 npoints =tree->Draw("sa",cutA+cut1+cut2);
114 vsa.ResizeTo(npoints); vsa.SetElements(tree->GetV1());
116 npoints =tree->Draw("ca",cutA+cut1+cut2);
117 vca.ResizeTo(npoints); vca.SetElements(tree->GetV1());
119 npoints =tree->Draw("dz",cutA+cut1+cut2);
120 vdz.ResizeTo(npoints); vdz.SetElements(tree->GetV1());
122 npoints =tree->Draw("dtl",cutA+cut1+cut2);
123 vdtl.ResizeTo(npoints); vdtl.SetElements(tree->GetV1());
125 npoints =tree->Draw("ez",cutA+cut1+cut2);
126 vez.ResizeTo(npoints); vez.SetElements(tree->GetV1());
128 npoints =tree->Draw("etl*0.5",cutA+cut1+cut2);
129 vetl.ResizeTo(npoints); vetl.SetElements(tree->GetV1());
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);
142 TString fstringTl1="";
143 fstringTl1+="(side)++";
144 fstringTl1+="(tl)++";
145 fstringTl1+="side*(tl)++";
148 TString fstringZ1="";
149 fstringZ1+="(side)++";
150 fstringZ1+="(-side*(250.0-side*z))++";
151 fstringZ1+="((250.-side*z))++";
155 fstringZ1+="side*sa++";
156 fstringZ1+="side*ca++";
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);
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);
169 cut1 = "abs(dz-fitZ1)<2 && abs(dtl-fitTl1)<0.02";
176 TString fstringTl2="";
177 fstringTl2+="(side)++";
178 fstringTl2+="(tl)++";
179 fstringTl2+="(side*tl)++";
181 fstringTl2+="(sa)++";
182 fstringTl2+="(side*sa)++";
183 fstringTl2+="(ca)++";
184 fstringTl2+="(side*ca)++";
186 fstringTl2+="(sa*tl)++";
187 fstringTl2+="(side*sa*tl)++";
188 fstringTl2+="(ca*tl)++";
189 fstringTl2+="(side*ca*tl)++";
194 TString fstringZ2="";
195 fstringZ2+="(side)++";
196 fstringZ2+="(-side*(250.0-side*z))++";
197 fstringZ2+="((250.-side*z))++";
201 fstringZ2+="side*sa++";
202 fstringZ2+="side*ca++";
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);
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";
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");
232 // 3 - side offset P3
233 // 4 - side offset P1
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
239 for (Int_t i=0;i<npoints;i++){
240 for (Int_t jc =0;jc<10;jc++) xxx[jc]=0;
243 xxx[1] = -vside[i]*(250.0-vside[i]*vz[i]);
245 fitter1.AddPoint(xxx,vdz[i],vez[i]);
246 for (Int_t jc =0;jc<10;jc++) xxx[jc]=0;
249 fitter1.AddPoint(xxx,vdtl[i],vetl[i]);
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
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+=")";
262 fstrP3f1+=fitParam1[0];fstrP3f1+="+tl*(";fstrP3f1+=fitParam1[2];
263 fstrP3f1+=")+side*("; fstrP3f1+=fitParam1[3]; fstrP3f1+=")";
265 tree->SetAlias("fP1f1",fstrP1f1->Data());
266 tree->SetAlias("fP3f1",fstrP3f1->Data());
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");
284 // 3 - side offset P3
285 // 4 - side offset P1
287 // 5 - vdr difference -A side - C side
288 // 6 - vdr gradient -sin(alpha)*tl dependence
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
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
297 for (Int_t i=0;i<npoints;i++){
298 for (Int_t jc =0;jc<10;jc++) xxx[jc]=0;
301 xxx[1] = -vside[i]*(250.0-vside[i]*vz[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]);
307 for (Int_t jc =0;jc<10;jc++) xxx[jc]=0;
310 xxx[4] = vside[i]*vtl[i];
311 xxx[5] = vtl[i]*vsa[i];
312 fitter2.AddPoint(xxx,vdtl[i],vetl[i]);
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
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
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+=")";
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+=")";
337 tree->SetAlias("fP1f2",fstrP1f2->Data());
338 tree->SetAlias("fP3f2",fstrP3f2->Data());
363 TString fstringP3="";
367 fstringP3+="side*tl++";
368 fstringP3+="side*sa++";
369 fstringP3+="side*ca++";
371 TString *strP3 = toolkit.FitPlane(tree,"dtl:etl",fstringP3->Data(), cutA+cut38009, chi2,npoints,fitParam,covMatrix);
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);
379 TString fstringP3="";
382 fstringP3+="tl*(run<40000)++";
383 fstringP3+="tl*(abs(run-45600)<200)++";
384 fstringP3+="tl*(run>46000)++";
388 fstringP3+="sa*tl++";
389 fstringP3+="ca*tl++";
391 fstringP3+="side*tl++";
392 fstringP3+="side*sa++";
393 fstringP3+="side*ca++";
394 fstringP3+="side*sa*tl++";
395 fstringP3+="side*ca*tl++";
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);
407 TString fstringP1="";
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)++";
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)++";
426 fstringP1+="sa*tl++";
427 fstringP1+="ca*tl++";
429 fstringP1+="side*sa++";
430 fstringP1+="side*ca++";
431 fstringP1+="side*sa*tl++";
432 fstringP1+="side*ca*tl++";
434 TCut cuttl="abs(dtl-fitP3)<0.02"
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);