AliPhysics  vAN-20151012 (2287573)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
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 
38 Int_t DrawOccupancy(Long_t run, TString period, TString pass, TString fTrigger, TString system, TFile* f, TFile* fout, AliEMCALGeometry* geom, Int_t SavePlots);
39 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);
40 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);
41 
42 TH2F* FormatRunHisto(TH2F* aHisto,const char* title,const char* YTitle="");
43 TH2F* HistoPerMod(TH2F* name,const char* title);
44 TH2* AutoZoom(TH2* H,Option_t* aType="all", Int_t EntryMin=0);
45 int FindNumberOfSM(TFile* f, TString fTrigger,TString period);
46 
47 TString QAPATH;
48 TString QAPATHF = "./";
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 //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
65 Double_t pi0massP2(Double_t *x, Double_t *par)
66 {
67  Double_t gaus;
68 
69  if (par[2] != 0.) gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) / (2*par[2]*par[2]) );
70 
71  else gaus = 99999999.;
72 
73  Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
74 
75  return gaus+back;
76 
77 }
78 
79 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
80 Double_t pi0massP1(Double_t *x, Double_t *par)
81 {
82  Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) / (2*par[2]*par[2]) );
83 
84  Double_t back = par[3] + par[4]*x[0];
85 
86  return gaus+back;
87 
88 }
89 
90 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
91 Double_t fitE(Double_t *x, Double_t *par)
92 {
93 
94  Double_t levy;
95 
96  levy = par[0] * TMath::Exp( -par[1]/x[0]) * TMath::Power(x[0], -par[2]) ;
97 
98  return levy;
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 
105  QAPATH = TString(gSystem->Getenv("QAPATH"));
106  if(QAPATH.IsNull()) QAPATH = QAPATHF;
107  if(! QAPATH.BeginsWith("./")) { QAPATH = QAPATH + RunId + "/";}
108 
109  AliLog::SetGlobalLogLevel(AliLog::kError);
110  TFile *f = new TFile(filename);
111  AliLog::SetGlobalLogLevel(AliLog::kInfo);
112 
113  if (f->IsZombie()) {Error(__FUNCTION__,Form("Error openning the input file %s",filename)); return -1;}
114 
115  TList* TriggersList = new TList();
116 
117  if (fTrigger=="")
118  {
119  TPMERegexp name_re("CaloQA_\\w+");
120  TObjLink* link = f->GetListOfKeys()->FirstLink();
121 
122  while (link)
123  {
124  TString name = link->GetObject()->GetName();
125  if (name_re.Match(name))
126  {
127  TriggersList->Add(link->GetObject());
128  if(TString(filename).Contains("barrel") && ! name.Contains("default")) TriggersList->Remove(link->GetObject());
129  if(TString(filename).Contains("outer") && ! name.Contains("EMC")) TriggersList->Remove(link->GetObject());
130  }
131  link = link->Next();
132  }
133  } else {TriggersList->Add(new TObjString(fTrigger.Data()));}
134 
135  if(!TriggersList->GetEntries()) {Error(__FUNCTION__,"No trigger found!"); return -2;}
136 
137  int nSM = FindNumberOfSM(f,((TObjString*)TriggersList->Last())->GetString(),period);
138  if (nSM<0) {Error(__FUNCTION__,"Could not find the number of super modules!"); return -3;}
139  Info(__FUNCTION__,Form("%i super modules were discuvered",nSM));
140  TString GeomName;
141  if (nSM <= 4) { nSM=4; GeomName = "EMCAL_FIRSTYEARv1";}
142  else if (nSM <= 10) { nSM=10; GeomName = "EMCAL_COMPLETEv1";}
143  else if (nSM <= 12) { nSM=12; GeomName = "EMCAL_COMPLETE12SMv1";}
144  else {nSM = 20; GeomName = "EMCAL_COMPLETE12SMv1_DCAL_8SM";}
145 
146  AliEMCALGeometry *geom = new AliEMCALGeometry(GeomName.Data(),"EMCAL");
147  Info(__FUNCTION__,Form("Using %i super modules and the Geometry %s",nSM,GeomName.Data()));
148 
149  TFile *fout = new TFile(TString( QAPATH + period+"_"+pass + fTrigger+"_"+ (Long_t)RunId.Atoi() +"_QAplots.root").Data(),"RECREATE");
150 
151  if((system.IsNull()) && (period.EndsWith("h"))) {system = "PbPb";}
152 
153  Int_t ret=0;
154  TIter next(TriggersList);
155  while (TObject *obj = next())
156  {
157  fTrigger= TString(obj->GetName());
158  ret -= DrawOccupancy(RunId.Atoi(),period,pass,fTrigger,system,f,fout,geom,SavePlots);
159  ret -= DrawRun(RunId.Atoi(),period,pass,fTrigger,f,fout,SavePlots,nSM,filter);
160  }
161  ret-= TrendingEMCALTree(RunId.Atoi(),fCalorimeter,system,period,pass,nSM,TriggersList,f,fout,SavePlots);
162 
163  f->Close();
164  fout->Close();
165  delete f;
166  delete geom;
167 
168  return ret;
169 }
170 
171 //----------------------------------------------------------------------------------------------------------------------------------------------------------
172 Int_t DrawOccupancy(Long_t run , TString period, TString pass, TString fTrigger,TString system, TFile* f,TFile* fout, AliEMCALGeometry* geom, Int_t SavePlots)
173 {
174 
175  set_plot_style();
176  gStyle->SetOptStat(0);
177  TH1::AddDirectory(kFALSE);
178  TH2D *hEnergyMapReal = new TH2D("hEnergyMapReal","",96,-48,48,120,-0,120);
179  TH2D *hOccupancyMapReal = new TH2D("hOccupancyMapReal","",96,-48,48,120,-0,120);
180  hEnergyMapReal->SetXTitle("eta (bin)");
181  hEnergyMapReal->SetYTitle("phi (bin)");
182  hEnergyMapReal->SetZTitle("E(GeV)/event");
183  hEnergyMapReal->GetYaxis()->SetTitleOffset(1.2);
184  hEnergyMapReal->GetZaxis()->SetLabelSize(0.02);
185  hEnergyMapReal->GetZaxis()->SetTitleOffset(1.36);
186 
187  hOccupancyMapReal->SetXTitle("eta (bin)");
188  hOccupancyMapReal->SetYTitle("phi (bin)");
189  hOccupancyMapReal->GetYaxis()->SetTitleOffset(1.2);
190  hOccupancyMapReal->GetZaxis()->SetLabelSize(0.02);
191 
192  Int_t nSupMod, nModule, nIphi, nIeta;
193  Int_t iphi, ieta;
194  Int_t realbineta=0;
195  Int_t realbinphi=0;
196 
197  //NO MASK
198  Int_t mask[1] = {2222222};
199 
200  TH2F *hCellAmplitude;
201  TH1F *hNEvents;
202  Int_t Events;
203  Int_t n=0;
204 
205  TString direct;
206  if(!fTrigger.Contains("QA")) {
207  direct = "CaloQA_";
208  }
209  direct += fTrigger;
210  Bool_t dirok = f->cd(direct);
211  if (!dirok) { Error(__FUNCTION__,Form("No input drectory %s",direct.Data())); return -1;}
212  TList *outputList = (TList*)gDirectory->Get(direct);
213  if(!outputList){ Error(__FUNCTION__,"No input list! "); return -1;}
214  outputList->SetOwner();
215 
216  fout->mkdir(Form("%s/%s/%ld/%s/%s",period.Data(),pass.Data(),run,"RunLevelQA",fTrigger.Data()));
217  fout->cd();
218  fout->Cd(Form("%s/%s/%ld/%s/%s",period.Data(),pass.Data(),run,"RunLevelQA",fTrigger.Data()));
219 
220  hNEvents =(TH1F *)outputList->FindObject("hNEvents");
221  if(!hNEvents){ Error(__FUNCTION__,Form("hNEvent histogram not found for trigger %s ",fTrigger.Data())); return -2;}
222  Events = (Int_t)hNEvents->GetEntries();
223  if(Events==0){ Error(__FUNCTION__,Form("No event in trigger %s",fTrigger.Data())); return -3;}
224 
225  TH1F* hE =(TH1F *)outputList->FindObject("EMCAL_hE");
226  if (!hE || (!hE->GetEntries())) { Error(__FUNCTION__,Form("hE Histogram not found or it is empty for trigger %s, EMCAL was not in this run?",fTrigger.Data())); return -4;}
227 
228  Double_t Eth=1;
229  if(system=="PbPb"){
230  Eth = 5.;
231  if (fTrigger.Contains("EMC")) Eth=20.;
232  }
233  if(system=="pp"){
234  Eth = 1.;
235  if (fTrigger.Contains("EMC")) Eth=5.;
236  }
237 
238  hCellAmplitude =(TH2F *)outputList->FindObject("EMCAL_hAmpId");
239 
240  for(Int_t i = 0; i < geom->GetNCells() ; i++){
241  Double_t Esum = 0;
242  Double_t Nsum = 0;
243 
244  for (Int_t j = 1; j <= hCellAmplitude->GetNbinsX(); j++)
245  {
246  Double_t E = hCellAmplitude->GetXaxis()->GetBinCenter(j);
247  Double_t N = hCellAmplitude->GetBinContent(j, i+1);
248 
249  if (E < 0.3) continue;
250 
251  if (E <= Eth) {
252  Esum += E*N;
253  Nsum += N;
254  }
255  }
256 
257  Int_t absId = i;
258  if(n!=0) {if(mask[n]<=mask[n-1]) Warning(__FUNCTION__,"The list of bad cells is not sorted !!");}
259  if(i==mask[n]){n++ ; continue; } // skip bad cells
260 
261  geom->GetCellIndex(absId, nSupMod, nModule, nIphi, nIeta);
262  geom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
263 
264  realbinphi = 120-(nSupMod/2)*24 -iphi -1; //
265  if (nSupMod%2==0) realbineta= 48-ieta -1;
266  if (nSupMod%2==1) realbineta= -ieta -1;
267 
268  hEnergyMapReal->Fill(realbineta,realbinphi,Esum/(Double_t)Events);
269  // hOccupancyMapReal->Fill(realbineta,realbinphi,Nsum/(Double_t)Events);
270  }
271 
272  cout <<" Run: " << run << " trigger: " << fTrigger << " N_events: "<<Events<<endl;
273 
274  TPMERegexp r("_\\w+");
275  TString Energy; Energy = QAPATH + "MapEnergy" + fTrigger(r) + ".pdf";
276  TString Energy2; Energy2 = QAPATH + "MapEnergy" + fTrigger(r) + ".png";
277  TString Entries; Entries = QAPATH + "MapEntries" + fTrigger(r) + ".pdf";
278  TString Entries2; Entries2 = QAPATH + "MapEntries" + fTrigger(r) + ".png";
279 
280  TCanvas *c1 = new TCanvas("Energymap","Energy Map",600,600);
281  c1->SetFillColor(0);
282  c1->SetGrid();
283  c1->SetRightMargin(0.14);
284  TString title = "run ";
285  title += run ;
286  if(fTrigger.Contains("EMC")) { title += " EMC ";} else {title += " MB ";}
287  title += " Summed energy map";
288  hEnergyMapReal->SetTitle(title);
289  AutoZoom(hEnergyMapReal,"miny")->DrawCopy("colz");
290  if(SavePlots==2) c1->SaveAs(Energy);
291  if(SavePlots) c1->SaveAs(Energy2);
292  c1->Write();
293  delete c1;
294 
295 
296  // TCanvas *c2 = new TCanvas("Occupancy","Occupancy Map",600,600);
297  // c2->SetFillColor(0);
298  // c2->SetGrid();
299  // c2->SetRightMargin(0.14);
300  // TString title2 = "run ";
301  // title2 += run ;
302  // if(fTrigger.Contains("EMC")) { title2 += " EMC ";} else { title2 += " MB ";}
303  // title2 += " Occupancy map";
304  // hOccupancyMapReal->SetTitle(title2);
305  // AutoZoom(hOccupancyMapReal,"miny")->DrawCopy("colz");
306  // if(SavePlots==2) c2->SaveAs(Entries);
307  // if(SavePlots) c2->SaveAs(Entries2);
308  // c2->Write();
309  // delete c2;
310 
311  if (outputList) {outputList->Delete();}
312 
313  delete hEnergyMapReal;
314  delete hOccupancyMapReal;
315 
316  return 0;
317 }
318 
319 //-----------------------------------------------------------------------------------------------------------------------
320 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)
321 {
322 
323  TString direct;
324  if(!fTrigger.Contains("QA")) {
325  direct = "CaloQA_";
326  }
327  direct += fTrigger;
328 
329  f->cd(direct);
330  if(!direct) { Error(__FUNCTION__,Form("No input directory %s",direct.Data())); return -1;}
331  TList *outputList = (TList*)gDirectory->Get(direct);
332  if(!outputList){ Error(__FUNCTION__,Form("No input list! %s",direct.Data())); return -2;}
333  outputList->SetOwner();
334  if (kFilter)
335  {
336  fout->mkdir(Form("%s/%s/%ld/%s",period.Data(),pass.Data(),run,fTrigger.Data()));
337  fout->cd();
338  fout->Cd(Form("%s/%s/%ld/%s",period.Data(),pass.Data(),run,fTrigger.Data()));
339  outputList->Write();
340  }
341  fout->cd();
342  fout->Cd(Form("%s/%s/%ld/%s/%s",period.Data(),pass.Data(),run,"RunLevelQA",fTrigger.Data()));
343 
344 
345  set_plot_style();
346  gStyle->SetPalette(1);
347  // gStyle->SetOptStat(1);
348  TH1::AddDirectory(kFALSE);
349  TString outfilename;
350  TString outfilename2;
351  const char* legend="";
352  TPMERegexp r("_\\w+");
353 
354  if (fTrigger.Contains("EMC")){ legend = Form(" Run %d EMC ",(int)run);}
355  else legend = Form(" Run %d MB ",(int)run);
356 
357  TH1F* hNEvents =(TH1F *)outputList->FindObject("hNEvents");
358  if(!hNEvents){ Error(__FUNCTION__,Form("hNEvent histogram not found for trigger %s ",fTrigger.Data())); return -3;}
359  Int_t Events = (Int_t)hNEvents->GetEntries();
360  if(Events==0){ Error(__FUNCTION__,Form("No event in trigger %s",fTrigger.Data())); return -4 ;}
361 
362  TH1F* hE =(TH1F *)outputList->FindObject("EMCAL_hE");
363  if (!hE || (!hE->GetEntries())) { Error(__FUNCTION__,Form("hE Histogram not found or it is empty for trigger %s, EMCAL was not in this run?",fTrigger.Data())); return -4;}
364 
366  TCanvas* c00 = new TCanvas("Occupancy map", "Occupancy", 600, 600);
367  c00->SetLogz();
368  c00->SetFillColor(0);
369  c00->SetBorderSize(0);
370  c00->SetFrameBorderMode(0);
371 
372  TH2F* hOccupancyMapReal =(TH2F *)outputList->FindObject("EMCAL_hEtaPhi");
373  if(!hOccupancyMapReal) { Error(__FUNCTION__,Form("EMCAL_hEtaPhi: Histogram for trigger %s not found!",fTrigger.Data())); return -5;}
374  FormatRunHisto(hOccupancyMapReal,Form("Occupancy %s",legend),"phi (bin)");
375  hOccupancyMapReal->SetXTitle("eta (bin)");
376 
377  AutoZoom(hOccupancyMapReal,"all")->DrawCopy("colz");
378  outfilename = QAPATH + "OccupancyMap" + fTrigger(r) + ".pdf" ;
379  outfilename2 = QAPATH + "OccupancyMap" + fTrigger(r) + ".png" ;
380 
381  if(SavePlots==2) c00->SaveAs(outfilename);
382  if(SavePlots) c00->SaveAs(outfilename2);
383  c00->Write();
384  delete c00;
385 
386  TCanvas* c0 = new TCanvas("Grid Cell", "Occupancy from Grid Cells", 600, 600);
387  c0->SetLogz();
388  c0->SetFillColor(0);
389  c0->SetBorderSize(0);
390  c0->SetFrameBorderMode(0);
391 
392  TH2F* hGridCells =(TH2F *)outputList->FindObject("EMCAL_hGridCells");
393  if(!hGridCells) { Error(__FUNCTION__,Form("EMCAL_hGridCells: Histogram for trigger %s not found!",fTrigger.Data())); return -5;}
394  FormatRunHisto(hGridCells,Form("Occupancy from Grid Cells %s",legend),"phi row columns");
395  hGridCells->SetXTitle("eta row columns");
396 
397  AutoZoom(hGridCells,"all")->DrawCopy("colz");
398  outfilename = QAPATH + "GridCells" + fTrigger(r) + ".pdf" ;
399  outfilename2 = QAPATH + "GridCells" + fTrigger(r) + ".png" ;
400 
401  if(SavePlots==2) c0->SaveAs(outfilename);
402  if(SavePlots) c0->SaveAs(outfilename2);
403  c0->Write();
404  delete c0;
405 
406 
407  //
408 
409  TCanvas* c1 = new TCanvas("TimeVsE", "Cluster Time Vs Energy", 600, 600);
410  c1->SetLogz();
411  c1->SetFillColor(0);
412  c1->SetBorderSize(0);
413  c1->SetFrameBorderMode(0);
414 
415  TH2F* hClusterTimeEnergy =(TH2F *)outputList->FindObject("EMCAL_hClusterTimeEnergy");
416  if(!hClusterTimeEnergy) { Error(__FUNCTION__,Form("EMCAL_hClusterTimeEnergy: Histogram for trigger %s not found!",fTrigger.Data())); return -5;}
417  FormatRunHisto(hClusterTimeEnergy,Form("Time Vs Energy%s",legend),"EMCAL ToF(ns)");
418 
419  AutoZoom(hClusterTimeEnergy,"maxx")->DrawCopy("colz");
420  outfilename = QAPATH + "TimeRun" + fTrigger(r) + ".pdf" ;
421  outfilename2 = QAPATH + "TimeRun" + fTrigger(r) + ".png" ;
422 
423  if(SavePlots==2) c1->SaveAs(outfilename);
424  if(SavePlots) c1->SaveAs(outfilename2);
425  c1->Write();
426  delete c1;
427 
428  TCanvas * c2 = new TCanvas("ClusterVsTrack ","Correlation calo Mult Vs Track Multiplicity", 600, 600);
429  c2->SetLogz();
430  c2->SetFillColor(0);
431  c2->SetBorderSize(0);
432  c2->SetFrameBorderMode(0);
433 
434  TH2F* hClusterVsTrack =(TH2F *)outputList->FindObject("EMCAL_hCaloTrackMNClusters");
435  FormatRunHisto(hClusterVsTrack,Form("N cluster Vs N track%s",legend));
436 
437  AutoZoom(hClusterVsTrack,"maxx,maxy",1)->DrawCopy("colz");
438  outfilename = QAPATH + "CaloTrackMult" + fTrigger(r) + ".pdf";
439  outfilename2 = QAPATH + "CaloTrackMult" + fTrigger(r) + ".png";
440  if(SavePlots==2) c2->SaveAs(outfilename);
441  if(SavePlots) c2->SaveAs(outfilename2);
442  c2->Write();
443  delete c2;
444 
445  TCanvas* c3 = new TCanvas("ClusterEVsTrack ","Correlation E calo Vs Track Multiplicity", 600, 600);
446  c3->SetLogz();
447  c3->SetFillColor(0);
448  c3->SetBorderSize(0);
449  c3->SetFrameBorderMode(0);
450 
451  TH2F* hClusterEVsTrack =(TH2F*)outputList->FindObject("EMCAL_hCaloTrackMEClusters");
452  FormatRunHisto(hClusterEVsTrack,Form("Sum E cluster Vs N track%s",legend));
453 
454  AutoZoom(hClusterEVsTrack,"maxx,maxy",1)->DrawCopy("colz");
455  outfilename = QAPATH + "ETrackMult" + fTrigger(r) + ".pdf";
456  outfilename2 = QAPATH + "ETrackMult" + fTrigger(r) + ".png";
457  if(SavePlots==2) c3->SaveAs(outfilename);
458  if(SavePlots) c3->SaveAs(outfilename2);
459  c3->Write();
460  delete c3;
461 
462  TCanvas* c4 = new TCanvas("ClusterEVsV0 ","Correlation E calo Vs V0 signal", 600, 600);
463  c4->SetLogz();
464  c4->SetFillColor(0);
465  c4->SetBorderSize(0);
466  c4->SetFrameBorderMode(0);
467 
468  TH2F* hClusterEVsV0S =(TH2F*)outputList->FindObject("EMCAL_hCaloV0SEClusters");
469  FormatRunHisto(hClusterEVsV0S,Form("Sum E cluster Vs V0 signal%s",legend));
470 
471  AutoZoom(hClusterEVsV0S,"maxx,maxy",1)->DrawCopy("colz");
472  outfilename = QAPATH +"EVsV0s" + fTrigger(r) + ".pdf";
473  outfilename2 = QAPATH +"EVsV0s" + fTrigger(r) + ".png";
474  if(SavePlots==2) c4->SaveAs(outfilename);
475  if(SavePlots) c4->SaveAs(outfilename2);
476  c4->Write();
477  delete c4;
478 
479  TCanvas* c5 = new TCanvas("CellsperCluster","Nb of cells per cluster for each SM", 600, 600);
480  c5->SetLogz();
481  c5->SetFillColor(0);
482  c5->SetBorderSize(0);
483  c5->SetFrameBorderMode(0);
484  Bool_t mod3=0; if (nSM%3) mod3=1;
485  c5->Divide(3,(nSM/3)+mod3);
486 
487  for (int ism = 0; ism < nSM; ism++)
488  {
489  c5->cd(ism+1);
490  gPad->SetLogz();
491  if(TString(Form("Nb of cells per cluster%s Mod %d",legend,ism)).Length() > 60) { Error(__FUNCTION__,"Title too long!"); return -6;}
492  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");
493  }
494 
495  outfilename = QAPATH + "CellsperClusterSM" + fTrigger(r) + ".pdf";
496  outfilename2 = QAPATH + "CellsperClusterSM" + fTrigger(r) + ".png";
497  if(SavePlots==2) c5->SaveAs(outfilename);
498  if(SavePlots) c5->SaveAs(outfilename2);
499  c5->Write();
500  delete c5;
501 
502  if (outputList) outputList->Delete();
503 
504  return 0;
505 }
506 
507 //----------------------------------------------------------------------------------------------------------------------------------
508 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)
509 {
510 
511  TString fTrigger="";
512  TString aCalorimeter;
513  if (n<=12) {aCalorimeter = fCalorimeter;} else {aCalorimeter = TString("EMCAL_and_DCAL");}
514  TDatime now;
515 
516  Double_t Nevent=0 ;
517  Double_t xe=0.5;
518  Double_t NTClusters=0; // total number of clusters
519  Double_t NTClustersRMS=0;
520  Double_t NCClusters=0; // number of charged clusters
521  Double_t NCClustersRMS=0;
522  Double_t NMatchClustersP=0;
523  Double_t NMatchClustersPRMS=0;
524 
525  Double_t CellMean=0;
526  Double_t CellRMS=0;
527  Double_t ClusterMean=0;
528  Double_t ClusterRMS=0;
529  Double_t EtotalMean=0;
530  Double_t EtotalRMS=0;
531 
532  Double_t CellPerClusterMean=0; //
533  Double_t CellPerClusterRMS=0; //
534 
535  Double_t mPDG = 134.9766;
536  Double_t Npi0EMCAL=0;
537  Double_t Npi0EMCALErr=0;
538  Double_t MeanPosEMCAL=0;
539  Double_t MeanPosEMCALErr=0;
540  Double_t WidthEMCAL=0;
541  Double_t WidthEMCALErr=0;
542  Double_t Chi2NdfPi0EMCAL=0;
543  Double_t NggEMCAL=0;
544  Double_t NggEMCALErr=0;
545  Double_t SignifEMCAL=0; // !S/(S+B)
546  Double_t SignifEMCALErr=0; // !S/(S+B)
547 
548  Double_t Npi0DCAL=0;
549  Double_t Npi0DCALErr=0;
550  Double_t MeanPosDCAL=0;
551  Double_t MeanPosDCALErr=0;
552  Double_t WidthDCAL=0;
553  Double_t WidthDCALErr=0;
554  Double_t Chi2NdfPi0DCAL=0;
555  Double_t NggDCAL=0;
556  Double_t NggDCALErr=0;
557  Double_t SignifDCAL=0; // !S/(S+B)
558  Double_t SignifDCALErr=0; // !S/(S+B)
559 
560  TFile* ftree = new TFile(Form("%s/trending.root",QAPATH.Data()),"RECREATE");
561 
562  TTree *tree = new TTree("trending","Trending QA Tree");
563  tree->Branch("fDate",&now);
564  tree->Branch("fCalorimeter",&aCalorimeter);
565  tree->Branch("system",&system);
566  tree->Branch("period",&period);
567  tree->Branch("pass",&pass);
568  tree->Branch("fTrigger",&fTrigger);
569  tree->Branch("run",&RunId,"run/I");
570  tree->Branch("xe",&xe,"xe/D");
571 
572  tree->Branch("Nevent",&Nevent,"Nevent/D");
573  tree->Branch("NMatchClustersP",&NMatchClustersP,"NMatchClustersP/D");
574  tree->Branch("NMatchClustersPRMS",&NMatchClustersPRMS,"NMatchClustersPRMS/D");
575  tree->Branch("CellMean",&CellMean,"CellMean/D");
576  tree->Branch("CellRMS",&CellRMS,"CellRMS/D");
577  tree->Branch("ClusterMean",&ClusterMean,"ClusterMean/D");
578  tree->Branch("ClusterRMS",&ClusterRMS,"ClusterRMS/D");
579  tree->Branch("EtotalMean",&EtotalMean,"EtotalMean/D");
580  tree->Branch("EtotalRMS",&EtotalRMS,"EtotalRMS/D");
581 
582  tree->Branch("CellPerClusterMean",&CellPerClusterMean,"CellPerClusterMean/D"); //
583  tree->Branch("CellPerClusterRMS",&CellPerClusterRMS,"CellPerClusterRMS/D"); //
584 
585  tree->Branch("Npi0EMCAL",&Npi0EMCAL,"Npi0EMCAL/D");
586  tree->Branch("Npi0EMCALErr",&Npi0EMCALErr,"Npi0EMCALErr/D");
587  tree->Branch("Npi0DCAL",&Npi0DCAL,"Npi0DCAL/D");
588  tree->Branch("Npi0DCALErr",&Npi0DCALErr,"Npi0DCALErr/D");
589  tree->Branch("MeanPosEMCAL",&MeanPosEMCAL,"MeanPosEMCAL/D");
590  tree->Branch("MeanPosEMCALErr",&MeanPosEMCALErr,"MeanPosEMCALErr/D");
591  tree->Branch("MeanPosDCAL",&MeanPosDCAL,"MeanPosDCAL/D");
592  tree->Branch("MeanPosDCALErr",&MeanPosDCALErr,"MeanPosDCALErr/D");
593  tree->Branch("WidthEMCAL",&WidthEMCAL,"WidthEMCAL/D");
594  tree->Branch("WidthEMCALErr",&WidthEMCALErr,"WidthEMCALErr/D");
595  tree->Branch("WidthDCAL",&WidthDCAL,"WidthDCAL/D");
596  tree->Branch("WidthDCALErr",&WidthDCALErr,"WidthDCALErr/D");
597  tree->Branch("Chi2NdfPi0EMCAL",&Chi2NdfPi0EMCAL,"Chi2NdfPi0EMCAL/D");
598  tree->Branch("Chi2NdfPi0DCAL",&Chi2NdfPi0DCAL,"Chi2NdfPi0DCAL/D");
599  tree->Branch("NggEMCAL",&NggEMCAL,"NggEMCAL/D");
600  tree->Branch("NggEMCALErr",&NggEMCALErr,"NggEMCALErr/D");
601  tree->Branch("NggDCAL",&NggDCAL,"NggDCAL/D");
602  tree->Branch("NggDCALErr",&NggDCALErr,"NggDCALErr/D");
603  tree->Branch("SignifEMCAL",&SignifEMCAL,"SignifEMCAL/D");
604  tree->Branch("SignifDCALErr",&SignifDCALErr,"SignifDCALErr/D");
605  tree->Branch("SignifEMCAL",&SignifEMCAL,"SignifEMCAL/D");
606  tree->Branch("SignifDCALErr",&SignifDCALErr,"SignifDCALErr/D");
607 
608  tree->Branch("nSM",&n,"nSM/I");
609 
610  int nMax = 22;
611  Double_t CellMeanSM[nMax];
612  Double_t CellRMSSM[nMax];
613  Double_t ClusterMeanSM[nMax];
614  Double_t ClusterTotSM[nMax];
615  Double_t ClusterRMSSM[nMax];
616  Double_t EtotalMeanSM[nMax]; //mean total energy deposited per event
617  Double_t EtotalRMSSM[nMax];
618  Double_t CellPerClusterMeanSM[nMax];
619  Double_t CellPerClusterRMSSM[nMax];
620  Double_t ECell1MeanSM[nMax]; //total energy deposited per event without 1 cell clusters
621  Double_t ECell1RMSSM[nMax];
622 
623  Double_t MeanPosSM[nMax];
624  Double_t MeanPosErrSM[nMax];
625  Double_t WidthSM[nMax];
626  Double_t WidthErrSM[nMax];
627  Double_t Npi0SM[nMax];
628  Double_t Npi0ErrSM[nMax];
629 
630  tree->Branch("CellMeanSM",CellMeanSM,TString::Format("CellMeanSM[%i]/D",nMax));
631  tree->Branch("CellRMSSM",CellRMSSM,TString::Format("CellRMSSM[%i]/D",nMax));
632  tree->Branch("ClusterMeanSM",ClusterMeanSM,TString::Format("ClusterMeanSM[%i]/D",nMax));
633  tree->Branch("ClusterTotSM",ClusterTotSM,TString::Format("ClusterTotSM[%i]/D",nMax));
634  tree->Branch("ClusterRMSSM",ClusterRMSSM,TString::Format("ClusterRMSSM[%i]/D",nMax));
635  tree->Branch("EtotalMeanSM",EtotalMeanSM,TString::Format("EtotalMeanSM[%i]/D",nMax));
636  tree->Branch("EtotalRMSSM",EtotalRMSSM,TString::Format("EtotalRMSSM[%i]/D",nMax));
637  tree->Branch("CellPerClusterMeanSM",CellPerClusterMeanSM,TString::Format("CellPerClusterMeanSM[%i]/D",nMax));
638  tree->Branch("CellPerClusterRMSSM",CellPerClusterRMSSM,TString::Format("CellPerClusterRMSSM[%i]/D",nMax));
639  tree->Branch("ECell1MeanSM",ECell1MeanSM,TString::Format("ECell1MeanSM[%i]/D",nMax));
640  tree->Branch("ECell1RMSSM",ECell1RMSSM,TString::Format("ECell1RMSSM[%i]/D",nMax));
641 
642  tree->Branch("MeanPosSM",MeanPosSM,TString::Format("MeanPosSM[%i]/D",nMax));
643  tree->Branch("MeanPosErrSM",MeanPosErrSM,TString::Format("MeanPosErrSM[%i]/D",nMax));
644  tree->Branch("WidthSM",WidthSM,TString::Format("WidthSM[%i]/D",nMax));
645  tree->Branch("WidthErrSM",WidthErrSM,TString::Format("WidthErrSM[%i]/D",nMax));
646  tree->Branch("Npi0SM",Npi0SM,TString::Format("Npi0SM[%i]/D",nMax));
647  tree->Branch("Npi0ErrSM",Npi0ErrSM,TString::Format("Npi0ErrSM[%i]/D",nMax));
648 
649  TF1* fitMass = new TF1("fitMass",pi0massP2,100,250,6);
650  fitMass->SetParName(0,"A");
651  fitMass->SetParName(1,"m_{0}");
652  fitMass->SetParName(2,"sigma");
653  fitMass->SetParName(3,"a_{0}");
654  fitMass->SetParName(4,"a_{1}");
655  fitMass->SetParName(5,"a_{2}");
656  fitMass->SetParLimits(0, 1.e-5,1.e5);
657  fitMass->SetParLimits(1, 0.11, 0.16); //
658  fitMass->SetParLimits(2, 0.001,0.06);
659 
660  TList* outputList;
661 
662  TH1F* fhNEvents;
663  TH1F* fhE;
664  TH1F* fhNClusters = 0x0;
665  TH1F* fhECharged;
666  TH1F* fhNCells = 0x0;
667 
668  TH1F* NCells[n];
669  TH1F* NClusters[n];
670  TH2F* NCellsPerCluster[n];
671  TH1F* E[n];
672 
673  TH2F* fhIM;
674  TH2F* fhIMDCAL;
675  TH1F* fhMggEMCAL;
676  TH1F* fhMggDCAL;
677  TH2F* IM[n];
678  TH1F* MggSM[n];
679 
680  TPMERegexp r("_\\w+");
681  TIter next(TriggersList);
682  int ret = 0;
683  while (TObject *obj = next())
684  {
685  fTrigger= TString(obj->GetName());
686  TString namefile = QAPATH + period + "_" + pass + fTrigger(r).Data() + "_" + RunId + "_data.txt";
687  ofstream QAData(namefile, ios::app); // write checks at the end
688 
689  NTClusters=0;
690  NTClustersRMS=0;
691  NCClusters=0;
692  NCClustersRMS=0;
693  NMatchClustersP=0;
694  NMatchClustersPRMS=0;
695  CellMean=0;
696  CellRMS=0;
697  ClusterMean=0;
698  ClusterRMS=0;
699  EtotalMean=0;
700  EtotalRMS=0;
701  CellPerClusterMean=0; //
702  CellPerClusterRMS=0; //
703 
704 //EMCAL
705  Npi0EMCAL=0;
706  Npi0EMCALErr=0;
707  MeanPosEMCAL=0;
708  MeanPosEMCALErr=0;
709  WidthEMCAL=0;
710  WidthEMCALErr=0;
711  Chi2NdfPi0EMCAL=0;
712  NggEMCAL=0;
713  NggEMCALErr=0;
714  SignifEMCAL=0;
715  SignifEMCALErr=0;
716 
717  //DCAL
718  Npi0DCAL=0;
719  Npi0DCALErr=0;
720  MeanPosDCAL=0;
721  MeanPosDCALErr=0;
722  WidthDCAL=0;
723  WidthDCALErr=0;
724  Chi2NdfPi0DCAL=0;
725  NggDCAL=0;
726  NggDCALErr=0;
727  SignifDCAL=0;
728  SignifDCALErr=0;
729 
730  memset (CellMeanSM, 0, sizeof (Double_t) * nMax);
731  memset (CellRMSSM, 0, sizeof (Double_t) * nMax);
732  memset (ClusterMeanSM, 0, sizeof (Double_t) * nMax);
733  memset (ClusterTotSM, 0, sizeof (Double_t) * nMax);
734  memset (ClusterRMSSM, 0, sizeof (Double_t) * nMax);
735  memset (EtotalMeanSM, 0, sizeof (Double_t) * nMax);
736  memset (EtotalRMSSM, 0, sizeof (Double_t) * nMax);
737  memset (CellPerClusterMeanSM, 0, sizeof (Double_t) * nMax);
738  memset (CellPerClusterRMSSM, 0, sizeof (Double_t) * nMax);
739  memset (ECell1MeanSM, 0, sizeof (Double_t) * nMax);
740  memset (ECell1RMSSM, 0, sizeof (Double_t) * nMax);
741 
742  memset (MeanPosSM, 0, sizeof (Double_t) * nMax);
743  memset (MeanPosErrSM, 0, sizeof (Double_t) * nMax);
744  memset (WidthSM, 0, sizeof (Double_t) * nMax);
745  memset (WidthErrSM, 0, sizeof (Double_t) * nMax);
746  memset (Npi0SM, 0, sizeof (Double_t) * nMax);
747  memset (Npi0ErrSM, 0, sizeof (Double_t) * nMax);
748 
749  TString dirname;
750  if(!fTrigger.Contains("QA")) {
751  dirname = "CaloQA_";
752  }
753  dirname += fTrigger;
754 
755  Bool_t dirok = f->cd(dirname);
756  if(!dirok) { Error(__FUNCTION__,Form("No input directory %s",dirname.Data())); tree->Fill(); ftree->cd(); tree->Write(); ret= -1; continue;}
757  outputList = (TList*)gDirectory->Get(dirname);
758  if(!outputList){ Error(__FUNCTION__,Form("No input list! %s",dirname.Data())); tree->Fill(); ftree->cd(); tree->Write(); ret=-2; continue;;}
759  outputList->SetOwner();
760 
761  // number of events
762  fhNEvents =(TH1F *)outputList->FindObject("hNEvents");
763  if(!fhNEvents){ Error(__FUNCTION__,Form("NEvent histogram not found for trigger %s",fTrigger.Data())); tree->Fill(); ftree->cd(); tree->Write(); ret=-3; continue;}
764  Nevent=fhNEvents->GetEntries();
765  if(Nevent==0) {Error(__FUNCTION__,Form("No event in trigger %s",fTrigger.Data())); tree->Fill(); ftree->cd(); tree->Write(); ret=-4; continue;}
766  if(Nevent<20) {Error(__FUNCTION__,Form("Less than 20 events in trigger %s",fTrigger.Data())); tree->Fill(); ftree->cd(); tree->Write(); ret=-5; continue;}
767  TH1F* hE =(TH1F *)outputList->FindObject("EMCAL_hE");
768  if (!hE || (!hE->GetEntries())) { Error(__FUNCTION__,Form("hE Histogram not found or it is empty for trigger %s, EMCAL was not in this run?",fTrigger.Data())); EtotalMean=0; tree->Fill(); ftree->cd(); tree->Write(); ret=-4; continue;}
769 
770  // first do clusters trending
771  fhE = (TH1F *)outputList->FindObject(fCalorimeter+"_hE");
772  Double_t energy = 0. ;
773  fhECharged = (TH1F *)outputList->FindObject(fCalorimeter+"_hECharged");
774 
775  for(Int_t ibin = fhE->FindBin(0.6) ; ibin <fhE->FindBin(50.) ; ibin++){
776  energy+=fhE->GetBinCenter(ibin)*fhE->GetBinContent(ibin);
777  }
778  if(fhE->Integral(fhE->FindBin(0.6), fhE->FindBin(50.))==0){Error(__FUNCTION__,Form("Not enough events")); tree->Fill(); ftree->cd(); tree->Write(); ret=-6; continue;}
779  EtotalMean=energy/fhE->Integral(fhE->FindBin(0.6), fhE->FindBin(50.)) ;
780  EtotalRMS=fhE->GetMeanError();
781 
782  TString nameNCell = Form("%s_hNCells_Mod",fCalorimeter.Data());
783  TString nameNCluster = Form("%s_hNClusters_Mod",fCalorimeter.Data());
784  TString nameE = Form("%s_hE_Mod",fCalorimeter.Data());
785  TH2F* hNCellsMod= (TH2F*)outputList->FindObject(nameNCell.Data());
786  TH2F* hNClusterMod=(TH2F*)outputList->FindObject(nameNCluster.Data());
787  TH2F* hEMod=(TH2F*)outputList->FindObject(nameE.Data());
788 
789  if (!hNCellsMod || !hNClusterMod || !hEMod) {Error(__FUNCTION__,"A requiered histogram was not found (the imput QAresult.root might be too old)!"); tree->Fill(); ftree->cd(); tree->Write(); ret=-7; continue;}
790 
791  TCanvas* c1 = new TCanvas("Pi0InvMassSM","Pi0 Invariant Mass for each SM", 600, 600);
792  gStyle->SetOptStat(1); // MG modif
793  gStyle->SetOptFit(1); // MG modif
794 
795  c1->SetFillColor(0);
796  c1->SetBorderSize(0);
797  c1->SetFrameBorderMode(0);
798  // c1->SetOptStat(1); // MG modif
799  Bool_t mod3=0; if (n%3) mod3=1;
800  c1->Divide(3,n/3+mod3);
801 
802  //per sm trending
803  TString nameNCellPerCluster;
804  for(Int_t ism = 0 ; ism < n ; ism++){
805  cout << "#########################"<< endl;
806  cout << " Super Module " << ism << " Run " << RunId << endl;
807  // first do clusters trending
808  nameNCellPerCluster = Form("%s_hNCellsPerCluster_Mod%d",fCalorimeter.Data(),ism);
809  NCellsPerCluster[ism] = (TH2F*)outputList->FindObject(nameNCellPerCluster.Data());
810  if(! (TH2F*)outputList->FindObject(nameNCellPerCluster.Data()) ) { Error(__FUNCTION__,Form("NCellsPerCluster histogram not found for super module %i",ism));ret=-8; continue;}
811  NCellsPerCluster[ism] = (TH2F*)outputList->FindObject(nameNCellPerCluster.Data());
812 
813  NCells[ism] = (TH1F*)hNCellsMod->ProjectionX(Form("NCells%d",ism),ism+1,ism+1,"");
814  NClusters[ism] = (TH1F*)hNClusterMod->ProjectionX(Form("NClusters%d",ism),ism+1,ism+1,"");
815  E[ism] = (TH1F*)hEMod->ProjectionX(Form("E%d",ism),ism+1,ism+1,"");
816  CellMeanSM[ism]=NCells[ism]->GetMean();
817 
818  CellRMSSM[ism]=NCells[ism]->GetMeanError();
819  Int_t binmax;
820  binmax = NClusters[ism]->GetNbinsX();
821  for(Int_t ibin = 1 ;ibin <=binmax ;ibin++)
822  { ClusterTotSM[ism]+=(NClusters[ism]->GetBinCenter(ibin))*(NClusters[ism]->GetBinContent(ibin));
823  }
824  if (NClusters[ism]->GetEntries()>0)
825  ClusterTotSM[ism]= ClusterTotSM[ism]/NClusters[ism]->GetEntries();
826 
827  // ClusterTotSM[ism]=NClusters[ism]->GetMean();
828  ClusterMeanSM[ism]=NClusters[ism]->GetMean();
829  ClusterRMSSM[ism]=NClusters[ism]->GetMeanError();
830  CellPerClusterMeanSM[ism]=NCellsPerCluster[ism]->GetMean(2);
831  CellPerClusterRMSSM[ism]=NCellsPerCluster[ism]->GetMeanError(2);
832 
833  ECell1MeanSM[ism] =NCellsPerCluster[ism]->ProjectionX("",2,50,"")->Integral(5,50)/(Nevent);
834  ECell1RMSSM[ism] =NCellsPerCluster[ism]->ProjectionX("",2,50,"")->GetMeanError();
835  Double_t energySM = 0. ;
836  for(Int_t ibin = E[ism]->FindBin(0.6) ; ibin <E[ism]->FindBin(50.) ; ibin++){
837  energySM+=E[ism]->GetBinCenter(ibin)*(E[ism]->GetBinContent(ibin));
838  }
839  if(E[ism]->Integral(E[ism]->FindBin(0.6),E[ism]->FindBin(50.))==0){Error(__FUNCTION__,Form("Energy: Not enough events/SM")); continue;}
840  EtotalMeanSM[ism]=energySM/(E[ism]->Integral(E[ism]->FindBin(0.6),E[ism]->FindBin(50.)));
841 
842  EtotalRMSSM[ism]=E[ism]->GetMeanError();
843 
844  if(ism==0) {
845  fhNCells = (TH1F*)NCells[ism]->Clone("NCells");
846  fhNClusters = (TH1F*)NClusters[ism]->Clone("NClusters");
847  }
848  else {
849  fhNCells->Add(NCells[ism],1);
850  fhNClusters->Add(NClusters[ism],1);
851  }
852 
853  //Pi0
854  c1->cd(ism+1);
855  TString namePair = Form("%s_hIM_Mod%d",fCalorimeter.Data(),ism);
856  IM[ism] = (TH2F*)outputList->FindObject(namePair.Data());
857  IM[ism]->Sumw2();
858 
859  TString projname = Form("SM_%d",ism);
860  MggSM[ism] = (TH1F *)IM[ism]->ProjectionY(projname.Data(), 2, 20, "") ;//MG modif
861 
862 
863  if(MggSM[ism]->GetEntries()>100) {
864  fitMass->SetParameter(0, MggSM[ism]->GetBinContent(MggSM[ism]->GetMaximumBin()));
865  fitMass->SetParameter(1, mPDG); //
866  fitMass->SetParameter(2, 15.); //
867  fitMass->SetParameter(3,0.);
868  fitMass->SetParameter(4,MggSM[ism]->GetBinContent(MggSM[ism]->FindBin(0.11)));
869  fitMass->SetParameter(5,MggSM[ism]->GetBinContent(MggSM[ism]->FindBin(0.20)));
870 
871  if(MggSM[ism]->GetEntries()<1000){ MggSM[ism]->Rebin(5);} else {MggSM[ism]->Rebin();}
872  MggSM[ism]->Fit("fitMass", "WL R +","",0.05, 0.30);
873  MggSM[ism]->SetTitle(Form("Pi0 Mass for super module %i",ism));
874  MggSM[ism]->SetTitleSize(0.1);
875  MggSM[ism]->SetXTitle("Pi0 Mass");
876  MggSM[ism]->SetYTitle("Nb of entries");
877  MggSM[ism]->GetXaxis()->SetLabelSize(0.05);
878  MggSM[ism]->GetXaxis()->SetTitleSize(0.07);
879  MggSM[ism]->GetXaxis()->SetTitleOffset(0.68);
880  MggSM[ism]->GetYaxis()->SetLabelSize(0.05);
881  MggSM[ism]->GetYaxis()->SetTitleSize(0.06);
882  MggSM[ism]->GetYaxis()->SetTitleOffset(0.78);
883 
884  MeanPosSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(1)*1000;
885  MeanPosErrSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParError(1)*1000;
886  WidthSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(2)*1000;
887  WidthErrSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParError(2)*1000;
888  Npi0SM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(0)*(MggSM[ism]->GetFunction("fitMass")->GetParameter(2))*TMath::Sqrt(2*TMath::Pi())/(Nevent*MggSM[ism]->GetBinWidth(1));
889  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))
890  +(MggSM[ism]->GetFunction("fitMass")->GetParError(2)/MggSM[ism]->GetFunction("fitMass")->GetParameter(2))*(MggSM[ism]->GetFunction("fitMass")->GetParError(2)/MggSM[ism]->GetFunction("fitMass")->GetParameter(2)));
891  Npi0ErrSM[ism] = 0.; //
892 
893  }// end if enough events for Pi0 fit and trending
894  else { Info(__FUNCTION__,Form("Not enough events for Pi0 fit and trending for super module %i",ism));} ;
895  } //per SM loop
896 
897  // Now Pi0 global trending in EMCAL
898  TCanvas* c2 = new TCanvas("Pi0InvMassEMCAL","Pi0 Invariant Mass in EMCAL", 600, 600);
899  c2->SetFillColor(0);
900  c2->SetBorderSize(0);
901  c2->SetFrameBorderMode(0);
902  // c2->SetOptStat(1); // MG modif
903 
904  fhIM = (TH2F *)outputList->FindObject(fCalorimeter+"_hIM");
905  fhIM->Sumw2();
906  fhMggEMCAL = (TH1F *)fhIM->ProjectionY("MggEMCAL", 2, 20, "") ; // to modify projection range
907  if(fhMggEMCAL->GetEntries()==0) {Error(__FUNCTION__,"The Pi0 histogram in EMCAL is empty !"); tree->Fill(); ret=-8; continue;}
908  fitMass->SetParameter(0, 4500);
909  fitMass->SetParameter(1, mPDG);
910  fitMass->SetParameter(2, 0.01);
911  fitMass->SetParameter(3,0.);
912  fitMass->SetParameter(4,fhMggEMCAL->GetBinContent(fhMggEMCAL->FindBin(0.11)));
913  fitMass->SetParameter(5,fhMggEMCAL->GetBinContent(fhMggEMCAL->FindBin(0.20)));
914 
915  if(fhMggEMCAL->GetEntries()<5000){
916  fhMggEMCAL->Rebin(5);} //MG modif
917  else fhMggEMCAL->Rebin();
918 
919  fhMggEMCAL->Fit("fitMass", "L R +", "", 0.05, 0.20);
920 
921  fhMggEMCAL->SetTitle("Pi0 Mass in EMCAL");
922  fhMggEMCAL->SetTitleSize(0.1);
923  fhMggEMCAL->SetXTitle("Pi0 Mass");
924  fhMggEMCAL->SetYTitle("Nb of entries");
925  fhMggEMCAL->GetXaxis()->SetLabelSize(0.03);
926  fhMggEMCAL->GetXaxis()->SetTitleSize(0.03);
927  fhMggEMCAL->GetXaxis()->SetTitleOffset(1.3);
928  fhMggEMCAL->GetYaxis()->SetLabelSize(0.03);
929  fhMggEMCAL->GetYaxis()->SetTitleSize(0.03);
930  fhMggEMCAL->GetYaxis()->SetTitleOffset(1.3);
931 
932  MeanPosEMCAL = fhMggEMCAL->GetFunction("fitMass")->GetParameter(1)*1000;
933  MeanPosEMCALErr = fhMggEMCAL->GetFunction("fitMass")->GetParError(1)*1000;
934  WidthEMCAL = fhMggEMCAL->GetFunction("fitMass")->GetParameter(2)*1000;
935  WidthEMCALErr = fhMggEMCAL->GetFunction("fitMass")->GetParError(2)*1000;
936  Chi2NdfPi0EMCAL = fhMggEMCAL->GetFunction("fitMass")->GetChisquare()/fhMggEMCAL->GetFunction("fitMass")->GetNDF();
937  Npi0EMCAL = fhMggEMCAL->GetFunction("fitMass")->GetParameter(0)*fhMggEMCAL->GetFunction("fitMass")->GetParameter(2)*TMath::Sqrt(2*TMath::Pi())/(Nevent*fhMggEMCAL->GetBinWidth(10));
938  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));
939  Npi0EMCALErr = 0.; //
940  NggEMCAL = fhMggEMCAL->GetFunction("fitMass")->Integral(0.11, 0.16)/(Nevent*fhMggEMCAL->GetBinWidth(10));
941  NggEMCALErr = fhMggEMCAL->GetFunction("fitMass")->IntegralError(0.11, 0.16)/(fhMggEMCAL->Integral()*fhMggEMCAL->GetBinWidth(10));
942  SignifEMCAL = Npi0EMCAL/NggEMCAL;
943  SignifEMCALErr = TMath::Sqrt((Npi0EMCALErr/Npi0EMCAL*(Npi0EMCALErr/Npi0EMCAL)+(NggEMCALErr/NggEMCAL*(NggEMCALErr/NggEMCAL))));
944  SignifEMCALErr = SignifEMCAL*SignifEMCALErr;
945 
946  TCanvas* c3 = new TCanvas("Pi0InvMassDCAL","Pi0 Invariant Mass in DCAL", 600, 600);
947  c3->SetFillColor(0);
948  c3->SetBorderSize(0);
949  c3->SetFrameBorderMode(0);
950 
951  fhIMDCAL = (TH2F *)outputList->FindObject(fCalorimeter+"_hIMDCAL");
952  if(fhIMDCAL)
953  {
954 
955  // Now Pi0 global trending in DCAL
956  // TCanvas* c3 = new TCanvas("Pi0InvMassDCAL","Pi0 Invariant Mass in DCAL", 600, 600);
957  // c3->SetFillColor(0);
958  // c3->SetBorderSize(0);
959  // c3->SetFrameBorderMode(0);
960  // c2->SetOptStat(1); // MG modif
961 
962  // fhIMDCAL = (TH2F *)outputList->FindObject(fCalorimeter+"_hIMDCAL");
963  // if(!fhIMDCAL) continue;
964  fhIMDCAL->Sumw2();
965  fhMggDCAL = (TH1F *)fhIMDCAL->ProjectionY("MggDCAL", 2, 20, "") ; // to modify projection range
966  if(fhMggDCAL->GetEntries()==0) {Error(__FUNCTION__,"The Pi0 histogram in DCAL is empty !"); tree->Fill(); ret=-8; continue;}
967  fitMass->SetParameter(0, 4500);
968  fitMass->SetParameter(1, mPDG);
969  fitMass->SetParameter(2, 0.01);
970  fitMass->SetParameter(3,0.);
971  fitMass->SetParameter(4,fhMggDCAL->GetBinContent(fhMggDCAL->FindBin(0.11)));
972  fitMass->SetParameter(5,fhMggDCAL->GetBinContent(fhMggDCAL->FindBin(0.20)));
973 
974  if(fhMggDCAL->GetEntries()<5000){
975  fhMggDCAL->Rebin(5);} //MG modif
976  else fhMggDCAL->Rebin();
977 
978  fhMggDCAL->Fit("fitMass", "L R +", "", 0.05, 0.20);
979 
980  fhMggDCAL->SetTitle("Pi0 Mass in DCAL");
981  fhMggDCAL->SetTitleSize(0.1);
982  fhMggDCAL->SetXTitle("Pi0 Mass");
983  fhMggDCAL->SetYTitle("Nb of entries");
984  fhMggDCAL->GetXaxis()->SetLabelSize(0.03);
985  fhMggDCAL->GetXaxis()->SetTitleSize(0.03);
986  fhMggDCAL->GetXaxis()->SetTitleOffset(1.3);
987  fhMggDCAL->GetYaxis()->SetLabelSize(0.03);
988  fhMggDCAL->GetYaxis()->SetTitleSize(0.03);
989  fhMggDCAL->GetYaxis()->SetTitleOffset(1.3);
990 
991  MeanPosDCAL = fhMggDCAL->GetFunction("fitMass")->GetParameter(1)*1000;
992  MeanPosDCALErr = fhMggDCAL->GetFunction("fitMass")->GetParError(1)*1000;
993  WidthDCAL = fhMggDCAL->GetFunction("fitMass")->GetParameter(2)*1000;
994  WidthDCALErr = fhMggDCAL->GetFunction("fitMass")->GetParError(2)*1000;
995  Chi2NdfPi0DCAL = fhMggDCAL->GetFunction("fitMass")->GetChisquare()/fhMggDCAL->GetFunction("fitMass")->GetNDF();
996  Npi0DCAL = fhMggDCAL->GetFunction("fitMass")->GetParameter(0)*fhMggDCAL->GetFunction("fitMass")->GetParameter(2)*TMath::Sqrt(2*TMath::Pi())/(Nevent*fhMggDCAL->GetBinWidth(10));
997  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));
998  Npi0DCALErr = 0.; //
999  NggDCAL = fhMggDCAL->GetFunction("fitMass")->Integral(0.11, 0.16)/(Nevent*fhMggDCAL->GetBinWidth(10));
1000  NggDCALErr = fhMggDCAL->GetFunction("fitMass")->IntegralError(0.11, 0.16)/(fhMggDCAL->Integral()*fhMggDCAL->GetBinWidth(10));
1001  SignifDCAL = Npi0DCAL/NggDCAL;
1002  SignifDCALErr = TMath::Sqrt((Npi0DCALErr/Npi0DCAL*(Npi0DCALErr/Npi0DCAL)+(NggDCALErr/NggDCAL*(NggDCALErr/NggDCAL))));
1003  SignifDCALErr = SignifDCAL*SignifDCALErr;
1004 
1005  }
1006 
1007  cout<<"******************"<<endl;
1008  //end of global trending
1009  NTClusters=fhE->IntegralAndError(fhE->GetXaxis()->FindBin(0.),fhE->GetXaxis()->FindBin(200.),NTClustersRMS);
1010  cout<<"********************************** nb total de clusters ***************************************************"<<endl;
1011  NCClusters=fhECharged->IntegralAndError(fhECharged->GetXaxis()->FindBin(0.),fhECharged->GetXaxis()->FindBin(50.),NCClustersRMS);
1012  NMatchClustersP = NCClusters/NTClusters;
1013  cout<<" ici se trouve la valeur recherchée: "<<NMatchClustersP<<" "<<endl;
1014  NMatchClustersPRMS = NMatchClustersP* TMath::Sqrt(TMath::Power(NCClustersRMS/NCClusters,2)+TMath::Power(NTClustersRMS/NTClusters,2));
1015  ClusterMean=fhNClusters->GetMean();
1016  ClusterRMS=fhNClusters->GetMeanError();
1017  CellMean=fhNCells->GetMean();
1018  CellRMS=fhNCells->GetMeanError();
1019  tree->Fill();
1020 
1021  TString outfilename = QAPATH + "Pi0InvMassEMCAL" + fTrigger(r) + ".pdf";
1022  TString outfilename2 = QAPATH + "Pi0InvMassEMCAL" + fTrigger(r) + ".png";
1023  if(SavePlots==2) c2->SaveAs(outfilename);
1024  if(SavePlots) c2->SaveAs(outfilename2);
1025 
1026  outfilename = QAPATH + "Pi0InvMassSM" + fTrigger(r) + ".pdf";
1027  outfilename2 = QAPATH + "Pi0InvMassSM" + fTrigger(r) + ".png";
1028  if(SavePlots==2) c1->SaveAs(outfilename);
1029  if(SavePlots) c1->SaveAs(outfilename2);
1030 
1031  cout<<"************************************************** "<<n<<" **************************************************"<<endl;
1032 
1033  if(fhIMDCAL)
1034  {
1035  TString outfilename4 = QAPATH + "Pi0InvMassDCAL" + fTrigger(r) + ".pdf";
1036  TString outfilename5= QAPATH + "Pi0InvMassDCAL" + fTrigger(r) + ".png";
1037  if(SavePlots==2) c3->SaveAs(outfilename4);
1038  if(SavePlots) c3->SaveAs(outfilename5);
1039  }
1040 
1041 
1042  fout->cd();
1043  fout->Cd(Form("%s/%s/%ld/%s/%s",period.Data(),pass.Data(),RunId,"RunLevelQA",fTrigger.Data()));
1044  c2->Write();
1045  c1->Write();
1046  if(fhIMDCAL) c3->Write();
1047  delete c1;
1048  delete c2;
1049  if(fhIMDCAL) delete c3;
1050  if (outputList) outputList->Delete() ;
1051 
1052  QAData << RunId<<" "<< Nevent
1053  <<"\n";
1054 
1055  QAData.close();
1056 
1057  }
1058 
1059  ftree->cd();
1060  tree->Write();
1061  ftree->Close();
1062 
1063  return ret;
1064 
1065 }
1066 
1067 //-------------------------------------------------------------------------
1068 TH2F* FormatRunHisto(TH2F* aHisto,const char* title,const char* YTitle)
1069 {
1070 
1071  if(!aHisto) {Error(__FUNCTION__,Form("The histogram with title \"%s\" was not found!",title)); return new TH2F();}
1072  aHisto->SetStats(kFALSE);
1073  aHisto->SetTitle(title);
1074  aHisto->SetStats(kFALSE);
1075  aHisto->SetYTitle(YTitle);
1076  aHisto->GetYaxis()->SetTitleOffset(1.2);
1077  aHisto->GetYaxis()->SetLabelSize(0.03);
1078  aHisto->GetZaxis()->SetLabelSize(0.02);
1079 
1080  return aHisto;
1081 
1082 }
1083 
1084 //--------------------------------------------------------------------------------------------------------------
1085 TH2F* HistoPerMod(TH2F* hTmpPerMod,const char* title)
1086 {
1087 
1088  if(!hTmpPerMod) {Error(__FUNCTION__,Form("The histogram with title \"%s\" was not found!",title)); return new TH2F();}
1089  hTmpPerMod->SetStats(kFALSE);
1090  hTmpPerMod->SetTitle(title);
1091  hTmpPerMod->SetTitleSize(0.1);
1092  hTmpPerMod->GetXaxis()->SetTitleOffset(1.1);
1093  hTmpPerMod->GetXaxis()->SetTitleSize(0.05);
1094  hTmpPerMod->GetXaxis()->SetLabelSize(0.06);
1095  hTmpPerMod->GetYaxis()->SetTitleOffset(1.1);
1096  hTmpPerMod->GetYaxis()->SetTitleSize(0.05);
1097  hTmpPerMod->GetYaxis()->SetLabelSize(0.06);
1098  hTmpPerMod->GetZaxis()->SetLabelSize(0.04);
1099 
1100  return hTmpPerMod;
1101 
1102 }
1103 
1104 //---------------------------------------------------------------------------------------------------
1105 TH2* AutoZoom(TH2* H,Option_t* aType, Int_t EntryMin)
1106 {
1107 
1108  Int_t shiftX = (Int_t)(H->GetNbinsX()/30.);
1109  Int_t shiftY = (Int_t)(H->GetNbinsY()/30.);
1110 
1111  TString opt = aType;
1112  opt.ToLower();
1113 
1114  int minX = 0;
1115  int maxX = H->GetNbinsX();
1116  int New_minX = minX;
1117  int New_maxX = maxX;
1118 
1119  int minY = 0;
1120  int maxY = H->GetNbinsY();
1121  int New_minY = minY;
1122  int New_maxY = maxY;
1123 
1124  if (opt.Contains("all")) opt = TString("minx,maxx,miny,maxy");
1125 
1126  if (opt.Contains("maxx"))
1127  {
1128 
1129  for (New_maxX = maxX;New_maxX >=minX; New_maxX--)
1130  { Stat_t c = 0;
1131  for (int i_y = maxY; i_y >= minY;i_y--)
1132  { c = H->GetBinContent(New_maxX,i_y); if (c>EntryMin) break;}
1133  if (c>EntryMin) break;
1134  }
1135  }
1136 
1137  if (opt.Contains("maxy"))
1138  {
1139 
1140  for (New_maxY = maxY;New_maxY >=minY;New_maxY--)
1141  { Stat_t c = 0;
1142  for (int i_x=maxX; i_x>=minX;i_x--)
1143  { c = H->GetBinContent(i_x, New_maxY ); if (c>EntryMin) break;}
1144  if (c>EntryMin) break;
1145  }
1146 
1147  }
1148 
1149  if (opt.Contains("minx"))
1150  {
1151 
1152  for (New_minX = minX;New_minX <=maxX; New_minX++)
1153  { Stat_t c = 0;
1154  for (int i_y = minY; i_y <= maxY;i_y++)
1155  { c = H->GetBinContent(New_minX,i_y); if (c>EntryMin) break;}
1156  if (c>EntryMin) break;
1157  }
1158  }
1159 
1160  if (opt.Contains("miny"))
1161  {
1162  for (New_minY = minY;New_minY <=maxY;New_minY++)
1163  { Stat_t c = 0;
1164  for (int i_x=minX; i_x<=maxX;i_x++)
1165  { c = H->GetBinContent(i_x, New_minY ); if (c>EntryMin) break;}
1166  if (c>EntryMin) break;
1167  }
1168  }
1169 
1170  if (New_maxX!=-1 && New_maxY!=-1) H->GetXaxis()->SetRange(New_minX - shiftX , New_maxX + shiftX);
1171  if (New_maxX!=-1 && New_maxY!=-1) H->GetYaxis()->SetRange(New_minY - shiftY , New_maxY + shiftY);
1172 
1173  return H;
1174 
1175 }
1176 
1177 //----------------------------------------------------------------------------------------------------
1178 int FindNumberOfSM(TFile* f, TString fTrigger, TString period)
1179 {
1180 
1181  TString direct;
1182  if(!fTrigger.Contains("QA")) {
1183  direct = "CaloQA_";
1184  }
1185  direct += fTrigger;
1186 
1187  Int_t nSMt=-1;
1188  Int_t year = 2000 + TString(period(3,2)).Atoi();
1189  if ( year == 2010 ) { nSMt=4; }
1190  else if ( year == 2011 || year == 2012 ) { nSMt=10; }
1191  else if ( year == 2013 || year == 2014 ) { nSMt=12; }
1192  else { nSMt=20; }
1193 
1194  TList *outputList;
1195  Bool_t dirok = f->cd(direct);
1196  if (!dirok) { Error(__FUNCTION__,Form("No input directory %s, the number SMs will be returned based on the year!",direct.Data()));}
1197  else { outputList = (TList*)gDirectory->Get(direct);}
1198  if(!outputList) { Error(__FUNCTION__,"No input list, the number SMs will be returned based on the year! ");}
1199  else {
1200  outputList->SetOwner();
1201  TH2F* hNSM =(TH2F *)outputList->FindObject("EMCAL_hE_Mod");
1202  if (!hNSM || (!hNSM->GetEntries())) { Error(__FUNCTION__,"hNSM Histogram not found or it is empty, the number SMs will be returned based on the year!");}
1203  else {
1204  nSMt = hNSM->GetYaxis()->GetBinUpEdge(hNSM->FindLastBinAbove(0,2));
1205  }
1206  }
1207  if (outputList) {outputList->Delete();}
1208 
1209  return nSMt;
1210 
1211 }
TH2F * FormatRunHisto(TH2F *aHisto, const char *title, const char *YTitle="")
TString QAPATH
const char * title
Definition: MakeQAPdf.C:26
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)
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_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)
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)
energy
Double_t fitE(Double_t *x, Double_t *par)
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")