AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CreateEMCALRunQA.C
Go to the documentation of this file.
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <Riostream.h>
3 #include "AliEMCALGeometry.h"
4 #include <TColor.h>
5 #include <TStyle.h>
6 #include <TSystem.h>
7 #include <TDirectory.h>
8 #include <TF1.h>
9 #include <TFile.h>
10 #include <TH2F.h>
11 #include <TCanvas.h>
12 #include <TGraphErrors.h>
13 #include <TLegend.h>
14 #include <TTree.h>
15 #include <TPRegexp.h>
16 #include <TList.h>
17 #include <TObjString.h>
18 #include <TDatime.h>
19 #include <TError.h>
20 #include <AliLog.h>
21 #endif
22 
36 
37 Int_t DrawOccupancy(Long_t run, TString period, TString pass, TString fTrigger, TString system, TFile* f, TFile* fout, AliEMCALGeometry* geom, Int_t SavePlots);
38 Int_t DrawRun(Long_t run, TString period, TString pass, TString fTrigger, TFile *f, TFile* fout, Int_t SavePlots, Int_t nSM, Bool_t kFilter);
39 Int_t TrendingEMCALTree(Long_t RunId, TString fCalorimeter, TString system, TString period, TString pass, const int n, TList* TriggersList, TFile* f, TFile *fout, Int_t SavePlots);
40 
41 TH2F* FormatRunHisto(TH2F* aHisto, const char* title, const char* YTitle="");
42 TH2F* HistoPerMod(TH2F* name, const char* title);
43 TH2* AutoZoom(TH2* H, Option_t* aType="all", Int_t EntryMin=0);
44 int FindNumberOfSM(TFile* f, TString fTrigger, TString period);
45 
48 
49 //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
51 
52  const Int_t NRGBs = 5;
53  const Int_t NCont = 255;
54 
55  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
56  Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 };
57  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
58  Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 };
59  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
60  gStyle->SetNumberContours(NCont);
61 
62 }
63 
64 //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
66 
67  Double_t gaus;
68 
69  if(par[2] != 0.)
70  gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) / (2*par[2]*par[2]) );
71  else
72  gaus = 99999999.;
73 
74  Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
75 
76  return gaus+back;
77 
78 }
79 
80 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
82 
83  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) / (2*par[2]*par[2]) );
84  Double_t back = par[3] + par[4]*x[0];
85 
86  return gaus+back;
87 
88 }
89 
90 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
92 
93  Double_t levy;
94 
95  levy = par[0] * TMath::Exp( -par[1]/x[0]) * TMath::Power(x[0], -par[2]) ;
96 
97  return levy;
98 
99 }
100 
101 //--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
102 int CreateEMCALRunQA(const char* filename, TString RunId, TString period, TString pass, Int_t SavePlots = 0, Bool_t filter=0, TString fTrigger = "", TString system = "", TString fCalorimeter = "EMCAL"){
103 
104  QAPATH = TString(gSystem->Getenv("QAPATH"));
105  if(QAPATH.IsNull())
106  QAPATH = QAPATHF;
107 
108  if(! QAPATH.BeginsWith("./"))
109  QAPATH = QAPATH + RunId + "/";
110 
111  AliLog::SetGlobalLogLevel(AliLog::kError);
112  TFile *f = new TFile(filename);
113  AliLog::SetGlobalLogLevel(AliLog::kInfo);
114 
115  if(f->IsZombie()){
116  Error(__FUNCTION__, Form("Error openning the input file %s", filename));
117  return -1;
118  }
119 
120  TList* TriggersList = new TList();
121 
122  if(fTrigger==""){
123  TPMERegexp name_re("CaloQA_\\w+");
124  TObjLink* link = f->GetListOfKeys()->FirstLink();
125 
126  while(link){
127  TString name = link->GetObject()->GetName();
128  if(name_re.Match(name)){
129  TriggersList->Add(link->GetObject());
130  if(TString(filename).Contains("barrel") && ! name.Contains("default")) TriggersList->Remove(link->GetObject());
131  if(TString(filename).Contains("outer") && ! name.Contains("EMC")) TriggersList->Remove(link->GetObject());
132  }
133  link = link->Next();
134  }
135  }
136  else{
137  TriggersList->Add(new TObjString(fTrigger.Data()));
138  }
139 
140  if(!TriggersList->GetEntries()){
141  Error(__FUNCTION__, "No trigger found!");
142  return -2;
143  }
144 
145  int nSM = FindNumberOfSM(f, ((TObjString*)TriggersList->Last())->GetString(), period);
146 
147  if(nSM<0){
148  Error(__FUNCTION__, "Could not find the number of super modules!");
149  return -3;
150  }
151  Info(__FUNCTION__, Form("%i super modules were discovered", nSM));
152 
153  TString GeomName;
154  if(nSM <= 4) { nSM = 4; GeomName = "EMCAL_FIRSTYEARv1"; }
155  else if(nSM <= 10){ nSM = 10; GeomName = "EMCAL_COMPLETEv1"; }
156  else if(nSM <= 12){ nSM = 12; GeomName = "EMCAL_COMPLETE12SMv1"; }
157  else { nSM = 20; GeomName = "EMCAL_COMPLETE12SMv1_DCAL_8SM";}
158 
159  AliEMCALGeometry *geom = new AliEMCALGeometry(GeomName.Data(), "EMCAL");
160  Info(__FUNCTION__, Form("Using %i super modules and the Geometry %s", nSM, GeomName.Data()));
161 
162  TFile *fout = new TFile(TString( QAPATH + period+"_"+pass + fTrigger+"_"+ (Long_t)RunId.Atoi() +"_QAplots.root").Data(), "RECREATE");
163 
164  if((system.IsNull()) && (period.EndsWith("h")))
165  system = "PbPb";
166 
167  Int_t ret=0;
168  TIter next(TriggersList);
169  while(TObject *obj = next()){
170  fTrigger= TString(obj->GetName());
171  ret -= DrawOccupancy(RunId.Atoi(), period, pass, fTrigger, system, f, fout, geom, SavePlots);
172  ret -= DrawRun(RunId.Atoi(), period, pass, fTrigger, f, fout, SavePlots, nSM, filter);
173  }
174  ret-= TrendingEMCALTree(RunId.Atoi(), fCalorimeter, system, period, pass, nSM, TriggersList, f, fout, SavePlots);
175 
176  f->Close();
177  fout->Close();
178  delete f;
179  delete geom;
180 
181  return ret;
182 }
183 
184 //----------------------------------------------------------------------------------------------------------------------------------------------------------
185 Int_t DrawOccupancy(Long_t run, TString period, TString pass, TString fTrigger, TString system, TFile* f, TFile* fout, AliEMCALGeometry* geom, Int_t SavePlots){
186 
187  set_plot_style();
188  gStyle->SetOptStat(0);
189 
190  TH1::AddDirectory(kFALSE);
191  TH2D *hEnergyMapReal = new TH2D("hEnergyMapReal", "", 96, -48, 48, 120, -0, 120);
192  TH2D *hOccupancyMapReal = new TH2D("hOccupancyMapReal", "", 96, -48, 48, 120, -0, 120);
193 
194  hEnergyMapReal->SetXTitle("#eta (bin)");
195  hEnergyMapReal->SetYTitle("#varphi (bin)");
196  hEnergyMapReal->SetZTitle("E (GeV)/event");
197  hEnergyMapReal->GetYaxis()->SetTitleOffset(1.2);
198  hEnergyMapReal->GetZaxis()->SetLabelSize(0.02);
199  hEnergyMapReal->GetZaxis()->SetTitleOffset(1.36);
200 
201  hOccupancyMapReal->SetXTitle("#eta (bin)");
202  hOccupancyMapReal->SetYTitle("#varphi (bin)");
203  hOccupancyMapReal->GetYaxis()->SetTitleOffset(1.2);
204  hOccupancyMapReal->GetZaxis()->SetLabelSize(0.02);
205 
206  Int_t nSupMod, nModule, nIphi, nIeta;
207  Int_t iphi, ieta;
208  Int_t realbineta=0;
209  Int_t realbinphi=0;
210 
211  // NO MASK
212  Int_t mask[1] = {2222222};
213 
214  TH2F *hCellAmplitude;
215  TH1F *hNEvents;
216  Int_t Events;
217  Int_t n=0;
218 
219  TString direct;
220  if(!fTrigger.Contains("QA")){
221  direct = "CaloQA_";
222  }
223  direct += fTrigger;
224 
225  Bool_t dirok = f->cd(direct);
226  if(!dirok){
227  Error(__FUNCTION__, Form("No input drectory %s", direct.Data()));
228  return -1;
229  }
230 
231  TList *outputList = (TList*)gDirectory->Get(direct);
232  if(!outputList){
233  Error(__FUNCTION__, "No input list! ");
234  return -1;
235  }
236  outputList->SetOwner();
237 
238  fout->mkdir(Form("%s/%s/%ld/%s/%s", period.Data(), pass.Data(), run, "RunLevelQA", fTrigger.Data()));
239  fout->cd();
240  fout->Cd(Form("%s/%s/%ld/%s/%s", period.Data(), pass.Data(), run, "RunLevelQA", fTrigger.Data()));
241 
242  hNEvents = (TH1F *)outputList->FindObject("hNEvents");
243  if(!hNEvents){
244  Error(__FUNCTION__, Form("hNEvent histogram not found for trigger %s ", fTrigger.Data()));
245  return -2;
246  }
247 
248  Events = (Int_t)hNEvents->GetEntries();
249  if(Events==0){
250  Error(__FUNCTION__, Form("No event in trigger %s", fTrigger.Data()));
251  return -3;
252  }
253 
254  TH1F* hE = (TH1F *)outputList->FindObject("EMCAL_hE");
255  if(!hE || (!hE->GetEntries())){
256  Error(__FUNCTION__, Form("hE Histogram not found or it is empty for trigger %s, EMCal was not in this run?", fTrigger.Data()));
257  return -4;
258  }
259 
260  Double_t Eth=1;
261  if(system=="PbPb"){
262  Eth = 5.;
263  if(fTrigger.Contains("EMC"))
264  Eth=20.;
265  }
266  if(system=="pp"){
267  Eth = 1.;
268  if(fTrigger.Contains("EMC"))
269  Eth=5.;
270  }
271 
272  hCellAmplitude = (TH2F *)outputList->FindObject("EMCAL_hAmpId");
273 
274  for(Int_t i = 0; i < geom->GetNCells() ; i++){
275  Double_t Esum = 0;
276  Double_t Nsum = 0;
277 
278  for(Int_t j = 1; j <= hCellAmplitude->GetNbinsX(); j++){
279  Double_t E = hCellAmplitude->GetXaxis()->GetBinCenter(j);
280  Double_t N = hCellAmplitude->GetBinContent(j, i+1);
281 
282  if(E < 0.3)
283  continue;
284 
285  if(E <= Eth){
286  Esum += E*N;
287  Nsum += N;
288  }
289  }
290 
291  Int_t absId = i;
292  if(n!=0){
293  if(mask[n]<=mask[n-1])
294  Warning(__FUNCTION__, "The list of bad cells is not sorted !!");
295  }
296 
297  if(i==mask[n]){
298  n++ ;
299  continue;
300  } // skip bad cells
301 
302  geom->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
303  geom->GetCellPhiEtaIndexInSModule(nSupMod, nModule, nIphi, nIeta, iphi, ieta);
304 
305  realbinphi = 120-(nSupMod/2)*24 -iphi -1;
306  if(nSupMod%2==0)
307  realbineta= 48-ieta -1;
308  if(nSupMod%2==1)
309  realbineta= -ieta -1;
310 
311  hEnergyMapReal->Fill(realbineta, realbinphi, Esum/(Double_t)Events);
312  // hOccupancyMapReal->Fill(realbineta, realbinphi, Nsum/(Double_t)Events);
313  }
314 
315  cout <<" Run: " << run << " trigger: " << fTrigger << " N_events: "<<Events<<endl;
316 
317  TPMERegexp r("_\\w+");
318  TString Energy; Energy = QAPATH + "MapEnergy" + fTrigger(r) + ".pdf";
319  TString Energy2; Energy2 = QAPATH + "MapEnergy" + fTrigger(r) + ".png";
320  TString Entries; Entries = QAPATH + "MapEntries" + fTrigger(r) + ".pdf";
321  TString Entries2; Entries2 = QAPATH + "MapEntries" + fTrigger(r) + ".png";
322 
323  TCanvas *c1 = new TCanvas("Energymap", "Energy Map", 600, 600);
324  c1->SetFillColor(0);
325  c1->SetGrid();
326  c1->SetRightMargin(0.14);
327 
328  TString title = "run ";
329  title += run ;
330  if(fTrigger.Contains("EMC"))
331  title += " EMC ";
332  else
333  title += " MB ";
334  title += " Summed energy map";
335 
336  hEnergyMapReal->SetTitle(title);
337  AutoZoom(hEnergyMapReal, "miny")->DrawCopy("colz");
338 
339  if(nSupMod<=12){
340  if(SavePlots==2)
341  c1->SaveAs(Energy);
342  if(SavePlots)
343  c1->SaveAs(Energy2);
344  c1->Write();
345  }
346  delete c1;
347 
348  // TCanvas *c2 = new TCanvas("Occupancy", "Occupancy Map", 600, 600);
349  // c2->SetFillColor(0);
350  // c2->SetGrid();
351  // c2->SetRightMargin(0.14);
352  // TString title2 = "run ";
353  // title2 += run ;
354  // if(fTrigger.Contains("EMC")){ title2 += " EMC ";} else{ title2 += " MB ";}
355  // title2 += " Occupancy map";
356  // hOccupancyMapReal->SetTitle(title2);
357  // AutoZoom(hOccupancyMapReal, "miny")->DrawCopy("colz");
358  // if(SavePlots==2) c2->SaveAs(Entries);
359  // if(SavePlots) c2->SaveAs(Entries2);
360  // c2->Write();
361  // delete c2;
362 
363  if(outputList)
364  outputList->Delete();
365  delete hEnergyMapReal;
366  delete hOccupancyMapReal;
367 
368  return 0;
369 }
370 
371 //-----------------------------------------------------------------------------------------------------------------------
372 Int_t DrawRun(const Long_t run, TString period, TString pass, TString fTrigger, TFile *f, TFile *fout, Int_t SavePlots, Int_t nSM, Bool_t kFilter){
373 
374  TString direct;
375  if(!fTrigger.Contains("QA")){
376  direct = "CaloQA_";
377  }
378  direct += fTrigger;
379 
380  f->cd(direct);
381  if(!direct){
382  Error(__FUNCTION__, Form("No input directory %s", direct.Data()));
383  return -1;
384  }
385 
386  TList *outputList = (TList*)gDirectory->Get(direct);
387  if(!outputList){
388  Error(__FUNCTION__, Form("No input list! %s", direct.Data()));
389  return -2;
390  }
391  outputList->SetOwner();
392 
393  if(kFilter){
394  fout->mkdir(Form("%s/%s/%ld/%s", period.Data(), pass.Data(), run, fTrigger.Data()));
395  fout->cd();
396  fout->Cd(Form("%s/%s/%ld/%s", period.Data(), pass.Data(), run, fTrigger.Data()));
397  outputList->Write();
398  }
399  fout->cd();
400  fout->Cd(Form("%s/%s/%ld/%s/%s", period.Data(), pass.Data(), run, "RunLevelQA", fTrigger.Data()));
401 
402  set_plot_style();
403  gStyle->SetPalette(1);
404  set_plot_style();
405  // gStyle->SetOptStat(1);
406 
407  TH1::AddDirectory(kFALSE);
408  TString outfilename;
409  TString outfilename2;
410  const char* legend="";
411  TPMERegexp r("_\\w+");
412 
413  if(fTrigger.Contains("EMC"))
414  legend = Form(" Run %d EMC ", (int)run);
415  else
416  legend = Form(" Run %d MB ", (int)run);
417 
418  TH1F* hNEvents =(TH1F *)outputList->FindObject("hNEvents");
419  if(!hNEvents){
420  Error(__FUNCTION__, Form("hNEvent histogram not found for trigger %s ", fTrigger.Data()));
421  return -3;
422  }
423 
424  Int_t Events = (Int_t)hNEvents->GetEntries();
425  if(Events==0){
426  Error(__FUNCTION__, Form("No event in trigger %s", fTrigger.Data()));
427  return -4 ;
428  }
429 
430  TH1F* hE =(TH1F *)outputList->FindObject("EMCAL_hE");
431  if(!hE || (!hE->GetEntries())){
432  Error(__FUNCTION__, Form("hE Histogram not found or it is empty for trigger %s, EMCal was not in this run?", fTrigger.Data()));
433  return -4;
434  }
435 
436  TCanvas* c00 = new TCanvas("Occupancy map", "Occupancy", 600, 600);
437  c00->SetLogz();
438  c00->SetFillColor(0);
439  c00->SetBorderSize(0);
440  c00->SetFrameBorderMode(0);
441 
442  TH2F* hOccupancyMapReal =(TH2F *)outputList->FindObject("EMCAL_hEtaPhi");
443  if(!hOccupancyMapReal){
444  Error(__FUNCTION__, Form("EMCAL_hEtaPhi: Histogram for trigger %s not found!", fTrigger.Data()));
445  return -5;
446  }
447  FormatRunHisto(hOccupancyMapReal, Form("Occupancy%s", legend), "#varphi (bin)");
448  hOccupancyMapReal->SetXTitle("#eta (bin)");
449 
450  AutoZoom(hOccupancyMapReal, "all")->DrawCopy("colz");
451  outfilename = QAPATH + "OccupancyMap" + fTrigger(r) + ".pdf" ;
452  outfilename2 = QAPATH + "OccupancyMap" + fTrigger(r) + ".png" ;
453 
454  if(SavePlots==2)
455  c00->SaveAs(outfilename);
456  if(SavePlots)
457  c00->SaveAs(outfilename2);
458  c00->Write();
459  delete c00;
460 
461  TCanvas* c0 = new TCanvas("Grid Cell", "Occupancy from Grid Cells", 600, 600);
462  c0->SetLogz();
463  c0->SetFillColor(0);
464  c0->SetBorderSize(0);
465  c0->SetFrameBorderMode(0);
466 
467  TH2F* hGridCells =(TH2F *)outputList->FindObject("EMCAL_hGridCells");
468  if(!hGridCells){
469  Error(__FUNCTION__, Form("EMCAL_hGridCells: Histogram for trigger %s not found!", fTrigger.Data()));
470  return -5;
471  }
472  FormatRunHisto(hGridCells, Form("Occupancy from Grid Cells%s", legend), "#varphi row columns");
473  hGridCells->SetXTitle("#eta row columns");
474  AutoZoom(hGridCells, "all")->DrawCopy("colz");
475 
476  outfilename = QAPATH + "GridCells" + fTrigger(r) + ".pdf" ;
477  outfilename2 = QAPATH + "GridCells" + fTrigger(r) + ".png" ;
478 
479  if(SavePlots==2)
480  c0->SaveAs(outfilename);
481  if(SavePlots)
482  c0->SaveAs(outfilename2);
483  c0->Write();
484  delete c0;
485 
486  TCanvas* c1 = new TCanvas("TimeVsE", "Cluster Time vs. Energy", 600, 600);
487  c1->SetLogz();
488  c1->SetFillColor(0);
489  c1->SetBorderSize(0);
490  c1->SetFrameBorderMode(0);
491 
492  TH2F* hClusterTimeEnergy =(TH2F *)outputList->FindObject("EMCAL_hClusterTimeEnergy");
493  if(!hClusterTimeEnergy){
494  Error(__FUNCTION__, Form("EMCAL_hClusterTimeEnergy: Histogram for trigger %s not found!", fTrigger.Data()));
495  return -5;
496  }
497  FormatRunHisto(hClusterTimeEnergy, Form("Time vs. Energy%s", legend), "EMCal ToF (ns)");
498  AutoZoom(hClusterTimeEnergy, "maxx")->DrawCopy("colz");
499 
500  outfilename = QAPATH + "TimeRun" + fTrigger(r) + ".pdf" ;
501  outfilename2 = QAPATH + "TimeRun" + fTrigger(r) + ".png" ;
502 
503  if(SavePlots==2)
504  c1->SaveAs(outfilename);
505  if(SavePlots)
506  c1->SaveAs(outfilename2);
507  c1->Write();
508  delete c1;
509 
510  if(pass!="simu"){
511  TCanvas* c1a = new TCanvas("TimeVsAbsId", "Cell Id vs. time ", 600, 600);
512  c1a->SetLogz();
513  c1a->SetFillColor(0);
514  c1a->SetBorderSize(0);
515  c1a->SetFrameBorderMode(0);
516 
517  TH2F* hTimeAbsId =(TH2F *)outputList->FindObject("EMCAL_hTimeId");
518  if(!hTimeAbsId){
519  Error(__FUNCTION__, Form("EMCAL_hTimeId: Histogram for trigger %s not found!", fTrigger.Data()));
520  return -5;
521  }
522  FormatRunHisto(hTimeAbsId, Form("Cell ID vs. Time%s", legend), "Cell AbsID");
523  hTimeAbsId->RebinY(16);
524  hTimeAbsId->SetMinimum(32);
525  AutoZoom(hTimeAbsId, "maxx")->DrawCopy("colz");
526 
527  outfilename = QAPATH + "TimeAbsIdRun" + fTrigger(r) + ".pdf" ;
528  outfilename2 = QAPATH + "TimeAbsRun" + fTrigger(r) + ".png" ;
529 
530  if(SavePlots==2)
531  c1a->SaveAs(outfilename);
532  if(SavePlots)
533  c1a->SaveAs(outfilename2);
534  c1a->Write();
535  delete c1a;
536  }
537 
538  // TCanvas* c1aLG = new TCanvas("TimeVsAbsIdLG", "Cell Id vs. time LG ", 600, 600);
539  // c1aLG->SetLogz();
540  // c1aLG->SetFillColor(0);
541  // c1aLG->SetBorderSize(0);
542  // c1aLG->SetFrameBorderMode(0);
543 
544  // TH2F* hTimeAbsIdLG =(TH2F *)outputList->FindObject("EMCAL_hTimeIdLG");
545  // if(!hTimeAbsIdLG){
546  // Error(__FUNCTION__, Form("EMCAL_hTimeIdLG: Histogram for trigger %s not found!", fTrigger.Data()));
547  // return -5;
548  // }
549  // FormatRunHisto(hTimeAbsIdLG, Form("Cell ID vs. Time LG%s", legend), "cell AbsID");
550 
551  // AutoZoom(hTimeAbsIdLG, "maxx")->DrawCopy("colz");
552  // outfilename = QAPATH + "TimeAbsIdRunLG" + fTrigger(r) + ".pdf" ;
553  // outfilename2 = QAPATH + "TimeAbsRunLG" + fTrigger(r) + ".png" ;
554 
555  // if(SavePlots==2) c1a->SaveAs(outfilename);
556  // if(SavePlots) c1a->SaveAs(outfilename2);
557  // c1a->Write();
558  // delete c1a;
559 
560  TCanvas* c1b = new TCanvas("EVsAbsId", "Cell Id vs. Energy", 600, 600);
561  c1b->SetLogz();
562  c1b->SetFillColor(0);
563  c1b->SetBorderSize(0);
564  c1b->SetFrameBorderMode(0);
565 
566  TH2F* hEAbsId =(TH2F *)outputList->FindObject("EMCAL_hAmpId");
567  if(!hEAbsId){
568  Error(__FUNCTION__, Form("EMCAL_hAmpId: Histogram for trigger %s not found!", fTrigger.Data()));
569  return -5;
570  }
571  FormatRunHisto(hEAbsId, Form("Cell ID vs. Energy%s", legend), "Cell AbsID");
572  hEAbsId->RebinY(16);
573  hEAbsId->SetMinimum(32);
574  AutoZoom(hEAbsId, "maxx")->DrawCopy("colz");
575 
576  outfilename = QAPATH + "EAbsIdRun" + fTrigger(r) + ".pdf" ;
577  outfilename2 = QAPATH + "EAbsIdRun" + fTrigger(r) + ".png" ;
578 
579  if(SavePlots==2)
580  c1b->SaveAs(outfilename);
581  if(SavePlots)
582  c1b->SaveAs(outfilename2);
583  c1b->Write();
584  delete c1b;
585 
586  // TCanvas* c1bLG = new TCanvas("EVsAbsIdLG", "Cell Id vs. E LG", 600, 600);
587  // c1bLG->SetLogz();
588  // c1bLG->SetFillColor(0);
589  // c1bLG->SetBorderSize(0);
590  // c1bLG->SetFrameBorderMode(0);
591 
592  // TH2F* hEAbsIdLG =(TH2F *)outputList->FindObject("EMCAL_hAmpIdLG");
593  // if(!hEAbsIdLG){
594  // Error(__FUNCTION__, Form("EMCAL_hAmpIdLG: Histogram for trigger %s not found!", fTrigger.Data()));
595  // return -5;
596  // }
597  // FormatRunHisto(hEAbsIdLG, Form("Cell ID vs. Energy%s", legend), "Cell AbsID ");
598 
599  // AutoZoom(hEAbsIdLG, "maxx")->DrawCopy("colz");
600  // outfilename = QAPATH + "EAbsIdRunLG" + fTrigger(r) + ".pdf" ;
601  // outfilename2 = QAPATH + "EAbsIdRunLG" + fTrigger(r) + ".png" ;
602 
603  // if(SavePlots==2) c1bLG->SaveAs(outfilename);
604  // if(SavePlots) c1bLG->SaveAs(outfilename2);
605  // c1bLG->Write();
606  // delete c1bLG;
607 
608  if(pass!="muon_calo_pass1"){
609  TCanvas * c2 = new TCanvas("ClusterVsTrack ", "Correlation calo Mult vs. Track Multiplicity", 600, 600);
610  c2->SetLogz();
611  c2->SetFillColor(0);
612  c2->SetBorderSize(0);
613  c2->SetFrameBorderMode(0);
614 
615  TH2F* hClusterVsTrack =(TH2F *)outputList->FindObject("EMCAL_hCaloTrackMNClusters");
616  FormatRunHisto(hClusterVsTrack, Form("#it{N} clusters vs. #it{N} tracks%s", legend), "#it{N} clusters");
617  AutoZoom(hClusterVsTrack, "maxx, maxy", 1)->DrawCopy("colz");
618 
619  outfilename = QAPATH + "CaloTrackMult" + fTrigger(r) + ".pdf";
620  outfilename2 = QAPATH + "CaloTrackMult" + fTrigger(r) + ".png";
621 
622  // if(pass!="muon_calo_pass1"){
623  if(SavePlots==2)
624  c2->SaveAs(outfilename);
625  if(SavePlots)
626  c2->SaveAs(outfilename2);
627  c2->Write();
628  delete c2;
629 
630  TCanvas* c3 = new TCanvas("ClusterEVsTrack ", "Correlation #it{E} calo vs. Track Multiplicity", 600, 600);
631  c3->SetLogz();
632  c3->SetFillColor(0);
633  c3->SetBorderSize(0);
634  c3->SetFrameBorderMode(0);
635 
636  TH2F* hClusterEVsTrack =(TH2F*)outputList->FindObject("EMCAL_hCaloTrackMEClusters");
637  FormatRunHisto(hClusterEVsTrack, Form("Sum #it{E} clusters vs. #it{N} tracks%s", legend), "#Sigma#it{E} clusters");
638  AutoZoom(hClusterEVsTrack, "maxx, maxy", 1)->DrawCopy("colz");
639 
640  outfilename = QAPATH + "ETrackMult" + fTrigger(r) + ".pdf";
641  outfilename2 = QAPATH + "ETrackMult" + fTrigger(r) + ".png";
642 
643  // if(pass!="muon_calo_pass1"){
644  if(SavePlots==2)
645  c3->SaveAs(outfilename);
646  if(SavePlots)
647  c3->SaveAs(outfilename2);
648  c3->Write();
649  // }
650  delete c3;
651  }
652 
653  TCanvas* c4 = new TCanvas("ClusterEVsV0 ", "Correlation #it{E} calo vs. V0 signal", 600, 600);
654  c4->SetLogz();
655  c4->SetFillColor(0);
656  c4->SetBorderSize(0);
657  c4->SetFrameBorderMode(0);
658 
659  TH2F* hClusterEVsV0S =(TH2F*)outputList->FindObject("EMCAL_hCaloV0SEClusters");
660  FormatRunHisto(hClusterEVsV0S, Form("Sum #it{E} clusters vs. V0 signal%s", legend), "#Sigma#it{E} clusters");
661  AutoZoom(hClusterEVsV0S, "maxx, maxy", 1)->DrawCopy("colz");
662 
663  outfilename = QAPATH +"EVsV0s" + fTrigger(r) + ".pdf";
664  outfilename2 = QAPATH +"EVsV0s" + fTrigger(r) + ".png";
665 
666  if(SavePlots==2)
667  c4->SaveAs(outfilename);
668  if(SavePlots)
669  c4->SaveAs(outfilename2);
670  c4->Write();
671  delete c4;
672 
673  TCanvas* c6 = new TCanvas("SumCellEvsSM", "sum cell energy vs. SM", 600, 600);
674  c6->SetLogz();
675  c6->SetFillColor(0);
676  c6->SetBorderSize(0);
677  c6->SetFrameBorderMode(0);
678 
679  TH2F* hSumCellEVsSM =(TH2F*)outputList->FindObject("EMCAL_hSumCellsAmp_Mod");
680  FormatRunHisto(hSumCellEVsSM, Form("Sum #it{E} cells vs. SM%s", legend), "SM number");
681  AutoZoom(hSumCellEVsSM, "maxx, maxy", 1)->DrawCopy("colz");
682 
683  outfilename = QAPATH +"SumCellEVsSM" + fTrigger(r) + ".pdf";
684  outfilename2 = QAPATH +"SumCellEVsSM" + fTrigger(r) + ".png";
685 
686  if(SavePlots==2)
687  c6->SaveAs(outfilename);
688  if(SavePlots)
689  c6->SaveAs(outfilename2);
690  c6->Write();
691  delete c6;
692 
693  TCanvas* c5 = new TCanvas("CellsperCluster", "Nb of cells per cluster for each SM", 600, 600);
694  c5->SetLogz();
695  c5->SetFillColor(0);
696  c5->SetBorderSize(0);
697  c5->SetFrameBorderMode(0);
698 
699  Bool_t mod3=0;
700  if(nSM%3)
701  mod3=1;
702  c5->Divide(3, (nSM/3)+mod3);
703 
704  for(int ism = 0; ism < nSM; ism++){
705  c5->cd(ism+1);
706  gPad->SetLogz();
707 
708  if(TString(Form("Nb cells per clus%s Mod %d", legend, ism)).Length() > 60){
709  Error(__FUNCTION__, "Title too long!");
710  return -6;
711  }
712 
713  AutoZoom(HistoPerMod((TH2F*)outputList->FindObject(Form("EMCAL_hNCellsPerCluster_Mod%i", ism)), Form("Nb of cells per cluster%s Mod %d", legend, ism)), "all", 1)->DrawCopy("colz");
714  }
715 
716  outfilename = QAPATH + "CellsperClusterSM" + fTrigger(r) + ".pdf";
717  outfilename2 = QAPATH + "CellsperClusterSM" + fTrigger(r) + ".png";
718 
719  if(SavePlots==2)
720  c5->SaveAs(outfilename);
721  if(SavePlots)
722  c5->SaveAs(outfilename2);
723  c5->Write();
724  delete c5;
725  if(outputList)
726  outputList->Delete();
727 
728  return 0;
729 }
730 
731 //----------------------------------------------------------------------------------------------------------------------------------
732 Int_t TrendingEMCALTree(Long_t RunId, TString fCalorimeter, TString system, TString period, TString pass, int n, TList* TriggersList, TFile* f, TFile *fout, Int_t SavePlots){
733 
734  TString fTrigger="";
735  TString aCalorimeter;
736  if(n<=12)
737  aCalorimeter = fCalorimeter;
738  else
739  aCalorimeter = TString("EMCAL_and_DCAL");
740 
741  TDatime now;
742 
743  Double_t Nevent=0 ;
744  Double_t xe=0.5;
745  Double_t NTClusters=0; // total number of clusters
746  Double_t NTClustersRMS=0;
747  Double_t NCClusters=0; // number of charged clusters
748  Double_t NCClustersRMS=0;
749  Double_t NMatchClustersP=0;
750  Double_t NMatchClustersPRMS=0;
751 
752  Double_t CellMean=0;
753  Double_t CellRMS=0;
754  Double_t ClusterMean=0;
755  Double_t ClusterRMS=0;
756  Double_t EtotalMean=0;
757  Double_t EtotalRMS=0;
758 
759  Double_t CellPerClusterMean=0;
760  Double_t CellPerClusterRMS=0;
761 
762  Double_t mPDG = 134.9766;
763  Double_t Npi0EMCAL=0;
764  Double_t Npi0EMCALErr=0;
765  Double_t MeanPosEMCAL=0;
766  Double_t MeanPosEMCALErr=0;
767  Double_t WidthEMCAL=0;
768  Double_t WidthEMCALErr=0;
769  Double_t Chi2NdfPi0EMCAL=0;
770  Double_t NggEMCAL=0;
771  Double_t NggEMCALErr=0;
772  Double_t SignifEMCAL=0; // !S/(S+B)
773  Double_t SignifEMCALErr=0; // !S/(S+B)
774 
775  Double_t Npi0DCAL=0;
776  Double_t Npi0DCALErr=0;
777  Double_t MeanPosDCAL=0;
778  Double_t MeanPosDCALErr=0;
779  Double_t WidthDCAL=0;
780  Double_t WidthDCALErr=0;
781  Double_t Chi2NdfPi0DCAL=0;
782  Double_t NggDCAL=0;
783  Double_t NggDCALErr=0;
784  Double_t SignifDCAL=0; // !S/(S+B)
785  Double_t SignifDCALErr=0; // !S/(S+B)
786 
787  TFile* ftree = new TFile(Form("%s/trending.root", QAPATH.Data()), "RECREATE");
788 
789  TTree *tree = new TTree("trending", "Trending QA Tree");
790  tree->Branch("fDate", &now);
791  tree->Branch("fCalorimeter", &aCalorimeter);
792  tree->Branch("system", &system);
793  tree->Branch("period", &period);
794  tree->Branch("pass", &pass);
795  tree->Branch("fTrigger", &fTrigger);
796  tree->Branch("run", &RunId, "run/I");
797  tree->Branch("xe", &xe, "xe/D");
798 
799  tree->Branch("Nevent", &Nevent, "Nevent/D");
800  tree->Branch("NMatchClustersP", &NMatchClustersP, "NMatchClustersP/D");
801  tree->Branch("NMatchClustersPRMS", &NMatchClustersPRMS, "NMatchClustersPRMS/D");
802  tree->Branch("CellMean", &CellMean, "CellMean/D");
803  tree->Branch("CellRMS", &CellRMS, "CellRMS/D");
804  tree->Branch("ClusterMean", &ClusterMean, "ClusterMean/D");
805  tree->Branch("ClusterRMS", &ClusterRMS, "ClusterRMS/D");
806  tree->Branch("EtotalMean", &EtotalMean, "EtotalMean/D");
807  tree->Branch("EtotalRMS", &EtotalRMS, "EtotalRMS/D");
808 
809  tree->Branch("CellPerClusterMean", &CellPerClusterMean, "CellPerClusterMean/D");
810  tree->Branch("CellPerClusterRMS", &CellPerClusterRMS, "CellPerClusterRMS/D");
811 
812  tree->Branch("Npi0EMCAL", &Npi0EMCAL, "Npi0EMCAL/D");
813  tree->Branch("Npi0EMCALErr", &Npi0EMCALErr, "Npi0EMCALErr/D");
814  tree->Branch("Npi0DCAL", &Npi0DCAL, "Npi0DCAL/D");
815  tree->Branch("Npi0DCALErr", &Npi0DCALErr, "Npi0DCALErr/D");
816  tree->Branch("MeanPosEMCAL", &MeanPosEMCAL, "MeanPosEMCAL/D");
817  tree->Branch("MeanPosEMCALErr", &MeanPosEMCALErr, "MeanPosEMCALErr/D");
818  tree->Branch("MeanPosDCAL", &MeanPosDCAL, "MeanPosDCAL/D");
819  tree->Branch("MeanPosDCALErr", &MeanPosDCALErr, "MeanPosDCALErr/D");
820  tree->Branch("WidthEMCAL", &WidthEMCAL, "WidthEMCAL/D");
821  tree->Branch("WidthEMCALErr", &WidthEMCALErr, "WidthEMCALErr/D");
822  tree->Branch("WidthDCAL", &WidthDCAL, "WidthDCAL/D");
823  tree->Branch("WidthDCALErr", &WidthDCALErr, "WidthDCALErr/D");
824  tree->Branch("Chi2NdfPi0EMCAL", &Chi2NdfPi0EMCAL, "Chi2NdfPi0EMCAL/D");
825  tree->Branch("Chi2NdfPi0DCAL", &Chi2NdfPi0DCAL, "Chi2NdfPi0DCAL/D");
826  tree->Branch("NggEMCAL", &NggEMCAL, "NggEMCAL/D");
827  tree->Branch("NggEMCALErr", &NggEMCALErr, "NggEMCALErr/D");
828  tree->Branch("NggDCAL", &NggDCAL, "NggDCAL/D");
829  tree->Branch("NggDCALErr", &NggDCALErr, "NggDCALErr/D");
830  tree->Branch("SignifEMCAL", &SignifEMCAL, "SignifEMCAL/D");
831  tree->Branch("SignifDCALErr", &SignifDCALErr, "SignifDCALErr/D");
832  tree->Branch("SignifEMCAL", &SignifEMCAL, "SignifEMCAL/D");
833  tree->Branch("SignifDCALErr", &SignifDCALErr, "SignifDCALErr/D");
834 
835  tree->Branch("nSM", &n, "nSM/I");
836 
837  int nMax = 22;
838  Double_t CellMeanSM[nMax];
839  Double_t CellRMSSM[nMax];
840  Double_t ClusterMeanSM[nMax];
841  Double_t ClusterTotSM[nMax];
842  Double_t ClusterRMSSM[nMax];
843  Double_t EtotalMeanSM[nMax]; // mean total energy deposited per event
844  Double_t EtotalRMSSM[nMax];
845  Double_t CellPerClusterMeanSM[nMax];
846  Double_t CellPerClusterRMSSM[nMax];
847  Double_t ECell1MeanSM[nMax]; // total energy deposited per event without 1-cell clusters
848  Double_t ECell1RMSSM[nMax];
849 
850  Double_t MeanPosSM[nMax];
851  Double_t MeanPosErrSM[nMax];
852  Double_t WidthSM[nMax];
853  Double_t WidthErrSM[nMax];
854  Double_t Npi0SM[nMax];
855  Double_t Npi0ErrSM[nMax];
856 
857  tree->Branch("CellMeanSM", CellMeanSM, TString::Format("CellMeanSM[%i]/D", nMax));
858  tree->Branch("CellRMSSM", CellRMSSM, TString::Format("CellRMSSM[%i]/D", nMax));
859  tree->Branch("ClusterMeanSM", ClusterMeanSM, TString::Format("ClusterMeanSM[%i]/D", nMax));
860  tree->Branch("ClusterTotSM", ClusterTotSM, TString::Format("ClusterTotSM[%i]/D", nMax));
861  tree->Branch("ClusterRMSSM", ClusterRMSSM, TString::Format("ClusterRMSSM[%i]/D", nMax));
862  tree->Branch("EtotalMeanSM", EtotalMeanSM, TString::Format("EtotalMeanSM[%i]/D", nMax));
863  tree->Branch("EtotalRMSSM", EtotalRMSSM, TString::Format("EtotalRMSSM[%i]/D", nMax));
864  tree->Branch("CellPerClusterMeanSM", CellPerClusterMeanSM, TString::Format("CellPerClusterMeanSM[%i]/D", nMax));
865  tree->Branch("CellPerClusterRMSSM", CellPerClusterRMSSM, TString::Format("CellPerClusterRMSSM[%i]/D", nMax));
866  tree->Branch("ECell1MeanSM", ECell1MeanSM, TString::Format("ECell1MeanSM[%i]/D", nMax));
867  tree->Branch("ECell1RMSSM", ECell1RMSSM, TString::Format("ECell1RMSSM[%i]/D", nMax));
868 
869  tree->Branch("MeanPosSM", MeanPosSM, TString::Format("MeanPosSM[%i]/D", nMax));
870  tree->Branch("MeanPosErrSM", MeanPosErrSM, TString::Format("MeanPosErrSM[%i]/D", nMax));
871  tree->Branch("WidthSM", WidthSM, TString::Format("WidthSM[%i]/D", nMax));
872  tree->Branch("WidthErrSM", WidthErrSM, TString::Format("WidthErrSM[%i]/D", nMax));
873  tree->Branch("Npi0SM", Npi0SM, TString::Format("Npi0SM[%i]/D", nMax));
874  tree->Branch("Npi0ErrSM", Npi0ErrSM, TString::Format("Npi0ErrSM[%i]/D", nMax));
875 
876  TF1* fitMass = new TF1("fitMass", pi0massP2, 100, 250, 6);
877  fitMass->SetParName(0, "A");
878  fitMass->SetParName(1, "m_{0}");
879  fitMass->SetParName(2, "sigma");
880  fitMass->SetParName(3, "a_{0}");
881  fitMass->SetParName(4, "a_{1}");
882  fitMass->SetParName(5, "a_{2}");
883  fitMass->SetParLimits(0, 1.e-5, 1.e5);
884  fitMass->SetParLimits(1, 0.11, 0.16);
885  fitMass->SetParLimits(2, 0.001, 0.06);
886 
887  TList* outputList;
888 
889  TH1F* fhNEvents;
890  TH1F* fhE;
891  TH1F* fhNClusters = 0x0;
892  TH1F* fhECharged;
893  TH1F* fhNCells = 0x0;
894 
895  TH1F* NCells[n];
896  TH1F* NClusters[n];
897  TH2F* NCellsPerCluster[n];
898  TH1F* E[n];
899 
900  TH2F* fhIM;
901  TH2F* fhIMDCAL;
902  TH1F* fhMggEMCAL;
903  TH1F* fhMggDCAL;
904  TH2F* IM[n];
905  TH1F* MggSM[n];
906 
907  TPMERegexp r("_\\w+");
908  TIter next(TriggersList);
909  int ret = 0;
910 
911  while(TObject *obj = next()){
912  fTrigger= TString(obj->GetName());
913  TString namefile = QAPATH + period + "_" + pass + fTrigger(r).Data() + "_" + RunId + "_data.txt";
914  ofstream QAData(namefile, ios::app); // write checks at the end
915 
916  NTClusters=0;
917  NTClustersRMS=0;
918  NCClusters=0;
919  NCClustersRMS=0;
920  NMatchClustersP=0;
921  NMatchClustersPRMS=0;
922  CellMean=0;
923  CellRMS=0;
924  ClusterMean=0;
925  ClusterRMS=0;
926  EtotalMean=0;
927  EtotalRMS=0;
928  CellPerClusterMean=0;
929  CellPerClusterRMS=0;
930 
931  // EMCal
932  Npi0EMCAL=0;
933  Npi0EMCALErr=0;
934  MeanPosEMCAL=0;
935  MeanPosEMCALErr=0;
936  WidthEMCAL=0;
937  WidthEMCALErr=0;
938  Chi2NdfPi0EMCAL=0;
939  NggEMCAL=0;
940  NggEMCALErr=0;
941  SignifEMCAL=0;
942  SignifEMCALErr=0;
943 
944  // DCal
945  Npi0DCAL=0;
946  Npi0DCALErr=0;
947  MeanPosDCAL=0;
948  MeanPosDCALErr=0;
949  WidthDCAL=0;
950  WidthDCALErr=0;
951  Chi2NdfPi0DCAL=0;
952  NggDCAL=0;
953  NggDCALErr=0;
954  SignifDCAL=0;
955  SignifDCALErr=0;
956 
957  memset (CellMeanSM, 0, sizeof (Double_t) * nMax);
958  memset (CellRMSSM, 0, sizeof (Double_t) * nMax);
959  memset (ClusterMeanSM, 0, sizeof (Double_t) * nMax);
960  memset (ClusterTotSM, 0, sizeof (Double_t) * nMax);
961  memset (ClusterRMSSM, 0, sizeof (Double_t) * nMax);
962  memset (EtotalMeanSM, 0, sizeof (Double_t) * nMax);
963  memset (EtotalRMSSM, 0, sizeof (Double_t) * nMax);
964  memset (CellPerClusterMeanSM, 0, sizeof (Double_t) * nMax);
965  memset (CellPerClusterRMSSM, 0, sizeof (Double_t) * nMax);
966  memset (ECell1MeanSM, 0, sizeof (Double_t) * nMax);
967  memset (ECell1RMSSM, 0, sizeof (Double_t) * nMax);
968 
969  memset (MeanPosSM, 0, sizeof (Double_t) * nMax);
970  memset (MeanPosErrSM, 0, sizeof (Double_t) * nMax);
971  memset (WidthSM, 0, sizeof (Double_t) * nMax);
972  memset (WidthErrSM, 0, sizeof (Double_t) * nMax);
973  memset (Npi0SM, 0, sizeof (Double_t) * nMax);
974  memset (Npi0ErrSM, 0, sizeof (Double_t) * nMax);
975 
976  TString dirname;
977  if(!fTrigger.Contains("QA")){
978  dirname = "CaloQA_";
979  }
980  dirname += fTrigger;
981 
982  Bool_t dirok = f->cd(dirname);
983  if(!dirok){
984  Error(__FUNCTION__, Form("No input directory %s", dirname.Data()));
985  tree->Fill(); ftree->cd();
986  tree->Write();
987  ret= -1;
988  continue;
989  }
990 
991  outputList = (TList*)gDirectory->Get(dirname);
992  if(!outputList){
993  Error(__FUNCTION__, Form("No input list! %s", dirname.Data()));
994  tree->Fill(); ftree->cd();
995  tree->Write();
996  ret=-2;
997  continue;;
998  }
999  outputList->SetOwner();
1000 
1001  // number of events
1002  fhNEvents =(TH1F *)outputList->FindObject("hNEvents");
1003  if(!fhNEvents){
1004  Error(__FUNCTION__, Form("NEvent histogram not found for trigger %s", fTrigger.Data()));
1005  tree->Fill(); ftree->cd();
1006  tree->Write(); ret=-3;
1007  continue;
1008  }
1009 
1010  Nevent=fhNEvents->GetEntries();
1011 
1012  if(Nevent==0){
1013  Error(__FUNCTION__, Form("No event in trigger %s", fTrigger.Data()));
1014  tree->Fill(); ftree->cd();
1015  tree->Write();
1016  ret=-4;
1017  continue;
1018  }
1019  if(Nevent<20){
1020  Error(__FUNCTION__, Form("Less than 20 events in trigger %s", fTrigger.Data()));
1021  tree->Fill(); ftree->cd();
1022  tree->Write();
1023  ret=-5;
1024  continue;
1025  }
1026 
1027  TH1F* hE =(TH1F *)outputList->FindObject("EMCAL_hE");
1028  if(!hE || (!hE->GetEntries())){
1029  Error(__FUNCTION__, Form("hE Histogram not found or it is empty for trigger %s, EMCal was not in this run?", fTrigger.Data())); EtotalMean=0;
1030  tree->Fill(); ftree->cd();
1031  tree->Write();
1032  ret=-4;
1033  continue;
1034  }
1035 
1036  // first do clusters trending
1037  fhE = (TH1F *)outputList->FindObject(fCalorimeter+"_hE");
1038  fhECharged = (TH1F *)outputList->FindObject(fCalorimeter+"_hECharged");
1039 
1040  Double_t energy = 0. ;
1041 
1042  for(Int_t ibin = fhE->FindBin(0.6) ; ibin <fhE->FindBin(50.) ; ibin++){
1043  energy+=fhE->GetBinCenter(ibin)*fhE->GetBinContent(ibin);
1044  }
1045 
1046  if(fhE->Integral(fhE->FindBin(0.6), fhE->FindBin(50.))==0){
1047  Error(__FUNCTION__, Form("Not enough events"));
1048  tree->Fill(); ftree->cd();
1049  tree->Write();
1050  ret=-6;
1051  continue;
1052  }
1053 
1054  EtotalMean=energy/fhE->Integral(fhE->FindBin(0.6), fhE->FindBin(50.)) ;
1055  EtotalRMS=fhE->GetMeanError();
1056 
1057  TString nameNCell = Form("%s_hNCells_Mod", fCalorimeter.Data());
1058  TString nameNCluster = Form("%s_hNClusters_Mod", fCalorimeter.Data());
1059  TString nameE = Form("%s_hE_Mod", fCalorimeter.Data());
1060  TH2F* hNCellsMod= (TH2F*)outputList->FindObject(nameNCell.Data());
1061  TH2F* hNClusterMod=(TH2F*)outputList->FindObject(nameNCluster.Data());
1062  TH2F* hEMod=(TH2F*)outputList->FindObject(nameE.Data());
1063 
1064  if(!hNCellsMod || !hNClusterMod || !hEMod){
1065  Error(__FUNCTION__, "A requiered histogram was not found (the imput QAresult.root might be too old)!");
1066  tree->Fill(); ftree->cd();
1067  tree->Write();
1068  ret=-7;
1069  continue;
1070  }
1071 
1072  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1073  // EM addition (Nov. 2016) to observe the PAR Time issue
1074  TString trig_title = "";
1075  if(fTrigger.Contains("EMC"))
1076  trig_title += "EMC";
1077  else
1078  trig_title += "MB";
1079 
1080  // BC/4 = 0
1081  TCanvas* c_TimeSM_BC0 = new TCanvas("TimeSM_BC0", "Cell Time (for BC/4=0) for each SM", 1200, 1200);
1082  c_TimeSM_BC0->SetFillColor(0);
1083  c_TimeSM_BC0->SetBorderSize(0);
1084  c_TimeSM_BC0->SetFrameBorderMode(0);
1085  Bool_t mod3_BC0=0;
1086  if(n%3)
1087  mod3_BC0=1;
1088  c_TimeSM_BC0->Divide(3, n/3+mod3_BC0);
1089 
1090  TH2F* hTimeSM_BC0 = (TH2F*)outputList->FindObject("EMCAL_hTimePerSM_BC0");
1091  TH1F* hTimeSM_BC0_proj[n];
1092  TString TimeSM_BC0_projName;
1093 
1094  if(hTimeSM_BC0){
1095  cout << "---------------------" << endl;
1096  cout << " BC/4=0 | Run " << RunId << endl;
1097  for(Int_t ism = 0 ; ism < n ; ism++){
1098  c_TimeSM_BC0->cd(ism+1);
1099  gPad->SetTopMargin(0.1);
1100  gPad->SetBottomMargin(0.15);
1101  // gPad->SetLogy();
1102 
1103  TimeSM_BC0_projName = Form("SM_%d", ism);
1104  hTimeSM_BC0_proj[ism] = (TH1F *)hTimeSM_BC0->ProjectionX(TimeSM_BC0_projName.Data(), ism+1, ism+1, "") ;
1105  hTimeSM_BC0_proj[ism]->Draw("hist");
1106 
1107  hTimeSM_BC0_proj[ism]->SetTitle(Form("#it{t}_{cell} for super module %i (BC/4=0) %s", ism, trig_title.Data()));
1108  hTimeSM_BC0_proj[ism]->SetTitleSize(0.1);
1109  hTimeSM_BC0_proj[ism]->SetXTitle("#it{t}_{cell} (ns)");
1110  hTimeSM_BC0_proj[ism]->SetYTitle("Nb of entries");
1111  hTimeSM_BC0_proj[ism]->SetStats(kTRUE);
1112 
1113  hTimeSM_BC0_proj[ism]->GetXaxis()->SetRangeUser(400, 900);
1114  hTimeSM_BC0_proj[ism]->GetXaxis()->SetLabelSize(0.05);
1115  hTimeSM_BC0_proj[ism]->GetXaxis()->SetLabelOffset(0.015);
1116  hTimeSM_BC0_proj[ism]->GetXaxis()->SetTitleSize(0.075);
1117  hTimeSM_BC0_proj[ism]->GetXaxis()->SetTitleOffset(0.9);
1118 
1119  hTimeSM_BC0_proj[ism]->GetYaxis()->SetLabelSize(0.05);
1120  hTimeSM_BC0_proj[ism]->GetYaxis()->SetTitleSize(0.06);
1121  hTimeSM_BC0_proj[ism]->GetYaxis()->SetTitleOffset(0.6);
1122  }
1123 
1124  TString outfilename_BC0 = QAPATH + "TimeSM_BC0" + fTrigger(r) + ".pdf";
1125  TString outfilename2_BC0 = QAPATH + "TimeSM_BC0" + fTrigger(r) + ".png";
1126  if(SavePlots==2)
1127  c_TimeSM_BC0->SaveAs(outfilename_BC0);
1128  if(SavePlots)
1129  c_TimeSM_BC0->SaveAs(outfilename2_BC0);
1130  }
1131 
1132  // BC/4 = 1
1133  TCanvas* c_TimeSM_BC1 = new TCanvas("TimeSM_BC1", "Cell Time (for BC/4=1) for each SM", 1200, 1200);
1134  c_TimeSM_BC1->SetFillColor(0);
1135  c_TimeSM_BC1->SetBorderSize(0);
1136  c_TimeSM_BC1->SetFrameBorderMode(0);
1137  Bool_t mod3_BC1=0;
1138  if(n%3)
1139  mod3_BC1=1;
1140  c_TimeSM_BC1->Divide(3, n/3+mod3_BC1);
1141 
1142  TH2F* hTimeSM_BC1 = (TH2F*)outputList->FindObject("EMCAL_hTimePerSM_BC1");
1143  TH1F* hTimeSM_BC1_proj[n];
1144  TString TimeSM_BC1_projName;
1145 
1146  if(hTimeSM_BC1){
1147  cout << "---------------------" << endl;
1148  cout << " BC/4=1 | Run " << RunId << endl;
1149  for(Int_t ism = 0 ; ism < n ; ism++){
1150  c_TimeSM_BC1->cd(ism+1);
1151  gPad->SetTopMargin(0.1);
1152  gPad->SetBottomMargin(0.15);
1153  // gPad->SetLogy();
1154 
1155  TimeSM_BC1_projName = Form("SM_%d", ism);
1156  hTimeSM_BC1_proj[ism] = (TH1F *)hTimeSM_BC1->ProjectionX(TimeSM_BC1_projName.Data(), ism+1, ism+1, "") ;
1157  hTimeSM_BC1_proj[ism]->Draw("hist");
1158 
1159  hTimeSM_BC1_proj[ism]->SetTitle(Form("#it{t}_{cell} for super module %i (BC/4=1) %s", ism, trig_title.Data()));
1160  hTimeSM_BC1_proj[ism]->SetTitleSize(0.1);
1161  hTimeSM_BC1_proj[ism]->SetXTitle("#it{t}_{cell} (ns)");
1162  hTimeSM_BC1_proj[ism]->SetYTitle("Nb of entries");
1163  hTimeSM_BC1_proj[ism]->SetStats(kTRUE);
1164 
1165  hTimeSM_BC1_proj[ism]->GetXaxis()->SetRangeUser(400, 900);
1166  hTimeSM_BC1_proj[ism]->GetXaxis()->SetLabelSize(0.05);
1167  hTimeSM_BC1_proj[ism]->GetXaxis()->SetLabelOffset(0.015);
1168  hTimeSM_BC1_proj[ism]->GetXaxis()->SetTitleSize(0.075);
1169  hTimeSM_BC1_proj[ism]->GetXaxis()->SetTitleOffset(0.9);
1170 
1171  hTimeSM_BC1_proj[ism]->GetYaxis()->SetLabelSize(0.05);
1172  hTimeSM_BC1_proj[ism]->GetYaxis()->SetTitleSize(0.06);
1173  hTimeSM_BC1_proj[ism]->GetYaxis()->SetTitleOffset(0.6);
1174  }
1175 
1176  TString outfilename_BC1 = QAPATH + "TimeSM_BC1" + fTrigger(r) + ".pdf";
1177  TString outfilename2_BC1 = QAPATH + "TimeSM_BC1" + fTrigger(r) + ".png";
1178  if(SavePlots==2)
1179  c_TimeSM_BC1->SaveAs(outfilename_BC1);
1180  if(SavePlots)
1181  c_TimeSM_BC1->SaveAs(outfilename2_BC1);
1182  }
1183 
1184  // BC/4 = 2
1185  TCanvas* c_TimeSM_BC2 = new TCanvas("TimeSM_BC2", "Cell Time (for BC/4=2) for each SM", 1200, 1200);
1186  c_TimeSM_BC2->SetFillColor(0);
1187  c_TimeSM_BC2->SetBorderSize(0);
1188  c_TimeSM_BC2->SetFrameBorderMode(0);
1189  Bool_t mod3_BC2=0;
1190  if(n%3)
1191  mod3_BC2=1;
1192  c_TimeSM_BC2->Divide(3, n/3+mod3_BC2);
1193 
1194  TH2F* hTimeSM_BC2 = (TH2F*)outputList->FindObject("EMCAL_hTimePerSM_BC2");
1195  TH1F* hTimeSM_BC2_proj[n];
1196  TString TimeSM_BC2_projName;
1197 
1198  if(hTimeSM_BC2){
1199  cout << "---------------------" << endl;
1200  cout << " BC/4=2 | Run " << RunId << endl;
1201  for(Int_t ism = 0 ; ism < n ; ism++){
1202  c_TimeSM_BC2->cd(ism+1);
1203  gPad->SetTopMargin(0.1);
1204  gPad->SetBottomMargin(0.15);
1205  // gPad->SetLogy();
1206 
1207  TimeSM_BC2_projName = Form("SM_%d", ism);
1208  hTimeSM_BC2_proj[ism] = (TH1F *)hTimeSM_BC2->ProjectionX(TimeSM_BC2_projName.Data(), ism+1, ism+1, "") ;
1209  hTimeSM_BC2_proj[ism]->Draw("hist");
1210 
1211  hTimeSM_BC2_proj[ism]->SetTitle(Form("#it{t}_{cell} for super module %i (BC/4=2) %s", ism, trig_title.Data()));
1212  hTimeSM_BC2_proj[ism]->SetTitleSize(0.1);
1213  hTimeSM_BC2_proj[ism]->SetXTitle("#it{t}_{cell} (ns)");
1214  hTimeSM_BC2_proj[ism]->SetYTitle("Nb of entries");
1215  hTimeSM_BC2_proj[ism]->SetStats(kTRUE);
1216 
1217  hTimeSM_BC2_proj[ism]->GetXaxis()->SetRangeUser(400, 900);
1218  hTimeSM_BC2_proj[ism]->GetXaxis()->SetLabelSize(0.05);
1219  hTimeSM_BC2_proj[ism]->GetXaxis()->SetLabelOffset(0.015);
1220  hTimeSM_BC2_proj[ism]->GetXaxis()->SetTitleSize(0.075);
1221  hTimeSM_BC2_proj[ism]->GetXaxis()->SetTitleOffset(0.9);
1222 
1223  hTimeSM_BC2_proj[ism]->GetYaxis()->SetLabelSize(0.05);
1224  hTimeSM_BC2_proj[ism]->GetYaxis()->SetTitleSize(0.06);
1225  hTimeSM_BC2_proj[ism]->GetYaxis()->SetTitleOffset(0.6);
1226  }
1227 
1228  TString outfilename_BC2 = QAPATH + "TimeSM_BC2" + fTrigger(r) + ".pdf";
1229  TString outfilename2_BC2 = QAPATH + "TimeSM_BC2" + fTrigger(r) + ".png";
1230  if(SavePlots==2)
1231  c_TimeSM_BC2->SaveAs(outfilename_BC2);
1232  if(SavePlots)
1233  c_TimeSM_BC2->SaveAs(outfilename2_BC2);
1234  }
1235 
1236  // BC/4 = 3
1237  TCanvas* c_TimeSM_BC3 = new TCanvas("TimeSM_BC3", "Cell Time (for BC/4=3) for each SM", 1200, 1200);
1238  c_TimeSM_BC3->SetFillColor(0);
1239  c_TimeSM_BC3->SetBorderSize(0);
1240  c_TimeSM_BC3->SetFrameBorderMode(0);
1241  Bool_t mod3_BC3=0;
1242  if(n%3)
1243  mod3_BC3=1;
1244  c_TimeSM_BC3->Divide(3, n/3+mod3_BC3);
1245 
1246  TH2F* hTimeSM_BC3 = (TH2F*)outputList->FindObject("EMCAL_hTimePerSM_BC3");
1247  TH1F* hTimeSM_BC3_proj[n];
1248  TString TimeSM_BC3_projName;
1249 
1250  if(hTimeSM_BC3){
1251  cout << "---------------------" << endl;
1252  cout << " BC/4=3 | Run " << RunId << endl;
1253  for(Int_t ism = 0 ; ism < n ; ism++){
1254  c_TimeSM_BC3->cd(ism+1);
1255  gPad->SetTopMargin(0.1);
1256  gPad->SetBottomMargin(0.15);
1257  // gPad->SetLogy();
1258 
1259  TimeSM_BC3_projName = Form("SM_%d", ism);
1260  hTimeSM_BC3_proj[ism] = (TH1F *)hTimeSM_BC3->ProjectionX(TimeSM_BC3_projName.Data(), ism+1, ism+1, "") ;
1261  hTimeSM_BC3_proj[ism]->Draw("hist");
1262 
1263  hTimeSM_BC3_proj[ism]->SetTitle(Form("#it{t}_{cell} for super module %i (BC/4=3) %s", ism, trig_title.Data()));
1264  hTimeSM_BC3_proj[ism]->SetTitleSize(0.1);
1265  hTimeSM_BC3_proj[ism]->SetXTitle("#it{t}_{cell} (ns)");
1266  hTimeSM_BC3_proj[ism]->SetYTitle("Nb of entries");
1267  hTimeSM_BC3_proj[ism]->SetStats(kTRUE);
1268 
1269  hTimeSM_BC3_proj[ism]->GetXaxis()->SetRangeUser(400, 900);
1270  hTimeSM_BC3_proj[ism]->GetXaxis()->SetLabelSize(0.05);
1271  hTimeSM_BC3_proj[ism]->GetXaxis()->SetLabelOffset(0.015);
1272  hTimeSM_BC3_proj[ism]->GetXaxis()->SetTitleSize(0.075);
1273  hTimeSM_BC3_proj[ism]->GetXaxis()->SetTitleOffset(0.9);
1274 
1275  hTimeSM_BC3_proj[ism]->GetYaxis()->SetLabelSize(0.05);
1276  hTimeSM_BC3_proj[ism]->GetYaxis()->SetTitleSize(0.06);
1277  hTimeSM_BC3_proj[ism]->GetYaxis()->SetTitleOffset(0.6);
1278  }
1279 
1280  TString outfilename_BC3 = QAPATH + "TimeSM_BC3" + fTrigger(r) + ".pdf";
1281  TString outfilename2_BC3 = QAPATH + "TimeSM_BC3" + fTrigger(r) + ".png";
1282  if(SavePlots==2)
1283  c_TimeSM_BC3->SaveAs(outfilename_BC3);
1284  if(SavePlots)
1285  c_TimeSM_BC3->SaveAs(outfilename2_BC3);
1286  }
1287 
1288  //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1289 
1290  TCanvas* c1 = new TCanvas("Pi0InvMassSM", "Pi0 Invariant Mass for each SM", 600, 600);
1291  gStyle->SetOptStat(1); // MG modif
1292  gStyle->SetOptFit(1); // MG modif
1293  c1->SetFillColor(0);
1294  c1->SetBorderSize(0);
1295  c1->SetFrameBorderMode(0);
1296  // c1->SetOptStat(1); // MG modif
1297 
1298  Bool_t mod3=0;
1299  if(n%3)
1300  mod3=1;
1301  c1->Divide(3, n/3+mod3);
1302 
1303  // per sm trending
1304  TString nameNCellPerCluster;
1305  for(Int_t ism = 0 ; ism < n ; ism++){
1306  cout << "#########################"<< endl;
1307  cout << " Super Module " << ism << " Run " << RunId << endl;
1308 
1309  // first do clusters trending
1310  nameNCellPerCluster = Form("%s_hNCellsPerCluster_Mod%d", fCalorimeter.Data(), ism);
1311  NCellsPerCluster[ism] = (TH2F*)outputList->FindObject(nameNCellPerCluster.Data());
1312 
1313  if(! (TH2F*)outputList->FindObject(nameNCellPerCluster.Data()) ){
1314  Error(__FUNCTION__, Form("NCellsPerCluster histogram not found for super module %i", ism));
1315  ret=-8;
1316  continue;
1317  }
1318 
1319  NCellsPerCluster[ism] = (TH2F*)outputList->FindObject(nameNCellPerCluster.Data());
1320  NCells[ism] = (TH1F*)hNCellsMod->ProjectionX(Form("NCells%d", ism), ism+1, ism+1, "");
1321  NClusters[ism] = (TH1F*)hNClusterMod->ProjectionX(Form("NClusters%d", ism), ism+1, ism+1, "");
1322  E[ism] = (TH1F*)hEMod->ProjectionX(Form("E%d", ism), ism+1, ism+1, "");
1323 
1324  CellMeanSM[ism]=NCells[ism]->GetMean();
1325  CellRMSSM[ism]=NCells[ism]->GetMeanError();
1326 
1327  Int_t binmax;
1328  binmax = NClusters[ism]->GetNbinsX();
1329  for(Int_t ibin = 1 ;ibin <=binmax ;ibin++){
1330  ClusterTotSM[ism]+=(NClusters[ism]->GetBinCenter(ibin))*(NClusters[ism]->GetBinContent(ibin));
1331  }
1332  if(NClusters[ism]->GetEntries()>0)
1333  ClusterTotSM[ism]= ClusterTotSM[ism]/NClusters[ism]->GetEntries();
1334 
1335  // ClusterTotSM[ism]=NClusters[ism]->GetMean();
1336  ClusterMeanSM[ism]=NClusters[ism]->GetMean();
1337  ClusterRMSSM[ism]=NClusters[ism]->GetMeanError();
1338  CellPerClusterMeanSM[ism]=NCellsPerCluster[ism]->GetMean(2);
1339  CellPerClusterRMSSM[ism]=NCellsPerCluster[ism]->GetMeanError(2);
1340 
1341  ECell1MeanSM[ism] =NCellsPerCluster[ism]->ProjectionX("", 2, 50, "")->Integral(5, 50)/(Nevent);
1342  ECell1RMSSM[ism] =NCellsPerCluster[ism]->ProjectionX("", 2, 50, "")->GetMeanError();
1343 
1344  Double_t energySM = 0. ;
1345  for(Int_t ibin = E[ism]->FindBin(0.6) ; ibin <E[ism]->FindBin(50.) ; ibin++){
1346  energySM+=E[ism]->GetBinCenter(ibin)*(E[ism]->GetBinContent(ibin));
1347  }
1348 
1349  if(E[ism]->Integral(E[ism]->FindBin(0.6), E[ism]->FindBin(50.))==0){
1350  Error(__FUNCTION__, Form("Energy: Not enough events/SM"));
1351  continue;
1352  }
1353 
1354  EtotalMeanSM[ism]=energySM/(E[ism]->Integral(E[ism]->FindBin(0.6), E[ism]->FindBin(50.)));
1355  EtotalRMSSM[ism]=E[ism]->GetMeanError();
1356 
1357  if(ism==0){
1358  fhNCells = (TH1F*)NCells[ism]->Clone("NCells");
1359  fhNClusters = (TH1F*)NClusters[ism]->Clone("NClusters");
1360  }
1361  else{
1362  fhNCells->Add(NCells[ism], 1);
1363  fhNClusters->Add(NClusters[ism], 1);
1364  }
1365 
1366  // Pi0
1367  c1->cd(ism+1);
1368  TString namePair = Form("%s_hIM_Mod%d", fCalorimeter.Data(), ism);
1369  IM[ism] = (TH2F*)outputList->FindObject(namePair.Data());
1370  IM[ism]->Sumw2();
1371 
1372  TString projname = Form("SM_%d", ism);
1373  MggSM[ism] = (TH1F *)IM[ism]->ProjectionY(projname.Data(), 2, 20, "") ; // MG modif
1374 
1375  if(MggSM[ism]->GetEntries()>100){
1376  fitMass->SetParameter(0, MggSM[ism]->GetBinContent(MggSM[ism]->GetMaximumBin()));
1377  fitMass->SetParameter(1, mPDG);
1378  fitMass->SetParameter(2, 15.);
1379  fitMass->SetParameter(3, 0.);
1380  fitMass->SetParameter(4, MggSM[ism]->GetBinContent(MggSM[ism]->FindBin(0.11)));
1381  fitMass->SetParameter(5, MggSM[ism]->GetBinContent(MggSM[ism]->FindBin(0.20)));
1382 
1383  if(MggSM[ism]->GetEntries()<1000)
1384  MggSM[ism]->Rebin(5);
1385  else
1386  MggSM[ism]->Rebin();
1387 
1388  MggSM[ism]->Fit("fitMass", "WL R +", "", 0.1, 0.30);
1389  MggSM[ism]->SetTitle(Form("Pi0 Mass for super module %i", ism));
1390  MggSM[ism]->SetTitleSize(0.1);
1391  MggSM[ism]->SetXTitle("Pi0 Mass");
1392  MggSM[ism]->SetYTitle("Nb of entries");
1393  MggSM[ism]->GetXaxis()->SetLabelSize(0.05);
1394  MggSM[ism]->GetXaxis()->SetTitleSize(0.07);
1395  MggSM[ism]->GetXaxis()->SetTitleOffset(0.68);
1396  MggSM[ism]->GetYaxis()->SetLabelSize(0.05);
1397  MggSM[ism]->GetYaxis()->SetTitleSize(0.06);
1398  MggSM[ism]->GetYaxis()->SetTitleOffset(0.78);
1399 
1400  MeanPosSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(1)*1000;
1401  MeanPosErrSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParError(1)*1000;
1402  WidthSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(2)*1000;
1403  WidthErrSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParError(2)*1000;
1404  Npi0SM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(0)*(MggSM[ism]->GetFunction("fitMass")->GetParameter(2))*TMath::Sqrt(2*TMath::Pi())/(Nevent*MggSM[ism]->GetBinWidth(1));
1405  Npi0ErrSM[ism] = TMath::Sqrt((MggSM[ism]->GetFunction("fitMass")->GetParError(0)/MggSM[ism]->GetFunction("fitMass")->GetParameter(0))*(MggSM[ism]->GetFunction("fitMass")->GetParError(0)/MggSM[ism]->GetFunction("fitMass")->GetParameter(0))
1406  +(MggSM[ism]->GetFunction("fitMass")->GetParError(2)/MggSM[ism]->GetFunction("fitMass")->GetParameter(2))*(MggSM[ism]->GetFunction("fitMass")->GetParError(2)/MggSM[ism]->GetFunction("fitMass")->GetParameter(2)));
1407  Npi0ErrSM[ism] = 0.;
1408 
1409  } // end if enough events for Pi0 fit and trending
1410  else
1411  Info(__FUNCTION__, Form("Not enough events for Pi0 fit and trending for super module %i", ism));
1412  } // per SM loop
1413 
1414  // Now Pi0 global trending in EMCAL
1415  TCanvas* c2 = new TCanvas("Pi0InvMassEMCAL", "Pi0 Invariant Mass in EMCal", 600, 600);
1416  c2->SetFillColor(0);
1417  c2->SetBorderSize(0);
1418  c2->SetFrameBorderMode(0);
1419  // c2->SetOptStat(1); // MG modif
1420 
1421  fhIM = (TH2F *)outputList->FindObject(fCalorimeter+"_hIM");
1422  fhIM->Sumw2();
1423  fhMggEMCAL = (TH1F *)fhIM->ProjectionY("MggEMCAL", 2, 20, "") ; // to modify projection range
1424 
1425  if(fhMggEMCAL->GetEntries()==0){
1426  Error(__FUNCTION__, "The Pi0 histogram in EMCal is empty !"); tree->Fill();
1427  ret=-8;
1428  continue;
1429  }
1430 
1431  fitMass->SetParameter(0, 4500);
1432  fitMass->SetParameter(1, mPDG);
1433  fitMass->SetParameter(2, 0.01);
1434  fitMass->SetParameter(3, 0.);
1435  fitMass->SetParameter(4, fhMggEMCAL->GetBinContent(fhMggEMCAL->FindBin(0.11)));
1436  fitMass->SetParameter(5, fhMggEMCAL->GetBinContent(fhMggEMCAL->FindBin(0.20)));
1437 
1438  if(fhMggEMCAL->GetEntries()<5000)
1439  fhMggEMCAL->Rebin(5); // MG modif
1440  else
1441  fhMggEMCAL->Rebin();
1442 
1443  fhMggEMCAL->Fit("fitMass", "L R +", "", 0.05, 0.20);
1444 
1445  fhMggEMCAL->SetTitle("Pi0 Mass in EMCal");
1446  fhMggEMCAL->SetTitleSize(0.1);
1447  fhMggEMCAL->SetXTitle("Pi0 Mass");
1448  fhMggEMCAL->SetYTitle("Nb of entries");
1449  fhMggEMCAL->GetXaxis()->SetLabelSize(0.03);
1450  fhMggEMCAL->GetXaxis()->SetTitleSize(0.03);
1451  fhMggEMCAL->GetXaxis()->SetTitleOffset(1.3);
1452  fhMggEMCAL->GetYaxis()->SetLabelSize(0.03);
1453  fhMggEMCAL->GetYaxis()->SetTitleSize(0.03);
1454  fhMggEMCAL->GetYaxis()->SetTitleOffset(1.3);
1455 
1456  MeanPosEMCAL = fhMggEMCAL->GetFunction("fitMass")->GetParameter(1)*1000;
1457  MeanPosEMCALErr = fhMggEMCAL->GetFunction("fitMass")->GetParError(1)*1000;
1458  WidthEMCAL = fhMggEMCAL->GetFunction("fitMass")->GetParameter(2)*1000;
1459  WidthEMCALErr = fhMggEMCAL->GetFunction("fitMass")->GetParError(2)*1000;
1460  Chi2NdfPi0EMCAL = fhMggEMCAL->GetFunction("fitMass")->GetChisquare()/fhMggEMCAL->GetFunction("fitMass")->GetNDF();
1461  Npi0EMCAL = fhMggEMCAL->GetFunction("fitMass")->GetParameter(0)*fhMggEMCAL->GetFunction("fitMass")->GetParameter(2)*TMath::Sqrt(2*TMath::Pi())/(Nevent*fhMggEMCAL->GetBinWidth(10));
1462  Npi0EMCALErr = TMath::Sqrt((fhMggEMCAL->GetFunction("fitMass")->GetParError(0)/fhMggEMCAL->GetFunction("fitMass")->GetParameter(0))*(fhMggEMCAL->GetFunction("fitMass")->GetParError(0)/fhMggEMCAL->GetFunction("fitMass")->GetParameter(0))+(WidthEMCALErr/WidthEMCAL)*(WidthEMCALErr/WidthEMCAL));
1463  Npi0EMCALErr = 0.;
1464  NggEMCAL = fhMggEMCAL->GetFunction("fitMass")->Integral(0.11, 0.16)/(Nevent*fhMggEMCAL->GetBinWidth(10));
1465  NggEMCALErr = fhMggEMCAL->GetFunction("fitMass")->IntegralError(0.11, 0.16)/(fhMggEMCAL->Integral()*fhMggEMCAL->GetBinWidth(10));
1466  SignifEMCAL = Npi0EMCAL/NggEMCAL;
1467  SignifEMCALErr = TMath::Sqrt((Npi0EMCALErr/Npi0EMCAL*(Npi0EMCALErr/Npi0EMCAL)+(NggEMCALErr/NggEMCAL*(NggEMCALErr/NggEMCAL))));
1468  SignifEMCALErr = SignifEMCAL*SignifEMCALErr;
1469 
1470  TCanvas* c3 = new TCanvas("Pi0InvMassDCAL", "Pi0 Invariant Mass in DCal", 600, 600);
1471  c3->SetFillColor(0);
1472  c3->SetBorderSize(0);
1473  c3->SetFrameBorderMode(0);
1474 
1475  fhIMDCAL = (TH2F *)outputList->FindObject(fCalorimeter+"_hIMDCAL");
1476  if(fhIMDCAL){
1477  // Now Pi0 global trending in DCal
1478  // TCanvas* c3 = new TCanvas("Pi0InvMassDCAL", "Pi0 Invariant Mass in DCal", 600, 600);
1479  // c3->SetFillColor(0);
1480  // c3->SetBorderSize(0);
1481  // c3->SetFrameBorderMode(0);
1482  // c2->SetOptStat(1); // MG modif
1483 
1484  // fhIMDCAL = (TH2F *)outputList->FindObject(fCalorimeter+"_hIMDCAL");
1485  // if(!fhIMDCAL) continue;
1486  fhIMDCAL->Sumw2();
1487 
1488  fhMggDCAL = (TH1F *)fhIMDCAL->ProjectionY("MggDCAL", 2, 20, "") ; // to modify projection range
1489  if(fhMggDCAL->GetEntries()==0){
1490  Error(__FUNCTION__, "The Pi0 histogram in DCal is empty !"); tree->Fill();
1491  ret=-8;
1492  continue;
1493  }
1494 
1495  fitMass->SetParameter(0, 4500);
1496  fitMass->SetParameter(1, mPDG);
1497  fitMass->SetParameter(2, 0.01);
1498  fitMass->SetParameter(3, 0.);
1499  fitMass->SetParameter(4, fhMggDCAL->GetBinContent(fhMggDCAL->FindBin(0.11)));
1500  fitMass->SetParameter(5, fhMggDCAL->GetBinContent(fhMggDCAL->FindBin(0.20)));
1501 
1502  if(fhMggDCAL->GetEntries()<5000)
1503  fhMggDCAL->Rebin(5); // MG modif
1504  else
1505  fhMggDCAL->Rebin();
1506 
1507  fhMggDCAL->Fit("fitMass", "L R +", "", 0.05, 0.20);
1508 
1509  fhMggDCAL->SetTitle("Pi0 Mass in DCal");
1510  fhMggDCAL->SetTitleSize(0.1);
1511  fhMggDCAL->SetXTitle("Pi0 Mass");
1512  fhMggDCAL->SetYTitle("Nb of entries");
1513  fhMggDCAL->GetXaxis()->SetLabelSize(0.03);
1514  fhMggDCAL->GetXaxis()->SetTitleSize(0.03);
1515  fhMggDCAL->GetXaxis()->SetTitleOffset(1.3);
1516  fhMggDCAL->GetYaxis()->SetLabelSize(0.03);
1517  fhMggDCAL->GetYaxis()->SetTitleSize(0.03);
1518  fhMggDCAL->GetYaxis()->SetTitleOffset(1.3);
1519 
1520  MeanPosDCAL = fhMggDCAL->GetFunction("fitMass")->GetParameter(1)*1000;
1521  MeanPosDCALErr = fhMggDCAL->GetFunction("fitMass")->GetParError(1)*1000;
1522  WidthDCAL = fhMggDCAL->GetFunction("fitMass")->GetParameter(2)*1000;
1523  WidthDCALErr = fhMggDCAL->GetFunction("fitMass")->GetParError(2)*1000;
1524  Chi2NdfPi0DCAL = fhMggDCAL->GetFunction("fitMass")->GetChisquare()/fhMggDCAL->GetFunction("fitMass")->GetNDF();
1525  Npi0DCAL = fhMggDCAL->GetFunction("fitMass")->GetParameter(0)*fhMggDCAL->GetFunction("fitMass")->GetParameter(2)*TMath::Sqrt(2*TMath::Pi())/(Nevent*fhMggDCAL->GetBinWidth(10));
1526  Npi0DCALErr = TMath::Sqrt((fhMggDCAL->GetFunction("fitMass")->GetParError(0)/fhMggDCAL->GetFunction("fitMass")->GetParameter(0))*(fhMggDCAL->GetFunction("fitMass")->GetParError(0)/fhMggDCAL->GetFunction("fitMass")->GetParameter(0))+(WidthDCALErr/WidthDCAL)*(WidthDCALErr/WidthDCAL));
1527  Npi0DCALErr = 0.;
1528  NggDCAL = fhMggDCAL->GetFunction("fitMass")->Integral(0.11, 0.16)/(Nevent*fhMggDCAL->GetBinWidth(10));
1529  NggDCALErr = fhMggDCAL->GetFunction("fitMass")->IntegralError(0.11, 0.16)/(fhMggDCAL->Integral()*fhMggDCAL->GetBinWidth(10));
1530  SignifDCAL = Npi0DCAL/NggDCAL;
1531  SignifDCALErr = TMath::Sqrt((Npi0DCALErr/Npi0DCAL*(Npi0DCALErr/Npi0DCAL)+(NggDCALErr/NggDCAL*(NggDCALErr/NggDCAL))));
1532  SignifDCALErr = SignifDCAL*SignifDCALErr;
1533  }
1534 
1535  cout<<"******************"<<endl;
1536  // end of global trending
1537 
1538  // cout<<"********************************** Total number of clusters **********************************"<<endl;
1539  NTClusters=fhE->IntegralAndError(fhE->GetXaxis()->FindBin(0.), fhE->GetXaxis()->FindBin(200.), NTClustersRMS);
1540  NCClusters=fhECharged->IntegralAndError(fhECharged->GetXaxis()->FindBin(0.), fhECharged->GetXaxis()->FindBin(50.), NCClustersRMS);
1541 
1542  if(NTClusters!=0)
1543  NMatchClustersP = NCClusters/NTClusters;
1544  else
1545  NMatchClustersP=0;
1546  // cout<<" here is searched value "<<NMatchClustersP<<" "<<endl;
1547 
1548  if(NCClusters!=0)
1549  NMatchClustersPRMS = NMatchClustersP* TMath::Sqrt(TMath::Power(NCClustersRMS/NCClusters, 2)+TMath::Power(NTClustersRMS/NTClusters, 2));
1550  else
1551  NMatchClustersPRMS=1;
1552 
1553  ClusterMean=fhNClusters->GetMean();
1554  ClusterRMS=fhNClusters->GetMeanError();
1555  CellMean=fhNCells->GetMean();
1556  CellRMS=fhNCells->GetMeanError();
1557 
1558  tree->Fill();
1559 
1560  TString outfilename = QAPATH + "Pi0InvMassEMCAL" + fTrigger(r) + ".pdf";
1561  TString outfilename2 = QAPATH + "Pi0InvMassEMCAL" + fTrigger(r) + ".png";
1562  if(SavePlots==2)
1563  c2->SaveAs(outfilename);
1564  if(SavePlots)
1565  c2->SaveAs(outfilename2);
1566 
1567  outfilename = QAPATH + "Pi0InvMassSM" + fTrigger(r) + ".pdf";
1568  outfilename2 = QAPATH + "Pi0InvMassSM" + fTrigger(r) + ".png";
1569  if(SavePlots==2)
1570  c1->SaveAs(outfilename);
1571  if(SavePlots)
1572  c1->SaveAs(outfilename2);
1573  // cout<<"********************************** "<<n<<" **********************************"<<endl; // n = number of SMs
1574 
1575  if(fhIMDCAL){
1576  TString outfilename4 = QAPATH + "Pi0InvMassDCAL" + fTrigger(r) + ".pdf";
1577  TString outfilename5= QAPATH + "Pi0InvMassDCAL" + fTrigger(r) + ".png";
1578  if(SavePlots==2)
1579  c3->SaveAs(outfilename4);
1580  if(SavePlots)
1581  c3->SaveAs(outfilename5);
1582  }
1583 
1584  fout->cd();
1585  fout->Cd(Form("%s/%s/%ld/%s/%s", period.Data(), pass.Data(), RunId, "RunLevelQA", fTrigger.Data()));
1586 
1587  // Cell time distribution for each BC and SM
1588  if(hTimeSM_BC0) c_TimeSM_BC0->Write();
1589  if(hTimeSM_BC1) c_TimeSM_BC1->Write();
1590  if(hTimeSM_BC2) c_TimeSM_BC2->Write();
1591  if(hTimeSM_BC3) c_TimeSM_BC3->Write();
1592 
1593  c2->Write();
1594  c1->Write();
1595  if(fhIMDCAL)
1596  c3->Write();
1597 
1598  delete c1;
1599  delete c2;
1600  if(fhIMDCAL) delete c3;
1601  if(hTimeSM_BC0) delete c_TimeSM_BC0;
1602  if(hTimeSM_BC1) delete c_TimeSM_BC1;
1603  if(hTimeSM_BC2) delete c_TimeSM_BC2;
1604  if(hTimeSM_BC3) delete c_TimeSM_BC3;
1605  if(outputList) outputList->Delete();
1606 
1607  QAData << RunId << " " << Nevent
1608  << "\n";
1609 
1610  QAData.close();
1611  }
1612 
1613  ftree->cd();
1614  tree->Write();
1615  ftree->Close();
1616 
1617  return ret;
1618 
1619 }
1620 
1621 //-------------------------------------------------------------------------
1622 TH2F* FormatRunHisto(TH2F* aHisto, const char* title, const char* YTitle){
1623 
1624  if(!aHisto){
1625  Error(__FUNCTION__, Form("The histogram with title \"%s\" was not found!", title));
1626  return new TH2F();
1627  }
1628 
1629  aHisto->SetStats(kFALSE);
1630  aHisto->SetTitle(title);
1631  aHisto->SetStats(kFALSE);
1632  aHisto->SetYTitle(YTitle);
1633  aHisto->GetYaxis()->SetTitleOffset(1.2);
1634  aHisto->GetYaxis()->SetLabelSize(0.03);
1635  aHisto->GetZaxis()->SetLabelSize(0.02);
1636 
1637  return aHisto;
1638 
1639 }
1640 
1641 //--------------------------------------------------------------------------------------------------------------
1642 TH2F* HistoPerMod(TH2F* hTmpPerMod, const char* title){
1643 
1644  if(!hTmpPerMod){
1645  Error(__FUNCTION__, Form("The histogram with title \"%s\" was not found!", title));
1646  return new TH2F();
1647  }
1648 
1649  hTmpPerMod->SetStats(kFALSE);
1650  hTmpPerMod->SetTitle(title);
1651  hTmpPerMod->SetTitleSize(0.1);
1652  hTmpPerMod->GetXaxis()->SetTitleOffset(1.1);
1653  hTmpPerMod->GetXaxis()->SetTitleSize(0.05);
1654  hTmpPerMod->GetXaxis()->SetLabelSize(0.06);
1655  hTmpPerMod->GetYaxis()->SetTitleOffset(1.1);
1656  hTmpPerMod->GetYaxis()->SetTitleSize(0.05);
1657  hTmpPerMod->GetYaxis()->SetLabelSize(0.06);
1658  hTmpPerMod->GetZaxis()->SetLabelSize(0.04);
1659 
1660  return hTmpPerMod;
1661 
1662 }
1663 
1664 //---------------------------------------------------------------------------------------------------
1665 TH2* AutoZoom(TH2* H, Option_t* aType, Int_t EntryMin){
1666 
1667  Int_t shiftX = (Int_t)(H->GetNbinsX()/30.);
1668  Int_t shiftY = (Int_t)(H->GetNbinsY()/30.);
1669 
1670  TString opt = aType;
1671  opt.ToLower();
1672 
1673  int minX = 0;
1674  int maxX = H->GetNbinsX();
1675  int New_minX = minX;
1676  int New_maxX = maxX;
1677 
1678  int minY = 0;
1679  int maxY = H->GetNbinsY();
1680  int New_minY = minY;
1681  int New_maxY = maxY;
1682 
1683  if(opt.Contains("all")) opt = TString("minx, maxx, miny, maxy");
1684 
1685  if(opt.Contains("maxx")){
1686  for(New_maxX = maxX;New_maxX >=minX; New_maxX--){
1687  Stat_t c = 0;
1688  for(int i_y = maxY; i_y >= minY;i_y--){
1689  c = H->GetBinContent(New_maxX, i_y);
1690  if(c>EntryMin)
1691  break;
1692  }
1693  if(c>EntryMin)
1694  break;
1695  }
1696  }
1697 
1698  if(opt.Contains("maxy")){
1699  for(New_maxY = maxY;New_maxY >=minY;New_maxY--){
1700  Stat_t c = 0;
1701  for(int i_x=maxX; i_x>=minX;i_x--){
1702  c = H->GetBinContent(i_x, New_maxY );
1703  if(c>EntryMin)
1704  break;
1705  }
1706  if(c>EntryMin)
1707  break;
1708  }
1709  }
1710 
1711  if(opt.Contains("minx")){
1712  for(New_minX = minX;New_minX <=maxX; New_minX++){
1713  Stat_t c = 0;
1714  for(int i_y = minY; i_y <= maxY;i_y++){
1715  c = H->GetBinContent(New_minX, i_y);
1716  if(c>EntryMin)
1717  break;
1718  }
1719  if(c>EntryMin)
1720  break;
1721  }
1722  }
1723 
1724  if(opt.Contains("miny")){
1725  for(New_minY = minY;New_minY <=maxY;New_minY++){
1726  Stat_t c = 0;
1727  for(int i_x=minX; i_x<=maxX;i_x++){
1728  c = H->GetBinContent(i_x, New_minY );
1729  if(c>EntryMin)
1730  break;
1731  }
1732  if(c>EntryMin)
1733  break;
1734  }
1735  }
1736 
1737  if(New_maxX!=-1 && New_maxY!=-1)
1738  H->GetXaxis()->SetRange(New_minX - shiftX, New_maxX + shiftX);
1739  if(New_maxX!=-1 && New_maxY!=-1)
1740  H->GetYaxis()->SetRange(New_minY - shiftY, New_maxY + shiftY);
1741 
1742  return H;
1743 
1744 }
1745 
1746 //----------------------------------------------------------------------------------------------------
1747 int FindNumberOfSM(TFile* f, TString fTrigger, TString period){
1748 
1749  TString direct;
1750  if(!fTrigger.Contains("QA")){
1751  direct = "CaloQA_";
1752  }
1753  direct += fTrigger;
1754 
1755  Int_t nSMt=-1;
1756  Int_t year = 2000 + TString(period(3, 2)).Atoi();
1757  if ( year == 2010 ) { nSMt=4; }
1758  else if( year == 2011 || year == 2012 ){ nSMt=10; }
1759  else if( year == 2013 || year == 2014 ){ nSMt=12; }
1760  else { nSMt=20; }
1761 
1762  TList *outputList;
1763  Bool_t dirok = f->cd(direct);
1764 
1765  // ////
1766  // cout << "==============================================================" << endl;
1767  // cout << "direct = " << direct << endl;
1768  // cout << "dirok = " << dirok << endl;
1769  // cout << "==============================================================" << endl;
1770  // ////
1771 
1772  if(!dirok)
1773  Error(__FUNCTION__, Form("No input directory %s, the number SMs will be returned based on the year!", direct.Data()));
1774  else
1775  outputList = (TList*)gDirectory->Get(direct);
1776 
1777  if(!outputList)
1778  Error(__FUNCTION__, "No input list, the number SMs will be returned based on the year! ");
1779  else{
1780  outputList->SetOwner();
1781  TH2F* hNSM =(TH2F *)outputList->FindObject("EMCAL_hE_Mod");
1782  if(!hNSM || (!hNSM->GetEntries()))
1783  Error(__FUNCTION__, "hNSM Histogram not found or it is empty, the number SMs will be returned based on the year!");
1784  else
1785  nSMt = hNSM->GetYaxis()->GetBinUpEdge(hNSM->FindLastBinAbove(0, 2));
1786  }
1787 
1788  if(outputList)
1789  outputList->Delete();
1790 
1791  return nSMt;
1792 
1793 }
TH2F * FormatRunHisto(TH2F *aHisto, const char *title, const char *YTitle="")
const char * filename
Definition: TestFCM.C:1
double Double_t
Definition: External.C:58
TString QAPATH
Definition: External.C:236
const char * title
Definition: MakeQAPdf.C:26
energy
Definition: HFPtSpectrum.C:44
TSystem * gSystem
Int_t DrawOccupancy(Long_t run, TString period, TString pass, TString fTrigger, TString system, TFile *f, TFile *fout, AliEMCALGeometry *geom, Int_t SavePlots)
TCanvas * c
Definition: TestFitELoss.C:172
int FindNumberOfSM(TFile *f, TString fTrigger, TString period)
Double_t pi0massP1(Double_t *x, Double_t *par)
Double_t pi0massP2(Double_t *x, Double_t *par)
Int_t DrawRun(Long_t run, TString period, TString pass, TString fTrigger, TFile *f, TFile *fout, Int_t SavePlots, Int_t nSM, Bool_t kFilter)
TString QAPATHF
int Int_t
Definition: External.C:63
Int_t TrendingEMCALTree(Long_t RunId, TString fCalorimeter, TString system, TString period, TString pass, const int n, TList *TriggersList, TFile *f, TFile *fout, Int_t SavePlots)
Definition: External.C:228
TH2 * AutoZoom(TH2 *H, Option_t *aType="all", Int_t EntryMin=0)
void set_plot_style()
TH2F * HistoPerMod(TH2F *name, const char *title)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: External.C:220
Double_t fitE(Double_t *x, Double_t *par)
const char Option_t
Definition: External.C:48
bool Bool_t
Definition: External.C:53
TFile * fout
input train file
int CreateEMCALRunQA(const char *filename, TString RunId, TString period, TString pass, Int_t SavePlots=0, Bool_t filter=0, TString fTrigger="", TString system="", TString fCalorimeter="EMCAL")