PWGPP/TPC/macros/MakeReportTPC.C

   1 /*
   2    //gSystem->AddIncludePath("-I$ALICE_ROOT/PWGPP/TPC");
   3    //.L $ALICE_ROOT/PWGPP/TPC/macros/MakeReportTPC.C+
   4    //
   5    //MakeReportTPC();
   6
   7    gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
   8    gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
   9    AliXRDPROOFtoolkit tool;
  10    TChain * chain = tool.MakeChain("summaryTPCQA.txt","tpcQA",0,200000);
  11
  12    TTree *tree = chain->CopyTree("1");
  13 */
  14
  15
  16 #if !defined(__CINT__) || defined(__MAKECINT__)
  17 #include "TSystem.h"
  18 #include "TMath.h"
  19 #include "TVectorD.h"
  20 #include "TFile.h"
  21 #include "TGrid.h"
  22 #include "TF1.h"
  23 #include "TH1.h"
  24 #include "TH2.h"
  25 #include "TProfile.h"
  26 #include "THnSparse.h"
  27 #include "TTreeStream.h"
  28 #include "AliPerformanceTPC.h"
  29 #include "AliPerformanceDEdx.h"
  30 #endif
  31
  32 TObject *pTPCObject=0;
  33 TObject *pTPCObjectGain=0;
  34 TTreeSRedirector  *pcstream=0;
  35 void Init();
  36 void ReadObjects(const char * path=0);
  37 void MakeReportTPC(Int_t run, const char *path=0 );
  38 void AnalyzeNCL();
  39 void AnalyzeDrift();
  40 void AnalyzeGain();
  41
  42 void Init(){
  43   //
  44   //
  45   //
  46   gSystem->Load("libANALYSIS.so");
  47   gSystem->Load("libANALYSISalice.so");
  48   gSystem->Load("libTender.so");
  49   gSystem->Load("libCORRFW.so");
  50   gSystem->Load("libPWG0base.so");
  51   gSystem->Load("libPWG0dep.so");
  52   gSystem->Load("libPWG0selectors.so");
  53   gSystem->Load("libPWGPP.so");
  54   gSystem->Load("libTPCcalib");
  55 //   gSystem->Setenv("alien_CLOSE_SE","ALICE::GSI::SE");
  56 //   TGrid::Connect("alien://",0,0,"t");
  57 //   gSystem->Setenv("alien_CLOSE_SE","ALICE::GSI::SE");
  58 }
  59
  60
  61 void ReadObjects(const char * path){
  62   //
  63   //
  64   //
  65   TFile *f =0;
  66   if (path) f=TFile::Open(Form("%s/TPC.PerformanceQAQA.root",path));
  67   if (!path) f=TFile::Open("TPC.Performance.root");
  68   if (!f){
  69     printf("File %s not available\n", path);
  70     return;
  71   }
  72   TList *list = (TList*)f->Get("TPC");
  73   if (!list){
  74     printf("QA %s not available\n", path);
  75     return;
  76   }
  77   pTPCObject= (AliPerformanceTPC*)list->FindObject("AliPerformanceTPC");
  78   pTPCObjectGain = (AliPerformanceDEdx*)list->FindObject("AliPerformanceDEdxTPCInner");
  79 }
  80
  81
  82 void MakeReportTPC(Int_t run, const char *path ){
  83   //
  84   // make a tpcQA report
  85   //  typical variables exported to the tree
  86   //
  87   Init();
  88   ReadObjects(path);
  89
  90   pcstream= new TTreeSRedirector("TPC.PerformanceSummary.root");
  91   (*pcstream)<<"tpcQA"<<"run="<<run;
  92   AnalyzeNCL();
  93   AnalyzeDrift();
  94   AnalyzeGain();
  95   (*pcstream)<<"tpcQA"<<"\n";
  96   delete pcstream;
  97 }
  98
  99
 100 void AnalyzeNCL(){
 101   //
 102   // NCL statistic
 103   //
 104   // variables:
 105   static Double_t meanTPCnclF=0;
 106   static Double_t rmsTPCnclF=0;
 107   static Double_t slopeATPCnclF=0;
 108   static Double_t slopeCTPCnclF=0;
 109   static Double_t slopeATPCnclFErr=0;
 110   static Double_t slopeCTPCnclFErr=0;
 111   static Double_t meanTPCncl=0;
 112   static Double_t rmsTPCncl=0;
 113   static Double_t slopeATPCncl=0;
 114   static Double_t slopeCTPCncl=0;
 115   static Double_t slopeATPCnclErr=0;
 116   static Double_t slopeCTPCnclErr=0;
 117   AliPerformanceTPC * pTPC=  (AliPerformanceTPC *)pTPCObject;
 118
 119   TH1* his1D=0;
 120   TProfile* hprof=0;
 121   static TF1 *fpol1 = new TF1("fpol1","pol1");
 122   //
 123   // all clusters
 124   // eta cut - +-1
 125   // pt cut  - +-0.250 GeV
 126   pTPC->GetTPCTrackHisto()->GetAxis(5)->SetRangeUser(-1.,1.);
 127   pTPC->GetTPCTrackHisto()->GetAxis(7)->SetRangeUser(0.25,10);
 128   his1D = pTPC->GetTPCTrackHisto()->Projection(0);
 129   meanTPCncl= his1D->GetMean();
 130   rmsTPCncl= his1D->GetRMS();
 131   delete his1D;
 132   hprof = pTPC->GetTPCTrackHisto()->Projection(0,5)->ProfileX();
 133   hprof->Fit(fpol1,"QNR","QNR",0,0.8);
 134   slopeATPCncl= fpol1->GetParameter(1);
 135   slopeATPCnclErr= fpol1->GetParError(1);
 136   hprof->Fit(fpol1,"QNR","QNR",-0.8,0.0);
 137   slopeCTPCncl= fpol1->GetParameter(1);
 138   slopeCTPCnclErr= fpol1->GetParameter(1);
 139   delete hprof;
 140   //
 141   // findable clusters
 142   //
 143   pTPC->GetTPCTrackHisto()->GetAxis(2)->SetRangeUser(0.4,1.1);
 144   his1D = pTPC->GetTPCTrackHisto()->Projection(2);
 145   meanTPCnclF= his1D->GetMean();
 146   rmsTPCnclF= his1D->GetRMS();
 147   delete his1D;
 148   his1D = pTPC->GetTPCTrackHisto()->Projection(2,5)->ProfileX();
 149   his1D->Fit(fpol1,"QNR","QNR",0,0.8);
 150   slopeATPCnclF= fpol1->GetParameter(1);
 151   slopeATPCnclFErr= fpol1->GetParError(1);
 152   his1D->Fit(fpol1,"QNR","QNR",-0.8,0.0);
 153   slopeCTPCnclF= fpol1->GetParameter(1);
 154   slopeCTPCnclFErr= fpol1->GetParameter(1);
 155   delete his1D;
 156   printf("Cluster QA report\n");
 157   printf("meanTPCnclF=\t%f\n",meanTPCnclF);
 158   printf("rmsTPCnclF=\t%f\n",rmsTPCnclF);
 159   printf("slopeATPCnclF=\t%f\n",slopeATPCnclF);
 160   printf("slopeCTPCnclF=\t%f\n",slopeCTPCnclF);
 161   printf("meanTPCncl=\t%f\n",meanTPCncl);
 162   printf("rmsTPCncl=\t%f\n",rmsTPCncl);
 163   printf("slopeATPCncl=\t%f\n",slopeATPCncl);
 164   printf("slopeCTPCncl=\t%f\n",slopeCTPCncl);
 165   //
 166   // dump results to the tree
 167   //
 168   (*pcstream)<<"tpcQA"<<
 169     "meanTPCnclF="<<meanTPCnclF <<
 170     "rmsTPCnclF="<<rmsTPCnclF <<
 171     "slopeATPCnclF="<< slopeATPCnclF<<
 172     "slopeCTPCnclF="<< slopeCTPCnclF<<
 173     "slopeATPCnclFErr="<< slopeATPCnclFErr<<
 174     "slopeCTPCnclFErr="<< slopeCTPCnclFErr<<
 175     "meanTPCncl="<<meanTPCncl <<
 176     "rmsTPCncl="<< rmsTPCncl<<
 177     "slopeATPCncl="<< slopeATPCncl<<
 178     "slopeCTPCncl="<< slopeCTPCncl<<
 179     "slopeATPCnclErr="<< slopeATPCnclErr<<
 180     "slopeCTPCnclErr="<< slopeCTPCnclErr;
 181
 182
 183 }
 184
 185 void AnalyzeDrift(){
 186   //
 187   // Analyze drift velocity imperfections
 188   //
 189   // variables:
 190   static Double_t offsetdZA=0;
 191   static Double_t slopedZA=0;
 192   static Double_t offsetdZC=0;
 193   static Double_t slopedZC=0;
 194   static Double_t offsetdZAErr=0;
 195   static Double_t slopedZAErr=0;
 196   static Double_t offsetdZCErr=0;
 197   static Double_t slopedZCErr=0;
 198   static Double_t offsetdZAchi2=0;
 199   static Double_t slopedZAchi2=0;
 200   static Double_t offsetdZCchi2=0;
 201   static Double_t slopedZCchi2=0;
 202   AliPerformanceTPC * pTPC=  (AliPerformanceTPC *)pTPCObject;
 203
 204   TH1* his1D=0;
 205   TH2* his2D=0;
 206   static TF1 *fpol1 = new TF1("fpol1","pol1");
 207   TObjArray arrayFit;
 208   his2D = pTPC->GetTPCTrackHisto()->Projection(4,5);
 209   his2D->FitSlicesY(0,0,-1,10,"QNR",&arrayFit);
 210   delete his2D;
 211   his1D = (TH1*) arrayFit.At(1);
 212   his1D->Fit(fpol1,"QNRROB=0.8","QNR",-0.8,-0.1);
 213   offsetdZC=fpol1->GetParameter(0);
 214   slopedZC=fpol1->GetParameter(1);
 215   //offsetdZCchi2=fpol1->GetChisquare();
 216   slopedZCchi2=fpol1->GetChisquare();
 217   //
 218   his1D->Fit(fpol1,"QNRROB=0.8","QNR",0.1,0.8);
 219   offsetdZA=fpol1->GetParameter(0);
 220   slopedZA=fpol1->GetParameter(1);
 221   offsetdZAErr=fpol1->GetParError(0);
 222   slopedZAErr=fpol1->GetParError(1);
 223   offsetdZAchi2=fpol1->GetChisquare();
 224   slopedZAchi2=fpol1->GetChisquare();
 225   //
 226   printf("Drift velocity QA report\n");
 227   printf("offsetdZA\t%f\n",offsetdZA);
 228   printf("slopedZA\t%f\n",slopedZA);
 229   printf("offsetdZC\t%f\n",offsetdZC);
 230   printf("slopedZC\t%f\n",slopedZC);
 231   //
 232   // dump drift QA values
 233   //
 234   (*pcstream)<<"tpcQA"<<
 235     "offsetdZA="<< offsetdZA<<
 236     "slopedZA="<< slopedZA<<
 237     "offsetdZC="<< offsetdZC<<
 238     "slopedZC="<<slopedZC<<
 239     //
 240     "offsetdZAErr="<< offsetdZAErr<<
 241     "slopedZAErr="<< slopedZAErr<<
 242     "offsetdZCErr="<< offsetdZCErr<<
 243     "slopedZCErr="<<slopedZCErr<<
 244     //
 245     "offsetdZAchi2="<< offsetdZAchi2<<
 246     "slopedZAchi2="<< slopedZAchi2<<
 247     "offsetdZCchi2="<< offsetdZCchi2<<
 248     "slopedZCchi2="<<slopedZCchi2;
 249
 250 }
 251
 252
 253 void AnalyzeGain(){
 254   //
 255   // gain statistic
 256   //
 257   AliPerformanceDEdx * pTPCgain =  (AliPerformanceDEdx *)pTPCObjectGain;
 258   static TVectorD MIPvsSector(36);
 259   static TVectorD sector(36);
 260   static Float_t MIPmean = 0;
 261   static Float_t MIPresolution = 0;
 262   static Float_t attachSlopeC = 0;
 263   static Float_t attachSlopeA = 0;
 264
 265   MIPvsSector.Zero();
 266   //
 267   // select MIP particles
 268   //
 269   if (pTPCgain){
 270     pTPCgain->GetDeDxHisto()->GetAxis(7)->SetRangeUser(0.4,0.55);
 271     pTPCgain->GetDeDxHisto()->GetAxis(0)->SetRangeUser(35,60);
 272     pTPCgain->GetDeDxHisto()->GetAxis(6)->SetRangeUser(80,160);
 273     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(-1,1);
 274     //
 275     // MIP position and resolution
 276     //
 277     TF1 gausFit("gausFit","gaus");
 278     TH1 * hisProj1D = pTPCgain->GetDeDxHisto()->Projection(0);
 279     hisProj1D->Fit(&gausFit,"QN","QN");
 280     delete hisProj1D;
 281     MIPmean = gausFit.GetParameter(1);
 282     MIPresolution = 0;
 283     if (MIPmean!=0) MIPresolution = gausFit.GetParameter(2)/MIPmean;
 284     //
 285     // MIP position vs. dip angle (attachment)
 286     //
 287     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(-3,0);
 288     TH2* his2D = pTPCgain->GetDeDxHisto()->Projection(0,5);
 289     TF1 * fpol = new TF1("fpol","pol1");
 290     TObjArray arrayFit;
 291     his2D->FitSlicesY(0,0,-1,10,"QN",&arrayFit);
 292     delete his2D;
 293     TH1 * his1D = (TH1*) arrayFit.At(1);
 294     his1D->Fit(fpol,"QNROB=0.8","QNR",-1,0);
 295     attachSlopeC = fpol->GetParameter(1);
 296     //
 297     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(0,3);
 298     TH2* his2DA = pTPCgain->GetDeDxHisto()->Projection(0,5);
 299     TF1 * fpolA = new TF1("fpolA","pol1");
 300     TObjArray arrayFitA;
 301     his2DA->FitSlicesY(0,0,-1,10,"QN",&arrayFit);
 302     delete his2DA;
 303     TH1 * his1DA = (TH1*) arrayFit.At(1);
 304     his1DA->Fit(fpolA,"QNROB=0.8","QN",0,1);
 305     attachSlopeA = fpolA->GetParameter(1);
 306     //
 307     // MIP position vs. sector
 308     //
 309     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(-3,0); // C side
 310     for(Int_t i = 0; i < 18; i++) { // loop over sectors; correct mapping to be checked!
 311       TH1* his1D=0;
 312       Float_t phiLow = -TMath::Pi() + i*(20./360.)*(2*TMath::Pi());
 313       Float_t phiUp  = -TMath::Pi() + (i+1)*(20./360.)*(2*TMath::Pi());
 314       pTPCgain->GetDeDxHisto()->GetAxis(1)->SetRangeUser(phiLow,phiUp);
 315       his1D = pTPCgain->GetDeDxHisto()->Projection(0);
 316       TF1 gausFunc("gausFunc","gaus");
 317       his1D->Fit(&gausFunc, "QN");
 318       MIPvsSector(i) = gausFunc.GetParameter(1);
 319       sector(i)=i;
 320       delete his1D;
 321     }
 322     //
 323     pTPCgain->GetDeDxHisto()->GetAxis(5)->SetRangeUser(0,3); // A side
 324     for(Int_t i = 0; i < 18; i++) { // loop over sectors; correct mapping to be checked!
 325       TH1* his1D=0;
 326       Float_t phiLow = -TMath::Pi() + i*(20./360.)*(2*TMath::Pi());
 327       Float_t phiUp  = -TMath::Pi() + (i+1)*(20./360.)*(2*TMath::Pi());
 328       pTPCgain->GetDeDxHisto()->GetAxis(1)->SetRangeUser(phiLow,phiUp);
 329       his1D = pTPCgain->GetDeDxHisto()->Projection(0);
 330       TF1 gausFunc("gausFunc","gaus");
 331       his1D->Fit(&gausFunc, "QN");
 332       MIPvsSector(i+18) = gausFunc.GetParameter(1);
 333       sector(i+18)=i+18;
 334       delete his1D;
 335     }
 336   }
 337   //
 338   printf("Gain QA report\n");
 339   printf("MIP mean\t%f\n",MIPmean);
 340   printf("MIP resolution\t%f\n",MIPresolution);
 341   printf("MIPslopeA\t%f\n",attachSlopeA);
 342   printf("MIPslopeC\t%f\n",attachSlopeC);
 343   //
 344   (*pcstream)<<"tpcQA"<<
 345     "MIPattachSlopeC="<<attachSlopeC<<
 346     "MIPattachSlopeA="<<attachSlopeA<<
 347     "MIPresolution="<<MIPresolution<<
 348     "MIPvsSector.="<<&MIPvsSector<<
 349     "sector.="<<&sector<<
 350     "MIPmean="<<MIPmean;
 351
 352 }