]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/CalibMacros/CalibLaserExBscan.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliTPCcalibTimeGain.cxx - Adding the Gamma conversion selected electorns
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibLaserExBscan.C
CommitLineData
b1c27e4f 1/*
2//
3 .x ~/rootlogon.C
4 .L $ALICE_ROOT/TPC/CalibMacros/CalibLaserExBscan.C
5
6
7 // 0. Make a calibration
8 // 1. Make a laser scan list
9 // e.g in TPC workscape
10
11 // 2. Define a reference data e.g:
12 // for a in `cat laserScan.txt`; do echo `pwd`/../mergerunMag0.list/laserMean.root; done >laserScanRef.txt
13
14 Init(); // init chain
15 MakeAliases(); // make alaises for variables
16
17
18*/
19
20
21
22
23TCut cutE="eY.fElements<0.02&&Rr.eY.fElements<0.02"; // error cut 200 microns
24TCut cutN="nCl.fElements>10&&Rr.nCl.fElements>10"; // number of clusters cut
25TCut cutEd="(abs(Rr.Y.fElements)-Rr.X.fElements*0.155<-1)"; // edge cut
26TCut cutB="LTr.fVecLX.fElements>0"; // local x cut
27TCut cutR="(LTr.fVecLX.fElements%3)==1&&((abs(LTr.fVecLZ.fElements)+abs(LTr.fVecLY.fElements))%3)==1"; // cut - skip points
28TCut cutA=cutE+cutN+cutEd+cutB;
29//
30TCut cutAM5 =cutA+"LTr.fP[1]>0&&abs(bz+5)<0.1";
31TCut cutAP2 =cutA+"LTr.fP[1]>0&&abs(bz-2)<0.1";
32TCut cutAP5 =cutA+"LTr.fP[1]>0&&abs(bz-5)<0.1";
33//
34TCut cutCM5 =cutA+"LTr.fP[1]<0&&abs(bz+5)<0.1";
35TCut cutCP2 =cutA+"LTr.fP[1]<0&&abs(bz-2)<0.1";
36TCut cutCP5 =cutA+"LTr.fP[1]<0&&abs(bz-5)<0.1";
37
38
39TChain * chain =0;
40
41TVectorD fitParamB[6]; // magnetic fit param
42TVectorD fitParamBT[6]; // + electric field tilt
43TVectorD fitParamBT0[6]; // + electric field close to ROC
44TVectorD fitParamA[6]; // + electric field rotation
45//
46TMatrixD covarB[6]; // covariance matrices
47TMatrixD covarBT[6]; //
48TMatrixD covarBT0[6]; //
49TMatrixD covarA[6]; //
50//
51TMatrixD *errB[6]; // covariance matrices
52TMatrixD *errBT[6]; //
53TMatrixD *errBT0[6]; //
54TMatrixD *errA[6]; //
55//
56Double_t chi2B[6]; // chi2
57Double_t chi2BT[6]; // chi2
58Double_t chi2BT0[6]; //
59Double_t chi2A[6]; //
60
61
62void Init(){
63 //
64 //
65 //
66 gSystem->Load("libANALYSIS");
67 gSystem->Load("libTPCcalib");
68 gSystem->Load("libSTAT.so");
69 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
70 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
71 AliXRDPROOFtoolkit tool;
72 chain = tool.MakeChainRandom("laserScan.txt","Mean",0,10200);
73 chain->Lookup();
74 chainRef = tool.MakeChain("laserScanRef.txt","Mean",0,10200);
75 chainRef->Lookup();
76 chain->AddFriend(chainRef,"Rr");
77}
78
79
80
81void MakeAliases(){
82 //
83 // shortcuts for variables to fit
84 //
85 // bserved distrotions
86 chain->SetAlias("dy","(dY.fElements-Rr.dY.fElements)"); // y ~ (r-phi) distortion
87 chain->SetAlias("ctany","(dY.fElements-Rr.dY.fElements)/(250*dr)"); // mean distortion (angle) over drift length
88 //
89 chain->SetAlias("dr","(250.-abs(LTr.fP[1]))"); // drift length 0-250
90 chain->SetAlias("dr1","(1.-abs(LTr.fP[1]/250.))"); // drift length 0-250
91 chain->SetAlias("phiL","(LTr.fVecLY.fElements/LTr.fVecLX.fElements)"); // tann of local phi
92 chain->SetAlias("Esign","(1.-LTr.fSide*2.)"); // E field sign: +1 for A side : -1 for c side
93 //
94 // Br and Brfi
95 //
96 chain->SetAlias("ablx","(iblx.fElements/(ibz.fElements))");
97 chain->SetAlias("ably","(ibly.fElements/(ibz.fElements))");
98 chain->SetAlias("abr","(ibr.fElements/(ibz.fElements))");
99 chain->SetAlias("abrf","(ibrphi.fElements/(ibz.fElements))");
100
101
102 chain->SetAlias("r","(R.fElements-165.5)/80.3"); // local X - 0 at middle -1 first +1 last padrow
103 chain->SetAlias("ky","(kY.fElements)"); // local inclination angle
104 chain->SetAlias("sec","LTr.fVecSec.fElements"); // sector number
105 chain->SetAlias("ca","cos(LTr.fVecPhi.fElements+0.)"); // cos of global phi position
106 chain->SetAlias("sa","sin(LTr.fVecPhi.fElements+0.)"); // sin of global phi position
107}
108
109
110TMatrixD * MakeErrVector(TMatrixD & mat){
111 //
112 // get error vector
113 //
114 Int_t nrows=mat.GetNrows();
115 TMatrixD *err = new TMatrixD(nrows,1);
116 for (Int_t i=0; i<nrows;i++) (*err)(i,0)=TMath::Sqrt(mat(i,i));
117 return err;
118}
119
120
121void MakeFit(Int_t i, TCut cutI, TString aName){
122 //
123 //
124 //
125 Int_t ntracks=3000000;
126 TStatToolkit toolkit;
127 Double_t chi2=0;
128 Int_t npoints=0;
129 //
130 TString fstringB=""; // magnetic part
131 TString fstringT=""; // E file dtilting part
132 TString fstring0=""; // ROC part
133 TString fstringL=""; // linear part
134 //
135 // Magnetic field map part
136 //
137 // // 0
138 fstringB+="ablx*dr++"; // 1
139 fstringB+="ablx*ky*dr++"; // 2
140 fstringB+="ably*dr++"; // 3
141 fstringB+="ably*ky*dr++"; // 4
142 //
143 // Electric field tilting part
144 //
145 //
146 fstringT+="(dr1)++"; // 5 - ex
147 fstringT+="(dr1)*ky++"; // 6 - ex
148 fstringT+="(dr1)*phiL++"; // 7 - ey
149 fstringT+="(dr1)*phiL*ky++"; // 8 - ey
150 //
151 // E field close to the ROC - radius independent part
152 //
153 //
154 fstring0+="ca++"; // 9
155 fstring0+="sa++"; // 10
156 fstring0+="ky++"; // 11
157 fstring0+="ky*ca++"; // 12
158 fstring0+="ky*sa++"; // 13
159 //
160 // E field close to the ROC - radius dependent part
161 //
162 fstring0+="r++"; // 14
163 fstring0+="ca*r++"; // 15
164 fstring0+="sa*r++"; // 16
165 fstring0+="ky*r++"; // 17
166 fstring0+="ky*ca*r++"; // 18
167 fstring0+="ky*sa*r++"; // 19
168 //
169 // E field rotation - drift length linear part
170 //
171 fstringL+="ca*dr1++"; //20
172 fstringL+="sa*dr1++"; //21
173 fstringL+="ky*ca*dr1++"; //22
174 fstringL+="ky*sa*dr1++"; //23
175 //
176 TString fstringBT = fstringB +fstringT;
177 TString fstringBT0 = fstringBT+fstring0;
178 TString fstringA = fstringBT0+fstringL;
179 //
180 //
181 TCut cutF=cutI;
182 TString * strFit = 0;
183 //
184 //
185 cutF=cutR+cutI;
186 //
187 strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringB.Data(),cutF, chi2,npoints,fitParamB[i],covarB[i],1,0, ntracks);
188 chain->SetAlias(aName+"_FB",strFit->Data());
189 cutF=cutR+cutI + Form("abs(dy-%s)<0.2",(aName+"_FB").Data());
190 //
191 //
192 //
193 strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringB.Data(),cutF, chi2,npoints,fitParamB[i],covarB[i],1,0, ntracks);
194 chain->SetAlias(aName+"_FB",strFit->Data());
195 chi2B[i]=TMath::Sqrt(chi2/npoints);
196 printf("B fit sqrt(Chi2/npoints)=%f\n",chi2B[i]);
197 //
198 strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringBT.Data(),cutF, chi2,npoints,fitParamBT[i],covarBT[i],1,0, ntracks);
199 chain->SetAlias(aName+"_FBT",strFit->Data());
200 chi2BT[i]=TMath::Sqrt(chi2/npoints);
201 printf("BT fit sqrt(Chi2/npoints)=%f\n",chi2BT[i]);
202 //
203 strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringBT0.Data(),cutF, chi2,npoints,fitParamBT0[i],covarBT0[i],1,0, ntracks);
204 chain->SetAlias(aName+"_FBT0",strFit->Data());
205 chi2BT0[i]=TMath::Sqrt(chi2/npoints);
206 printf("BT0 Fit: sqrt(Chi2/npoints)=%f\n",chi2BT0[i]);
207 //
208 //
209 strFit = TStatToolkit::FitPlane(chain,"dy:0.1", fstringA.Data(),cutF, chi2,npoints,fitParamA[i],covarA[i],1, 0, ntracks);
210 chain->SetAlias(aName+"_FA",strFit->Data());
211 chi2A[i]=TMath::Sqrt(chi2/npoints);
212 printf("A fit: sqrt(Chi2/npoints)=%f\n",chi2A[i]);
213 errB[i] =(MakeErrVector(covarB[i]));
214 errBT[i] =(MakeErrVector(covarBT[i]));
215 errBT0[i] =(MakeErrVector(covarBT0[i]));
216 errA[i] =(MakeErrVector(covarA[i]));
217}
218
219
220
221
222
223void MakeFitDy(){
224 //
225 //
226 //
227
228
229 MakeFit(0,cutAM5,"dyAM5");
230 MakeFit(1,cutAP2,"dyAP2");
231 MakeFit(2,cutAP5,"dyAP5");
232 //
233 MakeFit(3,cutCM5,"dyCM5");
234 MakeFit(4,cutCP2,"dyCP2");
235 // MakeFit(5,cutCP5,"dyAP5");
236 DumpFit();
237}
238
239
240void DrawPhi(){
241 //
242 //
243 //
244 chain->Draw("(dy-dyAM5_FA):LTr.fVecPhi.fElements>>hisAM5(60,-3.14,3.14,100,-0.2,0.2)",cutAM5,"");
245 chain->Draw("(dy-dyAP5_FA):LTr.fVecPhi.fElements>>hisAP5(60,-3.14,3.14,100,-0.2,0.2)",cutAP5,"");
246 chain->Draw("(dy-dyAP2_FA):LTr.fVecPhi.fElements>>hisAP2(60,-3.14,3.14,100,-0.2,0.2)",cutAP2,"");
247 hisAM5->FitSlicesY(0,0,-1,20);
248 hisAP5->FitSlicesY(0,0,-1,20);
249 hisAP2->FitSlicesY(0,0,-1,20);
250 hisAM5_1->SetMinimum(-0.2);
251 hisAM5_1->SetMaximum(0.2);
252 hisAM5_1->SetMarkerStyle(20);
253 hisAP5_1->SetMarkerStyle(21);
254 hisAP2_1->SetMarkerStyle(22);
255 hisAM5_1->SetMarkerColor(1);
256 hisAP5_1->SetMarkerColor(2);
257 hisAP2_1->SetMarkerColor(4);
258 hisAM5_1->Draw();
259 hisAP5_1->Draw("same");
260 hisAP2_1->Draw("same");
261}
262
263void MakeGraphs(){
264 Double_t bz[3] ={-5,2,5};
265 Double_t p0A[3]={fitParam[0][0],fitParam[1][0],fitParam[2][0]};
266 Double_t p1A[3]={fitParam[0][1],fitParam[1][1],fitParam[2][1]};
267 Double_t p2A[3]={fitParam[0][2],fitParam[1][2],fitParam[2][2]};
268 Double_t p3A[3]={fitParam[0][3],fitParam[1][3],fitParam[2][3]};
269 Double_t p4A[3]={fitParam[0][4],fitParam[1][4],fitParam[2][4]};
270 Double_t p5A[3]={fitParam[0][5],fitParam[1][5],fitParam[2][5]};
271 Double_t p6A[3]={fitParam[0][6],fitParam[1][6],fitParam[2][6]};
272 Double_t p7A[3]={fitParam[0][7],fitParam[1][7],fitParam[2][7]};
273 Double_t p8A[3]={fitParam[0][8],fitParam[1][8],fitParam[2][8]};
274 Double_t p9A[3]={fitParam[0][9],fitParam[1][9],fitParam[2][9]};
275 TGraph grA0(3,bz,p0A);
276 TGraph grA1(3,bz,p1A);
277 TGraph grA2(3,bz,p2A);
278 TGraph grA3(3,bz,p3A);
279 TGraph grA4(3,bz,p4A);
280 TGraph grA5(3,bz,p5A);
281 TGraph grA6(3,bz,p6A);
282 TGraph grA7(3,bz,p7A);
283 TGraph grA8(3,bz,p8A);
284 TGraph grA9(3,bz,p9A);
285 TF1 f1("f1","[0]*x/(1+([0]*x)^2)");
286 TF1 f2("f2","([0]*x)^2/(1+([0]*x)^2)");
287 TMatrixD matB(5,5);
288 TMatrixD matBT(6,5);
289 TMatrixD matBT0(6,5);
290 for (Int_t i=0; i<5;i++) for (Int_t j=0; j<5;j++) matB[i][j]=fitParamB[j][i];
291 for (Int_t i=0; i<5;i++) for (Int_t j=0; j<5;j++) matBT[i][j]=fitParamBT[j][i];
292 for (Int_t i=0; i<6;i++) for (Int_t j=0; j<5;j++) matB0[i][j]=fitParamBT0[j][i];
293}
294
295
296void DumpFit(){
297 //
298 //
299 //
300 TTreeSRedirector *pcstream = new TTreeSRedirector("exbFits.root");
301
302 for (Int_t i=0; i<5;i++){
303 Double_t bz=0;
304 Double_t side=0;
305 if (i==0) { bz=-5; side=1;}
306 if (i==1) { bz=2; side=1;}
307 if (i==2) { bz=5; side=1;}
308 if (i==3) { bz=-5; side=-1;}
309 if (i==4) { bz=2; side=-1;}
310 (*pcstream)<<"fit"<<
311 "bz="<<bz<<
312 "side="<<side<<
313 "pb.="<<&fitParamB[i]<<
314 "pbT.="<<&fitParamBT[i]<<
315 "pbT0.="<<&fitParamBT0[i]<<
316 "pbA.="<<&fitParamA[i]<<
317 "eb.="<<errB[i]<<
318 "ebT.="<<errBT[i]<<
319 "ebT0.="<<errBT0[i]<<
320 "ebA.="<<errA[i]<<
321 "\n";
322 }
323 delete pcstream;
324}
325
326void MakeFitPic(){
327 TFile f("exbFits.root");
328
329}