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