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