TPC/CalibMacros/CalibLaserExBscan.C

   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 /// ~~~
  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(){
  62   ///
  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(){
  80   /// shortcuts for variables to fit
  81   ///
  82   /// bserved distrotions
  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
 108 TMatrixD * MakeErrVector(TMatrixD & mat){
 109   /// get error vector
 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){
 119   ///
 120
 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(){
 220   ///
 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(){
 235   ///
 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(){
 290   ///
 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 }