Add centrality dependence in pA for muon tracking QA task and T0 trigger to normalize...
[u/mrichter/AliRoot.git] / PWGPP / MUON / lite / trigEffQA.C
1
2 #if !defined(__CINT__) || defined(__MAKECINT__)
3 // ROOT includes
4 #include "TFile.h"
5 #include "TH1.h"
6 #include "TH2.h"
7 #include "TGraphAsymmErrors.h"
8 #include "TSystem.h"
9 #include "Riostream.h"
10 #include "TCanvas.h"
11 #include "TStyle.h"
12 #include "TROOT.h"
13 #include "TLegend.h"
14 #include "TMath.h"
15 #include "TObjArray.h"
16 #include "TList.h"
17 #include "TObjString.h"
18 #include "TString.h"
19 #include "TGrid.h"
20 #include "TArrayD.h"
21 #include "TArrayI.h"
22 #include "TMap.h"
23 #include "TGridResult.h"
24
25 #include "AliCDBManager.h"
26 #include "AliCDBEntry.h"
27 #include "AliCDBPath.h"
28 #include "AliCDBStorage.h"
29 #include "AliMUONTriggerEfficiencyCells.h"
30 #include "AliMUONTriggerChamberEfficiency.h"
31 #include "AliMUONTriggerUtilities.h"
32 #include "AliMUONDigitMaker.h"
33 #include "AliMUONVDigit.h"
34 #include "AliMUONDigitStoreV2R.h"
35 #include "AliMUONCalibrationData.h"
36 #include "AliAnalysisTriggerScalers.h"
37 #include "AliCounterCollection.h"
38 #endif
39
40 const Int_t kNch = 4;
41 const Double_t kZero = 1.e-7; // Avoid problems when comparing to 0.
42
43 const Int_t kNCollType = 4; 
44 TString kSCollType[kNCollType] = {"-B-","-S-","-SC-","-SA-"}; 
45
46
47 //_____________________________________________________________________________
48 void SetMyStyle()
49 {
50   gStyle->SetCanvasColor(10);
51   gStyle->SetFrameFillColor(10);
52   gStyle->SetStatColor(10);
53   gStyle->SetFillColor(10);
54   gStyle->SetTitleFillColor(10);
55   
56   gStyle->SetTitleXSize(0.03);
57   gStyle->SetTitleXOffset(1.1);
58   gStyle->SetTitleYSize(0.03);
59   gStyle->SetTitleYOffset(1.9);
60   
61   gStyle->SetMarkerSize(0.7);
62   gStyle->SetHistLineWidth(2);
63   
64   gStyle->SetPadLeftMargin(0.12);
65   gStyle->SetPadRightMargin(0.04);
66   gStyle->SetPadBottomMargin(0.08);
67   gStyle->SetPadTopMargin(0.08);
68   
69   gROOT->ForceStyle();
70 }
71
72
73 //_____________________________________________________________________________
74 void SetRunAxisRange ( TAxis* axis )
75 {
76   for ( Int_t ibin=1; ibin<=axis->GetNbins(); ibin++ ) {
77     TString binLabel = axis->GetBinLabel(ibin);
78     if ( ! binLabel.IsNull()) continue;
79     axis->SetRange(1, ibin-1);
80     return;
81   }
82 }
83
84 //_____________________________________________________________________________
85 Int_t GetRunNumber(TString filePath)
86 {
87   TObjArray* array = filePath.Tokenize("/");
88   array->SetOwner();
89   TString auxString = "";
90   Int_t runNum = -1;
91   for ( Int_t ientry=0; ientry<array->GetEntries(); ientry++ ) {
92     auxString = array->At(ientry)->GetName();
93     if ( auxString.BeginsWith("000") ) {
94       runNum = auxString.Atoi();
95       break;
96     }
97   }
98   delete array;
99   
100   if ( runNum < 0 ) {
101     array = auxString.Tokenize("_");
102     array->SetOwner();
103     auxString = array->Last()->GetName();
104     auxString.ReplaceAll(".root","");
105     runNum = auxString.Atoi();
106     delete array;
107   }
108   
109   return runNum;
110 }
111
112 //_____________________________________________________________________________
113 Double_t* GetProdErr(Double_t* effErr, Int_t exclude, Int_t nFactors = kNch)
114 {
115   Double_t prod = 1.;
116   Double_t relErr = 0., relProdErrSquare = 0.;
117   for ( Int_t iprod=0; iprod<nFactors; iprod++ ) {
118     if ( iprod == exclude ) continue;
119     prod *= effErr[iprod];
120     relErr = ( effErr[iprod] > kZero ) ? effErr[iprod+nFactors]/effErr[iprod] : 0.;
121     relProdErrSquare += relErr*relErr;
122     //printf("%f +- %f  ", effErr[iprod], effErr[iprod+nFactors]); // alBER TO CUT
123   }
124   Double_t* prodErr = new Double_t[2];
125   prodErr[0] = prod;
126   prodErr[1] = prod*TMath::Sqrt(relProdErrSquare);
127   //printf("-> %f %f\n", prodErr[0], prodErr[1]); // REMEMBER TO CUT
128   return prodErr;
129 }
130
131
132 //_____________________________________________________________________________
133 Double_t* GetConditionalEffErr(Double_t* effErr1, Double_t* effErr2, Double_t* effErrBoth, Int_t exclude = -1)
134 {
135   Double_t* effErr = new Double_t[2*kNch];
136   for ( Int_t ich=0; ich<kNch; ich++ ) {
137     if ( ich == exclude ) {
138       effErr[ich] = ( effErr1[ich] < 1. ) ? ( effErr2[ich] - effErrBoth[ich] ) / ( 1. - effErr1[ich] ) : 0.;
139       effErr[ich+kNch] = 0;
140       if ( effErr1[ich] < 1. ) {
141         Double_t err2 = effErr2[ich+kNch] / ( 1. - effErr1[ich] );
142         Double_t errBoth = effErrBoth[ich+kNch] / ( 1. - effErr1[ich] );
143         Double_t err1 = effErr1[ich+kNch] * effErr[ich] / ( 1. - effErr1[ich] );
144         effErr[ich+kNch] = TMath::Sqrt(err2*err2 + errBoth*errBoth + err1*err1);
145       }
146     }
147     else {
148       effErr[ich] = ( effErr1[ich] > kZero ) ? effErrBoth[ich]/effErr1[ich] : 0.;
149       Double_t relErr1 = ( effErr1[ich] > kZero ) ? effErr1[ich+kNch]/effErr1[ich] : 0.;
150       Double_t relErrBoth = ( effErrBoth[ich] > kZero ) ? effErrBoth[ich+kNch]/effErrBoth[ich] : 0.;
151       effErr[ich+kNch] = effErr[ich] * TMath::Sqrt(relErr1*relErr1 + relErrBoth*relErrBoth);
152     }
153     //printf("%f  %f  %f -> %f\n", effErr1[ich], effErr2[ich], effErrBoth[ich], effErr[ich]); // REMEMBER TO CUT
154   } // loop on chambers
155   return effErr;
156 }
157
158
159 //_____________________________________________________________________________
160 Double_t* GetBinomial(Double_t* effErr1, Double_t* effErr2 = 0x0, Double_t* effErrBoth = 0x0)
161 {
162   Double_t effProd[4];
163   Double_t defaultEffErr[2] = {1.,0.};
164   Double_t* auxBinomial = 0x0;
165   Double_t* currEffErr44 = 0x0;
166   Double_t* effErrBinomial = new Double_t[2];
167   effErrBinomial[0] = 0.;
168   effErrBinomial[1] = 0.;
169   
170   for ( Int_t ich = -1; ich<kNch; ich++ ) {
171     Double_t* currEffErr = GetProdErr(effErr1, ich);
172     if ( ich >= 0 ) {
173       currEffErr[0] = currEffErr[0] - currEffErr44[0];
174       currEffErr[1] = TMath::Sqrt(currEffErr[1]*currEffErr[1] + currEffErr44[1]*currEffErr44[1]);
175     }
176     if ( effErr2 ) {
177       Double_t* auxEffErr = GetConditionalEffErr(effErr1, effErr2, effErrBoth, ich);
178       auxBinomial = GetBinomial(auxEffErr);
179       delete auxEffErr;
180     }
181     for ( Int_t ival=0; ival<2; ival++ ) {
182       effProd[2*ival] = currEffErr[ival];
183       effProd[2*ival+1] = ( effErr2 ) ? auxBinomial[ival] : defaultEffErr[ival];
184     }
185     if ( ich < 0 ) currEffErr44 = currEffErr;
186     else delete currEffErr;
187     delete auxBinomial;
188     
189     Double_t* effErr = GetProdErr(effProd, -1, 2);
190     //printf("%f * %f = %f\n", effProd[0], effProd[1], effErr[0]); // REMEMBER TO CUT
191     effErrBinomial[0] += effErr[0];
192     effErrBinomial[1] += effErr[1]*effErr[1];
193     delete effErr;
194   } // loop on chambers
195   
196   delete currEffErr44;
197   
198   effErrBinomial[1] = TMath::Sqrt(effErrBinomial[1]);
199   
200   return effErrBinomial;
201 }
202
203
204 //_____________________________________________________________________________
205 TH1* GetHisto(TString histoName, TFile* file, TList* histoList)
206 {
207   TH1* histo = 0x0;
208   if ( histoList )
209     histo = (TH1*)histoList->FindObject(histoName.Data());
210   else
211     histo = (TH1*)file->FindObjectAny(histoName.Data());
212   
213   return histo;
214 }
215
216 //_____________________________________________________________________________
217 Int_t GetEffIndex ( Int_t iel, Int_t icount, Int_t ich = -1 )
218 {
219   if ( iel == 0 ) return icount;
220   return 3 + 4*3*(iel-1) + 3*ich + icount;
221 }
222
223 //_____________________________________________________________________________
224 TList* GetOCDBList ( )
225 {
226   TString baseFolder = AliCDBManager::Instance()->GetDefaultStorage()->GetBaseFolder();
227   
228   TList* outList = new TList();
229   outList->SetOwner();
230   TString dirName[3] = {"GlobalTriggerCrateConfig", "RegionalTriggerConfig", "LocalTriggerBoardMasks"};
231   for ( Int_t idir=0; idir<3; idir++ ) {
232     TGridResult *res = gGrid->Ls(Form("%s/MUON/Calib/%s", baseFolder.Data(), dirName[idir].Data()));
233     if (!res) return 0x0;
234     for ( Int_t ires=0; ires<res->GetEntries(); ires++ ) {
235       TString currFile = static_cast<TMap*>(res->At(ires))->GetValue("name")->GetName();
236       outList->Add(new TObjString(currFile));
237     }
238   }
239   return outList;
240 }
241
242 //_____________________________________________________________________________
243 Bool_t IsOCDBChanged ( Int_t currRun, Int_t previousRun, TList* fileList )
244 {
245   if ( ! fileList ) return kTRUE;
246   for ( Int_t ifile=0; ifile<fileList->GetEntries(); ifile++ ) {
247     TString filename = static_cast<TObjString*>(fileList->At(ifile))->GetString();
248     filename.ReplaceAll("Run","");
249     TObjArray* array = filename.Tokenize("_");
250     Int_t firstRun = static_cast<TObjString*>(array->At(0))->GetString().Atoi();
251     Int_t lastRun = static_cast<TObjString*>(array->At(1))->GetString().Atoi();
252     delete array;
253     Bool_t isCurrRunInside = ( currRun >= firstRun && currRun <= lastRun );
254     Bool_t isPreviousRunInside = ( previousRun >= firstRun && previousRun <= lastRun );
255     if ( isCurrRunInside != isPreviousRunInside ) return kTRUE;
256   }
257   return kFALSE;
258 }
259
260 //_____________________________________________________________________________
261 void TrigEffTrending(TObjArray runNumArray, TObjArray fileNameArray, TList& outCanList, TList& outList)
262 {
263   TString elementName[3] = { "Chamber", "RPC", "Board" };
264   TString countTypeName[4] = { "allTracks", "bendPlane", "nonBendPlane", "bothPlanes" };
265
266   TString filename = "", effName = "", effTitle = "";
267
268   SetMyStyle();
269   Double_t effValues[3][2*kNch];
270   const Int_t kNgraphs = kNch*3*2+3;
271   TObjArray effList(kNgraphs);
272   const Int_t kNeffVsRun = kNgraphs+1;
273   TObjArray effVsRunList(kNeffVsRun);
274   
275   effName = "totalEffEvolution";
276   effTitle = "Multinomial probability of firing at least 3/4 chambers";
277   TH1D* totalEff = new TH1D(effName.Data(), effTitle.Data(), 0, 0., 1.);
278   effVsRunList.AddAt(totalEff, kNeffVsRun-1);
279
280   TString runNumString = "";
281   for ( Int_t irun=0; irun<runNumArray.GetEntries(); irun++ ) {
282     runNumString = runNumArray.At(irun)->GetName();
283
284     // Search corresponding file (for sorting)
285     for ( Int_t ifile=0; ifile<fileNameArray.GetEntries(); ifile++ ) {
286       filename = fileNameArray.At(ifile)->GetName();
287       if ( filename.Contains(runNumString.Data()) ) break;
288     }
289
290     if ( filename.Contains("alien://") && ! gGrid ) gGrid->Connect("alien://");
291     
292     //
293     // First get the list of efficiency graphs
294     //
295     
296     // Chamber efficiency
297     TFile* file = TFile::Open(filename.Data());
298     if ( ! file ) {
299       printf("Warning: cannot find %s\n", filename.Data());
300       continue;
301     }
302     
303     TList* trigEffList = (TList*)file->FindObjectAny("triggerChamberEff");
304     if ( ! trigEffList ) printf("Warning: histo list not found in %s. Check directly in file\n", filename.Data());
305     
306     TH1* histoDen = 0x0;
307     for ( Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts; icount++ ) {
308       effName = countTypeName[icount] + "Count" + elementName[0];
309       if ( icount == 0 ) {
310         histoDen = GetHisto(effName, file, trigEffList);
311         continue;
312       }
313       
314       TH1* histoNum = GetHisto(effName, file, trigEffList);
315       TGraphAsymmErrors* graph = new TGraphAsymmErrors(histoNum, histoDen);
316       effName.ReplaceAll("Count","Eff");
317       graph->SetName(effName.Data());
318       effList.AddAt(graph, GetEffIndex(0, icount-1));
319     }
320     file->Close();
321   
322     if ( ! histoDen ) {
323       printf("Error: cannot find histograms in file %s. Skip to next\n", filename.Data());
324       continue;
325     }
326   
327     // RPC/board efficiency
328     AliMUONTriggerChamberEfficiency trigChEff(filename);
329     for ( Int_t iel=1; iel<3; iel++ ) {
330       for ( Int_t ich=0; ich<kNch; ich++ ) {
331         for ( Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++ ) {
332           TObject* obj = trigChEff.GetEffObject(2-iel, icount, ich);
333           effList.AddAt(obj, GetEffIndex(iel, icount, ich));
334         }
335       }
336     }
337     
338     // Fill efficiency vs run
339     for ( Int_t iel=0; iel<3; iel++ ) {
340       for ( Int_t ich=0; ich<kNch; ich++ ) {
341         for ( Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++ ) {
342           TGraphAsymmErrors* graph = static_cast<TGraphAsymmErrors*>(effList.At(GetEffIndex(iel, icount, ich)));
343           Int_t nPoints = ( iel == 0 ) ? 1 : graph->GetN();
344           for ( Int_t ipoint=0; ipoint<nPoints; ipoint++ ) {
345             Int_t currPoint = ( iel == 0 ) ? ich : ipoint;
346             Double_t xpt, ypt;
347             graph->GetPoint(currPoint, xpt, ypt);
348             effValues[icount][ich] = ypt;
349             Int_t ihisto = GetEffIndex(iel,icount,ich);
350             TH2* effHisto = static_cast<TH2*>(effVsRunList.At(ihisto));
351             if ( ! effHisto ) {
352               effName = Form("effEvolution%s%s", countTypeName[icount+1].Data(), elementName[iel].Data());
353               effTitle = Form("Trigger chamber efficiency vs run");
354               if ( iel>0 ) {
355                 effName += Form("Ch%i", 11+ich);
356                 effTitle += Form(" for chamber %i", 11+ich);
357               }
358               effHisto = new TH2D(effName.Data(), effTitle.Data(), 0, 0., 1., graph->GetN(), xpt-0.5, xpt-0.5+(Double_t)graph->GetN());
359               effVsRunList.AddAt(effHisto, ihisto);
360             }
361             Int_t currBin = effHisto->Fill(runNumString.Data(), xpt, ypt);
362             Double_t err = 0.5*(graph->GetErrorYlow(ipoint) + graph->GetErrorYhigh(ipoint));
363             Int_t binx, biny, binz;
364             effHisto->GetBinXYZ(currBin, binx, biny, binz);
365             effHisto->SetBinError(binx, biny, err);
366             effValues[icount][ich+kNch] = err;
367           } // loop on points
368         } // loop on counts
369       } // loop on chambers
370       if ( iel > 0 ) continue;
371       Double_t* binomialEff = GetBinomial(effValues[0], effValues[1], effValues[2]);
372       Int_t currBin = totalEff->Fill(runNumString, binomialEff[0]);
373       // CAVEAT!!!!
374       // Well, error calculation of the binomial efficiency is a mess...
375       // Sometimes it happens that efficiency + error > 1.
376       // In that case reduce the error.
377       totalEff->SetBinError(currBin, TMath::Min(binomialEff[1], 1.-binomialEff[0]));
378       delete binomialEff;
379     } // loop on detection elements
380   } // loop on runs
381   
382   // Set correct range (do not show last empty bins)
383   for ( Int_t ihisto=0; ihisto<effVsRunList.GetEntries(); ihisto++ ) {
384     TH1* histo = static_cast<TH1*>(effVsRunList.At(ihisto));
385     SetRunAxisRange(histo->GetXaxis());
386     outList.Add(histo);
387     //histo->GetXaxis()->SetLabelSize(0.03);
388   }
389
390   TString canName = "totalEff";
391   TCanvas* can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
392   TH1* totEff = (TH1*)effVsRunList.At(kNeffVsRun-1);
393   totEff->GetYaxis()->SetRangeUser(0.9,1.05);
394   totEff->GetYaxis()->SetTitle("Probability to satisfy trigger conditions (3/4)");
395   totEff->SetStats(kFALSE);
396   totEff->DrawCopy();
397   outCanList.Add(can);
398
399   Int_t color[3] = {kBlack, kRed, kBlue};
400   Int_t markStyle[3] = {20, 24, 26};
401   TLegend* leg = 0x0;
402
403   for ( Int_t ich=0; ich<kNch; ich++ ) {
404     canName = Form("trigEffCh%i", 11+ich);
405     can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
406     can->SetGridy();
407     leg = new TLegend(0.6, 0.2, 0.9, 0.4);
408     leg->SetBorderSize(1);
409     //can->Divide(2,2);
410     TString drawOpt = "e";
411     for(Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++) {
412       //can->cd(icount+1);
413       TH2* histo = static_cast<TH2*>(effVsRunList.At(GetEffIndex(0, icount)));
414       TH1* chEff = histo->ProjectionX(Form("effEvolutionCh%i",11+ich), ich+1, ich+1);
415       chEff->SetTitle(Form("%s for chamber %i", histo->GetTitle(), 11+ich));
416       chEff->GetYaxis()->SetRangeUser(0.9,1.);
417       chEff->SetStats(kFALSE);
418       chEff->GetYaxis()->SetTitle("Trigger chamber efficiency");
419       TH1* copyEff = chEff->DrawCopy(drawOpt.Data());
420       copyEff->SetLineColor(color[icount]);
421       copyEff->SetMarkerColor(color[icount]);
422       copyEff->SetMarkerStyle(markStyle[icount]);
423       leg->AddEntry(copyEff, countTypeName[icount+1].Data(), "lp");
424       drawOpt = "esame";
425     } // loop on counts
426     leg->Draw("same");
427     outCanList.Add(can);
428   } // loop on chambers
429   
430   for ( Int_t iel=1; iel<3; iel++ ) {
431     for ( Int_t ich=0; ich<kNch; ich++ ) {
432       Int_t icount = AliMUONTriggerEfficiencyCells::kBothPlanesEff; // Just plot the efficiency for both
433 //      for ( Int_t icount=0; icount<AliMUONTriggerEfficiencyCells::kNcounts-1; icount++ ) {
434         canName = Form("trigEff%sCh%i", elementName[iel].Data(), 11+ich);
435         can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
436         can->SetRightMargin(0.14);
437         TH2* histo = static_cast<TH2*>(effVsRunList.At(GetEffIndex(iel, icount,ich)));
438         histo->SetStats(kFALSE);
439         histo->GetYaxis()->SetTitle(elementName[iel].Data());
440         histo->DrawCopy("COLZ");
441 //      } // loop on counts
442       outCanList.Add(can);
443     } // loop on chambers
444   } // loop on detection element type
445 }
446
447 //_____________________________________________________________________________
448 void MaskTrending ( TObjArray runNumArray, TString defaultStorage, TList& outCanList, TList& outList )
449 {
450  
451   if ( defaultStorage.Contains("alien://") || defaultStorage.Contains("raw://") ) {
452     if ( ! gGrid ) TGrid::Connect("alien://");
453     if ( ! gGrid ) {
454       printf("Error: Problem connetting to grid: nothing done");
455       return;
456     }
457   }
458   
459   
460   TObjArray maskedList(8);
461   TObjArray auxList(8);
462   auxList.SetOwner();
463   TString histoName = "", histoTitle = "";
464   for(Int_t icath=0; icath<2; icath++){
465     TString cathName = ( icath==0 ) ? "bendPlane" : "nonBendPlane";
466     for(Int_t ich=0; ich<kNch; ich++){
467       histoName = Form("%sMaskCh%i", cathName.Data(), 11+ich);
468       histoTitle = Form("Chamber %i - %s: fraction of masked channels", 11+ich, cathName.Data());
469       TH2* histo = new TH2D(histoName.Data(), histoTitle.Data(),0,0.,1., 234, 0.5, 234. + 0.5);
470       histo->GetYaxis()->SetTitle("Board Id");
471       histo->SetOption("COLZ");
472       Int_t imask = 2*ich + icath;
473       maskedList.AddAt(histo, imask);
474       auxList.AddAt(histo->Clone(Form("%s_aux",histoName.Data())), imask);
475     } // loop on chambers
476   } // loop on cathodes
477   
478   TArrayS xyPatternAll[2];
479   for(Int_t icath=0; icath<2; icath++){
480     xyPatternAll[icath].Set(kNch);
481     xyPatternAll[icath].Reset(0xFFFF);
482   }
483   
484   AliCDBManager::Instance()->SetDefaultStorage(defaultStorage.Data());
485   TList* ocdbFileList = 0x0;
486   Int_t previousRun = -1;
487   AliMUONDigitMaker* digitMaker = 0x0;
488   AliMUONDigitStoreV2R digitStore;
489   
490   AliMUONCalibrationData* calibData = 0x0;
491   AliMUONTriggerUtilities* trigUtilities = 0x0;
492   for ( Int_t irun=0; irun<runNumArray.GetEntries(); irun++ ) {
493     TString runNumString = runNumArray.At(irun)->GetName();
494     Int_t runNumber = runNumString.Atoi();
495     
496     if ( IsOCDBChanged(runNumber, previousRun, ocdbFileList) ) {
497       AliCDBManager::Instance()->SetRun(runNumber);
498       
499       if ( ! digitMaker ) {
500         digitMaker = new AliMUONDigitMaker(kFALSE);
501         // Create a store with all digits in trigger
502         for ( Int_t iboard=1; iboard<=234; iboard++ ) {
503           digitMaker->TriggerDigits(iboard, xyPatternAll, digitStore, kFALSE);
504         }
505       }
506       
507       if ( ! ocdbFileList ) ocdbFileList = GetOCDBList();
508   
509       delete calibData;
510       calibData = new AliMUONCalibrationData (runNumber);
511       delete trigUtilities;
512       trigUtilities = new AliMUONTriggerUtilities (calibData);
513     }
514     
515     previousRun = runNumber;
516     
517     TIter next(digitStore.CreateIterator());
518     AliMUONVDigit* dig = 0x0;
519     while ( ( dig = static_cast<AliMUONVDigit*>(next()) ) ) {
520       Int_t icath = dig->Cathode();
521       Int_t detElemId = dig->DetElemId();
522       Int_t ich = detElemId/100-11;
523       Int_t iboard = dig->ManuId();
524       Int_t imask = 2*ich + icath;
525       static_cast<TH2*>(auxList.At(imask))->Fill(runNumString.Data(),iboard,1.);
526       static_cast<TH2*>(maskedList.At(imask))->Fill(runNumString.Data(),iboard,(Double_t)trigUtilities->IsMasked(*dig));
527     }
528   } // loop on runs
529   delete calibData;
530   delete trigUtilities;
531   delete digitMaker;
532
533   TString canName = "";
534   for ( Int_t imask=0; imask<maskedList.GetEntries(); imask++ ) {
535     TH2* histo = static_cast<TH2*>(maskedList.At(imask));
536     histo->Divide(static_cast<TH2*>(auxList.At(imask)));
537     SetRunAxisRange(histo->GetXaxis());
538     outList.Add(histo);
539     
540     canName = Form("%sCan", histo->GetName());
541     TCanvas* can = new TCanvas(canName.Data(), canName.Data(), 200, 10, 600, 600);
542     can->SetRightMargin(0.14);
543     histo->SetStats(kFALSE);
544     histo->DrawCopy("COLZ");
545     outCanList.Add(can);
546   }
547 }
548 //_____________________________________________________________________________
549 TObjArray* BuildListOfT0Trigger(TObjArray* selectedTriggerArray, TString sT0Name = "C0TVX") {
550
551   if ( !selectedTriggerArray ) return 0;
552   
553   TObjArray* selectedList = new TObjArray();
554   selectedList->SetOwner();
555   
556   TString sT0TriggerName = "";
557   TString currTrigName = "";
558   Bool_t isCollType[kNCollType];
559
560   for ( Int_t iCol = 0; iCol < kNCollType; iCol++ ) isCollType[iCol] = kFALSE;
561
562   for (Int_t iTrig = 0; iTrig < selectedTriggerArray->GetEntries(); iTrig++){
563     currTrigName = ((TObjString*)selectedTriggerArray->At(iTrig))->GetName();
564     //search the collision type
565     for ( Int_t iCol = 0; iCol < kNCollType; iCol++ ) {
566       if ( currTrigName.Contains( kSCollType[iCol]  ) ) isCollType[iCol] = kTRUE;
567     }
568   }
569   for ( Int_t iCol = 0; iCol < kNCollType; iCol++ ) {
570     if ( !isCollType[iCol] ) continue;
571     sT0TriggerName = Form("%s%sNOPF-ALLNOTRD",sT0Name.Data(),kSCollType[iCol].Data());
572     selectedList->AddLast(new TObjString(sT0TriggerName.Data()));
573   }
574
575   return selectedList;
576 }
577
578 //_____________________________________________________________________________
579 TObjArray* BuildListOfTrigger(TString triggerListName, TString rejectPattern="OTHER TRUE PHI ANY EMC -ACE- -ABCE- WU MUP SPI SHM") {
580
581   TObjArray* selectedList = new TObjArray();
582   selectedList->SetOwner();
583   TObjArray* triggerArray = triggerListName.Tokenize(",");
584   TObjArray* rejectArray = rejectPattern.Tokenize(" ");
585   TString currTrigName = "";
586
587   for ( Int_t iTrig = 0; iTrig < triggerArray->GetEntries(); iTrig++ ){
588     currTrigName = ((TObjString*)triggerArray->At(iTrig))->GetName();
589     Bool_t isGoodTrig = kTRUE;
590     for ( Int_t irej=0; irej<rejectArray->GetEntries(); irej++ ) {
591       if (currTrigName.Contains(rejectArray->At(irej)->GetName()) ) {
592         isGoodTrig = kFALSE;
593         break;
594       }
595     }
596     if ( isGoodTrig) selectedList->AddLast(new TObjString(currTrigName.Data()));
597   }
598
599   delete triggerArray;
600   delete rejectArray;
601
602   return selectedList;
603
604 }
605
606 //_____________________________________________________________________________
607 TString FindT0TriggerName(TObjArray *T0TrigArray,TString sTrig, TString sT0Name = "C0TVX") {
608
609   TString name = "";
610   if ( !T0TrigArray ) return name;
611   TString currTrigName = "";
612
613   //search the collision type
614   for ( Int_t iCol = 0; iCol < kNCollType; iCol++ ) {
615     if ( sTrig.Contains( kSCollType[iCol] ) ) {
616       for ( Int_t iTrig = 0; iTrig < T0TrigArray->GetEntries(); iTrig++ ){
617         currTrigName = ((TObjString*) T0TrigArray->At(iTrig))->GetName();
618         if ( currTrigName.Contains( kSCollType[iCol] )  ) return currTrigName;       
619       }
620     }
621     
622   }
623   return name;
624 }
625   
626 //_____________________________________________________________________________
627 void ScalerTrending ( TObjArray runNumArray, TString mergedFileName, TString defaultStorage, TList& outCanList, TList& outList )
628 {
629
630   if ( defaultStorage.Contains("alien://") || defaultStorage.Contains("raw://") ) {
631     if ( ! gGrid ) TGrid::Connect("alien://");
632     if ( ! gGrid ) {
633       printf("Error: Problem connetting to grid: nothing done");
634       return;
635     }
636   }
637
638   //trigger count from ESDs
639   TFile *file = TFile::Open(mergedFileName.Data());
640   AliCounterCollection* ccol = (AliCounterCollection*)((TDirectoryFile*)file->FindObjectAny("MUON_QA"))->FindObjectAny("eventCounters");
641   
642   //Build the trigger list for trigger with muon only in readout and min. bias triggers
643   TString triggerListName = ccol->GetKeyWords("trigger");
644   TObjArray* selectedTriggerArray = BuildListOfTrigger(triggerListName);
645   if ( !selectedTriggerArray ) {
646     printf("No trigger selected from trigger list %s\n",triggerListName.Data());
647     return;
648   }
649   printf("Nr of triggers selected %i\n",selectedTriggerArray->GetEntries());
650
651   TObjArray* selectedT0TriggerArray = BuildListOfT0Trigger(selectedTriggerArray);
652   printf("Nr of T0 triggers selected %i\n",selectedT0TriggerArray->GetEntries());
653   
654   const Int_t nScaler = 3;
655   TString sScaler[nScaler] = {"L0B","L2A","L0BRATE"};
656   enum eScaler {kL0B = 0, kL2A=1, kL0BRATE=2};
657   Float_t maxScaler[nScaler] = {1e8,1e7,1e6};
658   TObjArray hFromQA;
659   TObjArray hFromScalers; 
660   TObjArray hOutput;
661
662   TString sHistName, sHistNameFull, sTitleName;
663   Int_t nRuns = runNumArray.GetEntries();
664   //
665   //Create histos for Scalers and QA 
666   //
667
668   //loop on recorded trigger list
669   for ( Int_t iTrig = 0; iTrig < selectedTriggerArray->GetEntries(); iTrig++ ) {
670     //loop on scaler list
671     for ( Int_t iScaler = 0; iScaler < nScaler; iScaler++ ) {
672       sHistName = ((TObjString*) selectedTriggerArray->At(iTrig))->GetName();
673       sHistName += "_";
674       sHistName += sScaler[iScaler];
675       sTitleName = sHistName;
676
677       //      printf("Trigger: %s\n",((TObjString*) selectedTriggerArray->At(iTrig))->GetName());
678            
679       sHistNameFull = "Scalers_";
680       sHistNameFull += sHistName;
681       TH1F* hCounterScalers = new TH1F(sHistNameFull,sTitleName,nRuns,1,(Double_t)nRuns);
682       hCounterScalers->SetDirectory(0);
683       hCounterScalers->SetMinimum(1);
684       hCounterScalers->SetMaximum(maxScaler[iScaler]);
685       hFromScalers.AddLast(hCounterScalers);
686       hOutput.AddLast(hCounterScalers);
687       
688       if ( !(sScaler[iScaler].Contains("L2A")) ) continue;
689
690       sHistNameFull = "QA_";
691       sHistNameFull += sHistName;
692       TH1F* hCounterQA = new TH1F(sHistNameFull,sTitleName,nRuns,1,(Double_t)nRuns);
693       hCounterQA->SetDirectory(0);
694       hFromQA.AddLast(hCounterQA);
695       
696     }
697   }
698   //loop on T0 trigger list (only for L0B)
699   for ( Int_t iTrig = 0; iTrig < selectedT0TriggerArray->GetEntries(); iTrig++ ) {
700
701     sHistName = ((TObjString*) selectedT0TriggerArray->At(iTrig))->GetName();
702     sHistName += "_";
703     sHistName += sScaler[0];
704     sTitleName = sHistName;
705
706     printf("Trigger: %s\n",((TObjString*) selectedT0TriggerArray->At(iTrig))->GetName());
707     
708     sHistNameFull = "Scalers_";
709     sHistNameFull += sHistName;
710     TH1F* hCounterScalers = new TH1F(sHistNameFull,sTitleName,nRuns,1,(Double_t)nRuns);
711     hCounterScalers->SetDirectory(0);
712     hCounterScalers->SetMinimum(1);
713     hCounterScalers->SetMaximum(maxScaler[0]);
714     hFromScalers.AddLast(hCounterScalers);
715     hOutput.AddLast(hCounterScalers);
716   }
717   //
718   //Fill histos for Scalers and QA
719   //
720   //loop on run list
721   for ( Int_t iRun = 0; iRun < runNumArray.GetEntries(); iRun++ ) {
722     
723     TString sRunNr = ((TObjString*)runNumArray.At(iRun))->GetString();
724     Int_t runNr = sRunNr.Atoi();
725     AliAnalysisTriggerScalers triggerScaler(runNr);
726     
727     //loop on trigger list
728     for ( Int_t iTrig = 0; iTrig < selectedTriggerArray->GetEntries(); iTrig++ ) {
729       
730       //loop on scaler list
731       for ( Int_t iScaler = 0; iScaler < nScaler; iScaler++ ) {
732         
733         //from Scalers
734         TGraph* graph = triggerScaler.PlotTrigger(selectedTriggerArray->At(iTrig)->GetName(),sScaler[iScaler].Data());
735         
736         sHistNameFull = Form("Scalers_%s_%s",selectedTriggerArray->At(iTrig)->GetName(),sScaler[iScaler].Data());
737         TH1* hist = (TH1*) hFromScalers.FindObject(sHistNameFull);
738         if ( !hist ) continue;
739         Double_t *tab = (Double_t*) graph->GetY();
740         if ( tab ) hist->SetBinContent(iRun+1,tab[0]);
741         else hist->SetBinContent(iRun+1,0);
742         hist->GetXaxis()->SetBinLabel(iRun+1,sRunNr.Data());
743         delete graph;
744
745         //from QA
746         if ( !(sScaler[iScaler].Contains("L2A")) ) continue;
747         TH1* histCounters = static_cast<TH1*>(ccol->Get("run",Form("run:%s/trigger:%s",sRunNr.Data(),selectedTriggerArray->At(iTrig)->GetName())));
748         
749         sHistNameFull = Form("QA_%s_%s",selectedTriggerArray->At(iTrig)->GetName(),sScaler[iScaler].Data());
750         hist = (TH1*) hFromQA.FindObject(sHistNameFull);
751         if (!hist) continue;
752         if ( histCounters ) hist->SetBinContent(iRun+1,histCounters->GetSumOfWeights());
753         else hist->SetBinContent(iRun+1,0);
754         hist->GetXaxis()->SetBinLabel(iRun+1,sRunNr.Data());
755         
756         delete histCounters;
757         
758       }//end loop on scaler list
759     }//end loop on trigger list
760
761     //loop on T0 trigger list
762     for ( Int_t iTrig = 0; iTrig < selectedT0TriggerArray->GetEntries(); iTrig++ ) {
763       
764       //from Scalers
765       TGraph* graph = triggerScaler.PlotTrigger(selectedT0TriggerArray->At(iTrig)->GetName(),sScaler[kL0B].Data());
766         
767       sHistNameFull = Form("Scalers_%s_%s",selectedT0TriggerArray->At(iTrig)->GetName(),sScaler[kL0B].Data());
768       TH1* hist = (TH1*) hFromScalers.FindObject(sHistNameFull);
769       if ( !hist ) continue;
770       Double_t *tab = (Double_t*) graph->GetY();
771       if ( tab ) hist->SetBinContent(iRun+1,tab[0]);
772       else hist->SetBinContent(iRun+1,0);
773       hist->GetXaxis()->SetBinLabel(iRun+1,sRunNr.Data());
774       delete graph;
775
776     }//end loop on T0 trigger list
777     
778   }//end loop on run list
779   
780   
781   //Set options for QA and Scalers histos
782    //loop on trigger list
783   for ( Int_t iTrig = 0; iTrig < selectedTriggerArray->GetEntries(); iTrig++ ) {
784     //loop on scaler list
785     for ( Int_t iScaler = 0; iScaler < nScaler; iScaler++ ) {
786       
787       sHistNameFull = Form("Scalers_%s_%s",((TObjString*) selectedTriggerArray->At(iTrig))->GetName(),sScaler[iScaler].Data());
788       TH1* histo = static_cast<TH1*> ( hFromScalers.FindObject(sHistNameFull) );
789       if (!histo) continue;
790       histo->LabelsOption("a");
791       histo->SetStats(kFALSE);
792       SetRunAxisRange(histo->GetXaxis());
793       
794       if ( !(sScaler[iScaler].Contains("L2A")) ) continue;
795       sHistNameFull = Form("QA_%s_%s",((TObjString*) selectedTriggerArray->At(iTrig))->GetName(),sScaler[iScaler].Data());
796       histo = static_cast<TH1*> ( hFromQA.FindObject(sHistNameFull) );
797       if (!histo) continue;
798       histo->LabelsOption("a");
799       histo->SetStats(kFALSE);
800       SetRunAxisRange(histo->GetXaxis());
801     }
802   }
803
804   
805   //Loop on histos from scalers and QA and create resulting histos from scalers
806   const Int_t nHisto = 3;
807   TString sHisto[nHisto] = {"L0BoverL0BC0TVX","L2AoverL0B","L2AQAoverSCALERS"};
808   TString sTitleHisto[nHisto] = {"L0B trigger / L0BC0TVX","L2A / L0B","L2A from QA / L2A from SCALERS"};
809   //  TString sHisto[nHisto] = {"L2AoverL0B","L2AQAoverSCALERS"};
810   
811   //loop on trigger list
812   for ( Int_t iTrig = 0; iTrig < selectedTriggerArray->GetEntries(); iTrig++ ) {
813
814     //C0TVX
815     TString sTrig = ( (TObjString*) selectedTriggerArray->At(iTrig) )->GetName();
816     TString sT0Trig = FindT0TriggerName(selectedT0TriggerArray,sTrig);
817
818     sHistNameFull = Form("Scalers_%s_L0B",sT0Trig.Data());
819
820     TH1* histo0 = static_cast<TH1*> ( hFromScalers.FindObject(sHistNameFull) );
821     if (!histo0) continue;
822
823     sHistNameFull = Form("Scalers_%s_L0B",((TObjString*) selectedTriggerArray->At(iTrig))->GetName());
824     TH1* histo1 = static_cast<TH1*> ( hFromScalers.FindObject(sHistNameFull) );
825     if (!histo1) continue;
826
827     sHistNameFull = Form("%s_%s",sHisto[0].Data(),((TObjString*) selectedTriggerArray->At(iTrig))->GetName());
828     TH1* histo10 = (TH1*) histo1->Clone(sHistNameFull);
829     histo10->SetTitle(sTitleHisto[0].Data());
830     histo10->Sumw2();
831     histo10->Divide(histo0);
832     histo10->SetMaximum(10);
833     histo10->SetMinimum(1e-5);
834     //outList.Add(histo10);
835     hOutput.AddLast(histo10);
836     //outList.Add(histo0);
837     //outList.Add(histo1);
838
839
840     //DEADTIME    
841     sHistNameFull = Form("Scalers_%s_L2A",((TObjString*) selectedTriggerArray->At(iTrig))->GetName());
842     TH1* histo2 = static_cast<TH1*> ( hFromScalers.FindObject(sHistNameFull) );
843     if (!histo2) continue;
844     
845     sHistNameFull = Form("%s_%s",sHisto[1].Data(),((TObjString*) selectedTriggerArray->At(iTrig))->GetName());
846     TH1* histo3 = (TH1*) histo2->Clone(sHistNameFull);
847     histo3->SetTitle(sTitleHisto[1]);
848     histo3->Sumw2();
849     histo3->Divide(histo1);
850     histo3->SetMaximum(1.2);
851     histo3->SetMinimum(1e-5);
852     //outList.Add(histo3);
853     hOutput.AddLast(histo3);
854     
855     //QA over Scalers
856     sHistNameFull = Form("QA_%s_L2A",((TObjString*) selectedTriggerArray->At(iTrig))->GetName());
857     TH1* histo4 = static_cast<TH1*> ( hFromQA.FindObject(sHistNameFull) );
858     if (!histo4) continue;
859     
860     sHistNameFull = Form("%s_%s",sHisto[2].Data(),((TObjString*) selectedTriggerArray->At(iTrig))->GetName());
861     TH1* histo5 = (TH1*) histo4->Clone(sHistNameFull);
862     histo5->SetTitle(sTitleHisto[2]);
863     histo5->Sumw2();
864     histo5->Divide(histo2);
865     histo5->SetMaximum(1.2);
866     histo5->SetMinimum(5e-1);
867     //outList.Add(histo5);
868     hOutput.AddLast(histo5);
869   }
870   
871   // Plot all on canvases (only canvases will be saved)
872   const Int_t nCanvases = nScaler + nHisto;
873   TString sCanvases[nCanvases];
874   for (Int_t iScaler = 0; iScaler < nScaler; iScaler++) sCanvases[iScaler] = sScaler[iScaler];
875   for (Int_t iHisto = 0; iHisto < nHisto; iHisto++) sCanvases[nScaler+iHisto] = sHisto[iHisto];
876   
877   //loop on canvases
878   for ( Int_t iCan = 0; iCan < nCanvases; iCan++) {
879     TCanvas* canvas = new TCanvas(sCanvases[iCan],sCanvases[iCan],200,10,600,600);
880     TLegend* leg  = new TLegend(0.72,0.7,0.9,0.85);
881     leg->SetBorderSize(1);
882     if ( iCan != 4 ) canvas->SetLogy();
883     TString optDraw = "e";
884     
885     //loop on trigger list
886     Int_t icolor = 1;
887     for ( Int_t iTrig = 0; iTrig < selectedTriggerArray->GetEntries(); iTrig++ ) {
888       
889       if ( iCan < nScaler ) sHistNameFull = Form("Scalers_%s_%s",selectedTriggerArray->At(iTrig)->GetName(),sCanvases[iCan].Data());
890       else sHistNameFull = Form("%s_%s",sCanvases[iCan].Data(),selectedTriggerArray->At(iTrig)->GetName());
891       TH1* histo1 = static_cast<TH1*> ( hOutput.FindObject(sHistNameFull) );
892       if (!histo1) continue;
893       
894       if ( icolor == 10 ) icolor++;
895       histo1->SetLineColor(icolor++);
896       histo1->Draw(optDraw);
897       optDraw = "esame";
898       
899       leg->AddEntry(histo1,selectedTriggerArray->At(iTrig)->GetName(),"l");
900     }
901     
902     leg->Draw();
903     outList.Add(canvas);
904     outCanList.Add(canvas);
905   }
906   
907   delete selectedTriggerArray;
908   delete selectedT0TriggerArray;
909   
910   file->Close();
911 }
912
913 //_____________________________________________________________________________
914 void trigEffQA(TString fileListName, TString outFilename = "", TString defaultStorage = "raw://", Bool_t doScalers = kFALSE)
915 {
916   ifstream inFile(fileListName.Data());
917   TObjArray fileNameArray, runNumArray;
918   fileNameArray.SetOwner();
919   runNumArray.SetOwner();
920   TString currString = "";
921   if (inFile.is_open()) {
922     while (! inFile.eof() ) {
923       currString.ReadLine(inFile,kTRUE); // Read line
924       if ( currString.IsNull() || ! currString.Contains(".root") ||
925           currString.BeginsWith("#") ) continue;
926       fileNameArray.AddLast(new TObjString(currString.Data()));
927       Int_t runNum = GetRunNumber(currString);
928       runNumArray.AddLast(new TObjString(Form("%i",runNum)));
929     }
930     inFile.close();
931   }
932   else {
933     printf("Fatal: cannot open input file %s\n",fileListName.Data());
934     return;
935   }
936   
937   runNumArray.Sort();
938   
939   TList outCanList, outList;
940   TrigEffTrending(runNumArray, fileNameArray, outCanList, outList);
941   if ( ! defaultStorage.IsNull() ) MaskTrending(runNumArray, defaultStorage, outCanList, outList);
942   if ( ! defaultStorage.IsNull() && doScalers ) ScalerTrending(runNumArray, "QAresults_Merged.root", defaultStorage, outCanList, outList);
943   
944   if ( outFilename.IsNull() ) return;
945   
946   TString outCanName = outFilename;
947   outCanName.ReplaceAll(".root",".pdf");
948   for ( Int_t ican=0; ican<outCanList.GetEntries(); ican++ ) {
949     TString canName = outCanName;
950     if ( ican == 0 ) canName.Append("("); // open pdf file
951     else if ( ican == outCanList.GetEntries()-1 ) canName.Append(")"); // close pdf file
952     static_cast<TCanvas*>(outCanList.At(ican))->Print(canName.Data());
953   }
954   // There is a bug when creating a pdf
955   // So create a ps and then convert via epstopdf
956   if ( outCanName.Contains(".ps") || outCanName.Contains(".eps") ) {
957     gSystem->Exec(Form("epstopdf %s", outCanName.Data()));
958     gSystem->Exec(Form("rm %s", outCanName.Data()));
959   }
960
961   TFile* outFile = new TFile(outFilename.Data(), "recreate");
962   outList.Write();
963   outFile->Close();
964 }