]>
Commit | Line | Data |
---|---|---|
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 | |
22 | TCut cutE="eY.fElements<0.02&&Rr.eY.fElements<0.02"; // error cut 200 microns | |
23 | TCut cutN="nCl.fElements>10&&Rr.nCl.fElements>10"; // number of clusters cut | |
24 | TCut cutEd="(abs(Rr.Y.fElements)-Rr.X.fElements*0.155<-1)"; // edge cut | |
25 | TCut cutB="LTr.fVecLX.fElements>0"; // local x cut | |
26 | TCut cutR="(LTr.fVecLX.fElements%3)==1&&((abs(LTr.fVecLZ.fElements)+abs(LTr.fVecLY.fElements))%3)==1"; // cut - skip points | |
27 | TCut cutA=cutE+cutN+cutEd+cutB; | |
28 | // | |
29 | TCut cutAM5 =cutA+"LTr.fP[1]>0&&abs(bz+5)<0.1"; | |
30 | TCut cutAP2 =cutA+"LTr.fP[1]>0&&abs(bz-2)<0.1"; | |
31 | TCut cutAP5 =cutA+"LTr.fP[1]>0&&abs(bz-5)<0.1"; | |
32 | // | |
33 | TCut cutCM5 =cutA+"LTr.fP[1]<0&&abs(bz+5)<0.1"; | |
34 | TCut cutCP2 =cutA+"LTr.fP[1]<0&&abs(bz-2)<0.1"; | |
35 | TCut cutCP5 =cutA+"LTr.fP[1]<0&&abs(bz-5)<0.1"; | |
36 | ||
37 | ||
38 | TChain * chain =0; | |
39 | ||
40 | TVectorD fitParamB[6]; // magnetic fit param | |
41 | TVectorD fitParamBT[6]; // + electric field tilt | |
42 | TVectorD fitParamBT0[6]; // + electric field close to ROC | |
43 | TVectorD fitParamA[6]; // + electric field rotation | |
44 | // | |
45 | TMatrixD covarB[6]; // covariance matrices | |
46 | TMatrixD covarBT[6]; // | |
47 | TMatrixD covarBT0[6]; // | |
48 | TMatrixD covarA[6]; // | |
49 | // | |
50 | TMatrixD *errB[6]; // covariance matrices | |
51 | TMatrixD *errBT[6]; // | |
52 | TMatrixD *errBT0[6]; // | |
53 | TMatrixD *errA[6]; // | |
54 | // | |
55 | Double_t chi2B[6]; // chi2 | |
56 | Double_t chi2BT[6]; // chi2 | |
57 | Double_t chi2BT0[6]; // | |
58 | Double_t chi2A[6]; // | |
59 | ||
60 | ||
61 | void 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 | ||
79 | void 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 | ||
108 | TMatrixD * 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 | ||
118 | void 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 | ||
219 | void 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 | ||
234 | void 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 | ||
256 | void 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 | ||
289 | void 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 | ||
318 | void MakeFitPic(){ | |
319 | TFile f("exbFits.root"); | |
320 | ||
321 | } |