PWG/muon/trigEffQA.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 // ROOT includes
   3 #include "TFile.h"
   4 #include "TH1.h"
   5 #include "TH2.h"
   6 #include "TGraphAsymmErrors.h"
   7 #include "TSystem.h"
   8 #include "Riostream.h"
   9 #include "TCanvas.h"
  10 #include "TStyle.h"
  11 #include "TROOT.h"
  12 #include "TLine.h"
  13 #include "TLegend.h"
  14 #include "TPaveLabel.h"
  15 #include "TMath.h"
  16 #include "TObjArray.h"
  17 #include "TObjString.h"
  18 #include "TString.h"
  19 #include "TGrid.h"
  20
  21 #include "AliCDBManager.h"
  22 #include "AliCDBEntry.h"
  23 #include "AliCDBPath.h"
  24 #include "AliMUONTriggerEfficiencyCells.h"
  25 #include "AliMUONTriggerChamberEfficiency.h"
  26 #endif
  27
  28 const Int_t kNch = 4;
  29 const Double_t kZero = 1.e-7; // Avoid problems when comparing to 0.
  30 void SetMyStyle();
  31 Int_t GetRunNumber(TString);
  32 //Double_t GetBinomial(Double_t*, Double_t* eff2 = 0x0, Double_t* effBoth = 0x0);
  33 Double_t* GetBinomial(Double_t*, Double_t* effErr2 = 0x0, Double_t* effErrBoth = 0x0);
  34 Double_t* GetProdErr(Double_t*, Int_t, Int_t nFactors=kNch);
  35 Double_t* GetConditionalEffErr(Double_t*, Double_t*, Double_t*, Int_t exclude=-1);
  36 TH1* GetHisto(TString, TFile*, TList*);
  37
  38 void trigEffQA(TString fileListName, TString outFilename = "")
  39 {
  40   ifstream inFile(fileListName.Data());
  41   TObjArray fileNameArray, runNumArray;
  42   fileNameArray.SetOwner();
  43   runNumArray.SetOwner();
  44   TString currString = "";
  45   if (inFile.is_open()) {
  46     while (! inFile.eof() ) {
  47       currString.ReadLine(inFile,kTRUE); // Read line
  48       if ( currString.IsNull() || ! currString.Contains(".root") ||
  49            currString.BeginsWith("#") ) continue;
  50       fileNameArray.AddLast(new TObjString(currString.Data()));
  51       Int_t runNum = GetRunNumber(currString);
  52       runNumArray.AddLast(new TObjString(Form("%i",runNum)));
  53     }
  54     inFile.close();
  55   }
  56   else {
  57     printf("Fatal: cannot open input file\n");
  58     return;
  59   }
  60
  61   runNumArray.Sort();
  62
  63   //TString elementName[2] = {"Board", "RPC"};
  64   TString countTypeName[3] = {"bendPlane", "nonBendPlane","bothPlanes"};
  65
  66   TObjArray chEffList(13);
  67
  68   TCanvas* can = 0x0;
  69   TString filename = "", histoName = "", histoTitle = "";
  70
  71   histoName = "totalEffEvolution";
  72   histoTitle = "Multinomial probability of firing at least 3/4 chambers";
  73   TH1* totalEff = new TH1D(histoName.Data(), histoTitle.Data(), 0, 0., 1.);
  74   chEffList.AddAtAndExpand(totalEff, 12);
  75
  76   SetMyStyle();
  77   Double_t effValues[3][2*kNch];
  78
  79   TString runNumString = "";
  80   for ( Int_t irun=0; irun<runNumArray.GetEntries(); irun++ ) {
  81     runNumString = runNumArray.At(irun)->GetName();
  82
  83     // Search corresponding file (for sorting)
  84     for ( Int_t ifile=0; ifile<fileNameArray.GetEntries(); ifile++ ) {
  85       filename = fileNameArray.At(ifile)->GetName();
  86       if ( filename.Contains(runNumString.Data()) ) break;
  87     }
  88
  89     if ( filename.Contains("alien://") && ! gGrid )
  90       gGrid->Connect("alien://");
  91
  92     TFile* file = TFile::Open(filename.Data());
  93     if ( ! file ) {
  94       printf("Warning: cannot find %s\n", filename.Data());
  95       continue;
  96     }
  97     TList* trigEffList = (TList*)file->FindObjectAny("triggerChamberEff");
  98     if ( ! trigEffList )
  99       printf("Warning: histo list not found in %s. Check directly in file\n", filename.Data());
 100
 101     TH1* histoDen = GetHisto("allTracksCountChamber", file, trigEffList);
 102     if ( ! histoDen ) {
 103       printf("Error: cannot find histograms in file %s. Skip to next\n", filename.Data());
 104       continue;
 105     }
 106     for(Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++) {
 107       histoName = countTypeName[icount] + "CountChamber";
 108       TH1* histoNum = GetHisto(histoName, file, trigEffList);
 109       TGraphAsymmErrors* graph = new TGraphAsymmErrors(histoNum, histoDen);
 110       Double_t xpt, ypt;
 111       for ( Int_t ich=0; ich<kNch; ich++ ) {
 112         graph->GetPoint(ich, xpt, ypt);
 113         effValues[icount][ich] = ypt;
 114         Int_t ihisto = kNch*icount+ich;
 115         TH1* chEff = (TH1*)chEffList.At(ihisto);
 116         if ( ! chEff ) {
 117           histoName = Form("effEvolution%sCh%i", countTypeName[icount].Data(), 11+ich);
 118           histoTitle = Form("Efficiency of trigger chamber %i", 11+ich);
 119           chEff = new TH1D(histoName.Data(), histoTitle.Data(), 0, 0., 1.);
 120           chEffList.AddAtAndExpand(chEff, ihisto);
 121         }
 122         Int_t currBin = chEff->Fill(runNumString.Data(), ypt);
 123         Double_t err = 0.5*(graph->GetErrorYlow(ich) + graph->GetErrorYhigh(ich));
 124         chEff->SetBinError(currBin, err);
 125         effValues[icount][ich+kNch] = err;
 126       } // loop on chambers
 127     } // loop on counts
 128     Double_t* binomialEff = GetBinomial(effValues[0], effValues[1], effValues[2]);
 129     Int_t currBin = totalEff->Fill(runNumString, binomialEff[0]);
 130     // CAVEAT!!!!
 131     // Well, error calculation of the binomial efficiency is a mess...
 132     // Sometimes it happens that efficiency + error > 1.
 133     // In that case reduce the error.
 134     totalEff->SetBinError(currBin, TMath::Min(binomialEff[1], 1.-binomialEff[0]));
 135     delete binomialEff;
 136   } // loop on runs
 137
 138   TString outCanName = outFilename;
 139   outCanName.ReplaceAll(".root",".pdf");
 140   TString canName = "";
 141   canName = "totalEff";
 142   can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
 143   if ( ! outCanName.IsNull() )
 144     can->Print(Form("%s[", outCanName.Data())); // Open pdf file
 145
 146   // Set correct range (do not show last empty bins)
 147   for ( Int_t ihisto=0; ihisto<chEffList.GetEntries(); ihisto++ ) {
 148     TH1* histo = (TH1*)chEffList.At(ihisto);
 149     histo->GetXaxis()->SetRange(1,fileNameArray.GetEntries());
 150     histo->GetXaxis()->SetLabelSize(0.04);
 151   }
 152
 153   TH1* totEff = (TH1*)chEffList.At(12);
 154   totEff->GetYaxis()->SetRangeUser(0.,1.1);
 155   totEff->GetYaxis()->SetTitle("Probability to satisfy trigger conditions (3/4)");
 156   totEff->SetStats(kFALSE);
 157   totEff->DrawCopy();
 158   if ( ! outCanName.IsNull() )
 159     can->Print(outCanName.Data());
 160
 161   Int_t color[3] = {kBlack, kRed, kBlue};
 162   Int_t markStyle[3] = {20, 24, 26};
 163   TLegend* leg = 0x0;
 164
 165   for ( Int_t ich=0; ich<kNch; ich++ ) {
 166     canName = Form("trigEffCh%i", 11+ich);
 167     can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
 168     can->SetGridy();
 169     leg = new TLegend(0.6, 0.2, 0.9, 0.4);
 170     leg->SetBorderSize(1);
 171     //can->Divide(2,2);
 172     TString drawOpt = "e";
 173     for(Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++) {
 174       //can->cd(icount+1);
 175       TH1* chEff = (TH1*)chEffList.At(4*icount+ich);
 176       chEff->GetYaxis()->SetRangeUser(0.48,1.05);
 177       chEff->SetStats(kFALSE);
 178       chEff->GetYaxis()->SetTitle("Trigger chamber efficiency");
 179       TH1* copyEff = chEff->DrawCopy(drawOpt.Data());
 180       copyEff->SetLineColor(color[icount]);
 181       copyEff->SetMarkerColor(color[icount]);
 182       copyEff->SetMarkerStyle(markStyle[icount]);
 183       leg->AddEntry(copyEff, countTypeName[icount].Data(), "lp");
 184       drawOpt = "esame";
 185     } // loop on counts
 186     leg->Draw("same");
 187     if ( ! outCanName.IsNull() )
 188       can->Print(outCanName.Data());
 189   } // loop on chambers
 190   if ( ! outCanName.IsNull() ) {
 191     can->Print(Form("%s]", outCanName.Data())); // close pdf file
 192     // There is a bug when creating a pdf
 193     // So create a ps and then convert via epstopdf
 194     if ( outCanName.Contains(".ps") || outCanName.Contains(".eps") ) {
 195       gSystem->Exec(Form("epstopdf %s", outCanName.Data()));
 196       gSystem->Exec(Form("rm %s", outCanName.Data()));
 197     }
 198   }
 199
 200   TFile* outFile = 0x0;
 201   if ( ! outFilename.IsNull() ) {
 202     outFile = new TFile(outFilename.Data(), "recreate");
 203     chEffList.Write();
 204     outFile->Close();
 205   }
 206 }
 207
 208
 209 //_____________________________________________________________________________
 210 void SetMyStyle()
 211 {
 212     gStyle->SetCanvasColor(10);
 213     gStyle->SetFrameFillColor(10);
 214     gStyle->SetStatColor(10);
 215     gStyle->SetFillColor(10);
 216     gStyle->SetTitleFillColor(10);
 217
 218     gStyle->SetTitleXSize(0.03);
 219     gStyle->SetTitleXOffset(1.1);
 220     gStyle->SetTitleYSize(0.03);
 221     gStyle->SetTitleYOffset(1.9);
 222
 223     gStyle->SetMarkerSize(0.7);
 224     gStyle->SetHistLineWidth(2);
 225
 226     gStyle->SetPadLeftMargin(0.12);
 227     gStyle->SetPadRightMargin(0.04);
 228     gStyle->SetPadBottomMargin(0.08);
 229     gStyle->SetPadTopMargin(0.08);
 230
 231     gROOT->ForceStyle();
 232 }
 233
 234 //_____________________________________________________________________________
 235 Int_t GetRunNumber(TString filePath)
 236 {
 237   TObjArray* array = filePath.Tokenize("/");
 238   array->SetOwner();
 239   TString auxString = "";
 240   Int_t runNum = -1;
 241   for ( Int_t ientry=0; ientry<array->GetEntries(); ientry++ ) {
 242     auxString = array->At(ientry)->GetName();
 243     if ( auxString.Contains("000") ) {
 244       runNum = auxString.Atoi();
 245       break;
 246     }
 247   }
 248   delete array;
 249
 250   if ( runNum < 0 ) {
 251     array = auxString.Tokenize("_");
 252     array->SetOwner();
 253     auxString = array->Last()->GetName();
 254     auxString.ReplaceAll(".root","");
 255     runNum = auxString.Atoi();
 256     delete array;
 257   }
 258
 259   return runNum;
 260 }
 261
 262
 263 // //_____________________________________________________________________________
 264 // Double_t GetBinomial(Double_t* eff1, Double_t* eff2, Double_t* effBoth)
 265 // {
 266 //   Double_t binomialEff = 0.;
 267 //   Double_t eff44 = 1.;
 268 //   Double_t auxEff[4];
 269 //   Double_t auxBinomial = 1.;
 270
 271 //   // All chamber efficient
 272 //   for ( Int_t ich=0; ich<4; ich++ ) {
 273 //     eff44 *= eff1[ich];
 274 //   }
 275 //   if ( eff2 ) {
 276 //     for ( Int_t ich=0; ich<4; ich++ ) {
 277 //       auxEff[ich] = ( eff1[ich] > kZero ) ? effBoth[ich]/eff1[ich] : 0.;
 278 //     }
 279 //     auxBinomial = GetBinomial(auxEff);
 280 //     eff44 *= auxBinomial;
 281 //   }
 282
 283 //   binomialEff += eff44;
 284
 285 //   // 1 chamber inefficient
 286 //   for ( Int_t ich=0; ich<4; ich++ ) {
 287 //     Double_t eff34 = 1.;
 288 //     for ( Int_t jch=0; jch<4; jch++ ) {
 289 //       eff34 *= ( ich == jch ) ? ( 1. - eff1[jch] ) : eff1[jch];
 290 //     }
 291 //     if ( eff2 ) {
 292 //       for ( Int_t jch=0; jch<4; jch++ ) {
 293 //         if ( ich == jch ) {
 294 //           auxEff[jch] = ( eff1[jch] < 1. ) ? ( eff2[jch] - effBoth[jch] ) / ( 1. - eff1[jch]) : 0.;
 295 //         }
 296 //         else {
 297 //           auxEff[jch] = ( eff1[ich] > kZero ) ? effBoth[ich]/eff1[ich] : 0.;
 298 //         }
 299 //         auxBinomial = GetBinomial(auxEff);
 300 //         eff34 *= auxBinomial;
 301 //       }
 302 //     }
 303
 304 //     binomialEff += eff34;
 305 //   }
 306
 307 //   return binomialEff;
 308 // }
 309
 310
 311 //_____________________________________________________________________________
 312 Double_t* GetBinomial(Double_t* effErr1, Double_t* effErr2, Double_t* effErrBoth)
 313 {
 314   Double_t effProd[4];
 315   Double_t defaultEffErr[2] = {1.,0.};
 316   Double_t* auxBinomial = 0x0;
 317   Double_t* currEffErr44 = 0x0;
 318   Double_t* effErrBinomial = new Double_t[2];
 319   effErrBinomial[0] = 0.;
 320   effErrBinomial[1] = 0.;
 321
 322   for ( Int_t ich = -1; ich<kNch; ich++ ) {
 323     Double_t* currEffErr = GetProdErr(effErr1, ich);
 324     if ( ich >= 0 ) {
 325       currEffErr[0] = currEffErr[0] - currEffErr44[0];
 326       currEffErr[1] = TMath::Sqrt(currEffErr[1]*currEffErr[1] + currEffErr44[1]*currEffErr44[1]);
 327     }
 328     if ( effErr2 ) {
 329       Double_t* auxEffErr = GetConditionalEffErr(effErr1, effErr2, effErrBoth, ich);
 330       auxBinomial = GetBinomial(auxEffErr);
 331       delete auxEffErr;
 332     }
 333     for ( Int_t ival=0; ival<2; ival++ ) {
 334       effProd[2*ival] = currEffErr[ival];
 335       effProd[2*ival+1] = ( effErr2 ) ? auxBinomial[ival] : defaultEffErr[ival];
 336     }
 337     if ( ich < 0 ) currEffErr44 = currEffErr;
 338     else delete currEffErr;
 339     delete auxBinomial;
 340
 341     Double_t* effErr = GetProdErr(effProd, -1, 2);
 342     //printf("%f * %f = %f\n", effProd[0], effProd[1], effErr[0]); // REMEMBER TO CUT
 343     effErrBinomial[0] += effErr[0];
 344     effErrBinomial[1] += effErr[1]*effErr[1];
 345     delete effErr;
 346   } // loop on chambers
 347
 348   delete currEffErr44;
 349
 350   effErrBinomial[1] = TMath::Sqrt(effErrBinomial[1]);
 351
 352   return effErrBinomial;
 353 }
 354
 355 //_____________________________________________________________________________
 356 Double_t* GetProdErr(Double_t* effErr, Int_t exclude, Int_t nFactors)
 357 {
 358   Double_t prod = 1.;
 359   Double_t relErr = 0., relProdErrSquare = 0.;
 360   for ( Int_t iprod=0; iprod<nFactors; iprod++ ) {
 361     if ( iprod == exclude ) continue;
 362     prod *= effErr[iprod];
 363     relErr = ( effErr[iprod] > kZero ) ? effErr[iprod+nFactors]/effErr[iprod] : 0.;
 364     relProdErrSquare += relErr*relErr;
 365     //printf("%f +- %f  ", effErr[iprod], effErr[iprod+nFactors]); // REMEMBER TO CUT
 366   }
 367   Double_t* prodErr = new Double_t[2];
 368   prodErr[0] = prod;
 369   prodErr[1] = prod*TMath::Sqrt(relProdErrSquare);
 370   //printf("-> %f %f\n", prodErr[0], prodErr[1]); // REMEMBER TO CUT
 371   return prodErr;
 372 }
 373
 374
 375 //_____________________________________________________________________________
 376 Double_t* GetConditionalEffErr(Double_t* effErr1, Double_t* effErr2, Double_t* effErrBoth, Int_t exclude)
 377 {
 378   Double_t* effErr = new Double_t[2*kNch];
 379   for ( Int_t ich=0; ich<kNch; ich++ ) {
 380     if ( ich == exclude ) {
 381       effErr[ich] = ( effErr1[ich] < 1. ) ? ( effErr2[ich] - effErrBoth[ich] ) / ( 1. - effErr1[ich] ) : 0.;
 382       effErr[ich+kNch] = 0;
 383       if ( effErr1[ich] < 1. ) {
 384         Double_t err2 = effErr2[ich+kNch] / ( 1. - effErr1[ich] );
 385         Double_t errBoth = effErrBoth[ich+kNch] / ( 1. - effErr1[ich] );
 386         Double_t err1 = effErr1[ich+kNch] * effErr[ich] / ( 1. - effErr1[ich] );
 387         effErr[ich+kNch] = TMath::Sqrt(err2*err2 + errBoth*errBoth + err1*err1);
 388       }
 389     }
 390     else {
 391       effErr[ich] = ( effErr1[ich] > kZero ) ? effErrBoth[ich]/effErr1[ich] : 0.;
 392       Double_t relErr1 = ( effErr1[ich] > kZero ) ? effErr1[ich+kNch]/effErr1[ich] : 0.;
 393       Double_t relErrBoth = ( effErrBoth[ich] > kZero ) ? effErrBoth[ich+kNch]/effErrBoth[ich] : 0.;
 394       effErr[ich+kNch] = effErr[ich] * TMath::Sqrt(relErr1*relErr1 + relErrBoth*relErrBoth);
 395     }
 396     //printf("%f  %f  %f -> %f\n", effErr1[ich], effErr2[ich], effErrBoth[ich], effErr[ich]); // REMEMBER TO CUT
 397   } // loop on chambers
 398   return effErr;
 399 }
 400
 401
 402 //_____________________________________________________________________________
 403 TH1* GetHisto(TString histoName, TFile* file, TList* histoList)
 404 {
 405   TH1* histo = 0x0;
 406   if ( histoList )
 407     histo = (TH1*)histoList->FindObject(histoName.Data());
 408   else
 409     histo = (TH1*)file->FindObjectAny(histoName.Data());
 410
 411   return histo;
 412 }