AliRoot Core  edcc906 (edcc906)
CalibLaserExBscan.C
Go to the documentation of this file.
1 
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
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
49 //
50 TMatrixD *errB[6]; // covariance matrices
51 TMatrixD *errBT[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(){
63 
64  gSystem->Load("libANALYSIS");
65  gSystem->Load("libTPCcalib");
66  gSystem->Load("libSTAT");
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(){
83 
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 
110 
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){
120 
121  Int_t ntracks=3000000;
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(){
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(){
236 
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(){
291 
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 }
TVectorD fitParamA[6]
TCut cutN
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
TOOLKIT for chain manipulation:
Summary of statistics functions.
TVectorD fitParamB[6]
TCut cutCM5
TROOT * gROOT
TMatrixD mat
Definition: AnalyzeLaser.C:9
static TChain * MakeChain(const char *fileIn, const char *treeName, const char *fName=0, Int_t maxFiles=-1, Int_t startFile=0, Int_t checkLevel=0)
Double_t chi2BT0[6]
TTreeSRedirector * pcstream
TChain * chain
void Init()
npoints
Definition: driftITSTPC.C:85
TCut cutAP2
TMatrixD covarB[6]
TCut cutAP5
TMatrixD covarA[6]
void DumpFit()
TMatrixD * errBT0[6]
void DrawPhi()
Double_t chi2
Definition: AnalyzeLaser.C:7
TVectorD fitParamBT[6]
TCut cutB
TCut cutCP5
TMatrixD * errA[6]
TCut cutA
TString * FitPlane(TTree *tree, const char *drawCommand, const char *formula, const char *cuts, Double_t &chi2, Int_t &npoints, TVectorD &fitParam, TMatrixD &covMatrix, Float_t frac=-1, Int_t start=0, Int_t stop=10000000, Bool_t fix0=kFALSE)
TCut cutCP2
void MakeAliases()
TMatrixD covarBT0[6]
TStatToolkit toolkit
Definition: gainCalib.C:18
TCut cutAM5
TCut cutEd
Double_t chi2BT[6]
TVectorD fitParamBT0[6]
TF1 * f
Definition: interpolTest.C:21
TMatrixD covarBT[6]
void MakeFitDy()
TCut cutE
TMatrixD * errBT[6]
void MakeFit(Int_t i, TCut cutI, TString aName)
Double_t chi2A[6]
TCut cutR
void MakeFitPic()
TMatrixD * errB[6]
TMatrixD * MakeErrVector(TMatrixD &mat)
class TVectorT< Double_t > TVectorD
TChain * chainRef
Definition: MakeGlobalFit.C:69
Double_t chi2B[6]
void MakeGraphs()
static TChain * MakeChainRandom(const char *fileIn, const char *treeName, const char *fName=0, Int_t maxFiles=-1, Int_t startFile=0, Int_t checkLevel=0)
class TMatrixT< Double_t > TMatrixD