TPC/CalibMacros/CalibLaserExBscan.C

   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
  23 TCut  cutE="eY.fElements<0.02&&Rr.eY.fElements<0.02";         // error cut 200 microns
  24 TCut  cutN="nCl.fElements>10&&Rr.nCl.fElements>10";           // number of clusters cut
  25 TCut  cutEd="(abs(Rr.Y.fElements)-Rr.X.fElements*0.155<-1)";   // edge cut
  26 TCut  cutB="LTr.fVecLX.fElements>0";                          // local x cut
  27 TCut  cutR="(LTr.fVecLX.fElements%3)==1&&((abs(LTr.fVecLZ.fElements)+abs(LTr.fVecLY.fElements))%3)==1";   // cut - skip points
  28 TCut  cutA=cutE+cutN+cutEd+cutB;
  29 //
  30 TCut  cutAM5  =cutA+"LTr.fP[1]>0&&abs(bz+5)<0.1";
  31 TCut  cutAP2  =cutA+"LTr.fP[1]>0&&abs(bz-2)<0.1";
  32 TCut  cutAP5  =cutA+"LTr.fP[1]>0&&abs(bz-5)<0.1";
  33 //
  34 TCut  cutCM5  =cutA+"LTr.fP[1]<0&&abs(bz+5)<0.1";
  35 TCut  cutCP2  =cutA+"LTr.fP[1]<0&&abs(bz-2)<0.1";
  36 TCut  cutCP5  =cutA+"LTr.fP[1]<0&&abs(bz-5)<0.1";
  37
  38
  39 TChain * chain =0;
  40
  41 TVectorD fitParamB[6];      // magnetic fit param
  42 TVectorD fitParamBT[6];     // + electric field tilt
  43 TVectorD fitParamBT0[6];    // + electric field close to ROC
  44 TVectorD fitParamA[6];      // + electric field rotation
  45 //
  46 TMatrixD covarB[6];         // covariance matrices
  47 TMatrixD covarBT[6];        //
  48 TMatrixD covarBT0[6];       //
  49 TMatrixD covarA[6];         //
  50 //
  51 TMatrixD *errB[6];         // covariance matrices
  52 TMatrixD *errBT[6];        //
  53 TMatrixD *errBT0[6];       //
  54 TMatrixD *errA[6];         //
  55 //
  56 Double_t chi2B[6];          // chi2
  57 Double_t chi2BT[6];         // chi2
  58 Double_t chi2BT0[6];        //
  59 Double_t chi2A[6];          //
  60
  61
  62 void Init(){
  63   //
  64   //
  65   //
  66   gSystem->Load("libANALYSIS");
  67   gSystem->Load("libTPCcalib");
  68   gSystem->Load("libSTAT");
  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
  81 void 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
 110 TMatrixD * 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
 121 void 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
 223 void 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
 240 void 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
 263 void 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
 296 void 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
 326 void MakeFitPic(){
 327   TFile f("exbFits.root");
 328
 329 }