PWGPP/EMCAL/QA/macros/CreateEMCALRunQA.C

   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
  23 // This macro produces runLevelQA for EMCAL from a QAresults.root file
  24 // Authors: Y. Mao, A. Mas, M. Germain & A.Shabetai  SUBATECH
  25 // re-factored for automatic QA processing A.SHABETAI
  26
  27 Int_t DrawOccupancy(Long_t run, TString period, TString pass, TString fTrigger, TString system, TFile* f, TFile* fout, AliEMCALGeometry* geom, Int_t SavePlots);
  28 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);
  29 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);
  30
  31 TH2F* FormatRunHisto(TH2F* aHisto,const char* title,const char* YTitle="");
  32 TH2F* HistoPerMod(TH2F* name,const char* title);
  33 TH2F* AutoZoom(TH2F* H,Option_t* aType="all", Int_t EntryMin=0);
  34 int FindNumberOfSM(TFile* f, TString fTrigger,TString period);
  35
  36 TString QAPATH;
  37 TString QAPATHF = "./";
  38 //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  39 void set_plot_style()
  40 {
  41   const Int_t NRGBs = 5;
  42   const Int_t NCont = 255;
  43
  44   Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
  45   Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
  46   Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
  47   Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
  48   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
  49   gStyle->SetNumberContours(NCont);
  50
  51 }
  52
  53 //-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  54 Double_t pi0massP2(Double_t *x, Double_t *par)
  55 {
  56   Double_t gaus;
  57
  58   if (par[2] != 0.) gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) / (2*par[2]*par[2]) );
  59
  60   else gaus = 99999999.;
  61
  62   Double_t back = par[3] + par[4]*x[0] + par[5]*x[0]*x[0];
  63
  64   return gaus+back;
  65
  66 }
  67
  68 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  69 Double_t pi0massP1(Double_t *x, Double_t *par)
  70 {
  71   Double_t gaus = par[0] * TMath::Exp( -(x[0]-par[1])*(x[0]-par[1]) / (2*par[2]*par[2]) );
  72
  73   Double_t back = par[3] + par[4]*x[0];
  74
  75   return gaus+back;
  76
  77 }
  78
  79 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  80 Double_t fitE(Double_t *x, Double_t *par)
  81 {
  82
  83   Double_t levy;
  84
  85   levy = par[0] * TMath::Exp( -par[1]/x[0]) * TMath::Power(x[0], -par[2]) ;
  86
  87   return levy;
  88 }
  89
  90 //--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
  91 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")
  92 {
  93
  94   QAPATH = TString(gSystem->Getenv("QAPATH"));
  95   if(QAPATH.IsNull()) QAPATH = QAPATHF;
  96   if(! QAPATH.BeginsWith("./")) { QAPATH = QAPATH + RunId + "/";}
  97
  98   AliLog::SetGlobalLogLevel(AliLog::kError);
  99   TFile *f = new TFile(filename);
 100   AliLog::SetGlobalLogLevel(AliLog::kInfo);
 101
 102   if (f->IsZombie()) {Error(__FUNCTION__,Form("Error openning the input file %s",filename)); return -1;}
 103
 104   TList* TriggersList = new TList();
 105
 106   if (fTrigger=="")
 107     {
 108       TPMERegexp name_re("CaloQA_\\w+");
 109       TObjLink* link = f->GetListOfKeys()->FirstLink();
 110
 111       while (link)
 112         {
 113           TString name = link->GetObject()->GetName();
 114           if (name_re.Match(name))
 115             {
 116               TriggersList->Add(link->GetObject());
 117               if(TString(filename).Contains("barrel") && ! name.Contains("default"))  TriggersList->Remove(link->GetObject());
 118               if(TString(filename).Contains("outer")  && ! name.Contains("EMC"))      TriggersList->Remove(link->GetObject());
 119             }
 120           link = link->Next();
 121         }
 122     } else {TriggersList->Add(new TObjString(fTrigger.Data()));}
 123
 124   if(!TriggersList->GetEntries()) {Error(__FUNCTION__,"No trigger found!"); return -2;}
 125
 126   int nSM = FindNumberOfSM(f,((TObjString*)TriggersList->Last())->GetString(),period);
 127   if (nSM<0) {Error(__FUNCTION__,"Could not find the number of super modules!"); return -3;}
 128   Info(__FUNCTION__,Form("%i super modules were discuvered",nSM));
 129   TString GeomName;
 130   if (nSM <= 6)         { nSM=6;  GeomName = "EMCAL_FIRSTYEARv1";}
 131   else if (nSM <= 10)   { nSM=10; GeomName = "EMCAL_COMPLETEv1";}
 132   else if (nSM <= 12)   { nSM=12; GeomName = "EMCAL_COMPLETE12SMv1";}
 133   else    {nSM = 20;              GeomName = "EMCAL_COMPLETE12SMv1_DCAL_8SM";}
 134
 135   AliEMCALGeometry *geom = new AliEMCALGeometry(GeomName.Data(),"EMCAL");
 136   Info(__FUNCTION__,Form("Using %i super modules and the Geometry %s",nSM,GeomName.Data()));
 137
 138   TFile *fout = new TFile(TString( QAPATH + period+"_"+pass + fTrigger+"_"+ (Long_t)RunId.Atoi() +"_QAplots.root").Data(),"RECREATE");
 139
 140   if((system.IsNull()) && (period.EndsWith("h"))) {system = "PbPb";}
 141
 142   Int_t ret=0;
 143   TIter next(TriggersList);
 144   while (TObject *obj = next())
 145     {
 146       fTrigger= TString(obj->GetName());
 147       ret -= DrawOccupancy(RunId.Atoi(),period,pass,fTrigger,system,f,fout,geom,SavePlots);
 148       ret -= DrawRun(RunId.Atoi(),period,pass,fTrigger,f,fout,SavePlots,nSM,filter);
 149     }
 150   ret-= TrendingEMCALTree(RunId.Atoi(),fCalorimeter,system,period,pass,nSM,TriggersList,f,fout,SavePlots);
 151
 152   f->Close();
 153   fout->Close();
 154   delete f;
 155   delete geom;
 156
 157   return ret;
 158 }
 159
 160 //----------------------------------------------------------------------------------------------------------------------------------------------------------
 161 Int_t DrawOccupancy(Long_t  run , TString period, TString pass, TString fTrigger,TString system, TFile* f,TFile* fout, AliEMCALGeometry* geom, Int_t SavePlots)
 162 {
 163
 164   set_plot_style();
 165   gStyle->SetOptStat(0);
 166   TH1::AddDirectory(kFALSE);
 167   TH2D *hEnergyMapReal = new TH2D("hEnergyMapReal","",96,-48,48,120,-0,120);
 168   TH2D *hOccupancyMapReal = new TH2D("hOccupancyMapReal","",96,-48,48,120,-0,120);
 169   hEnergyMapReal->SetXTitle("eta (bin)");
 170   hEnergyMapReal->SetYTitle("phi (bin)");
 171   hEnergyMapReal->SetZTitle("E(GeV)/event");
 172   hEnergyMapReal->GetYaxis()->SetTitleOffset(1.2);
 173   hEnergyMapReal->GetZaxis()->SetLabelSize(0.02);
 174   hEnergyMapReal->GetZaxis()->SetTitleOffset(1.36);
 175
 176   hOccupancyMapReal->SetXTitle("eta (bin)");
 177   hOccupancyMapReal->SetYTitle("phi (bin)");
 178   hOccupancyMapReal->GetYaxis()->SetTitleOffset(1.2);
 179   hOccupancyMapReal->GetZaxis()->SetLabelSize(0.02);
 180
 181   Int_t nSupMod, nModule, nIphi, nIeta;
 182   Int_t iphi, ieta;
 183   Int_t realbineta=0;
 184   Int_t realbinphi=0;
 185
 186   //NO MASK
 187   Int_t mask[1] = {2222222};
 188
 189   TH2F *hCellAmplitude;
 190   TH1F *hNEvents;
 191   Int_t Events;
 192   Int_t n=0;
 193
 194   TString direct;
 195   if(!fTrigger.Contains("QA")) {
 196     direct = "CaloQA_";
 197   }
 198   direct += fTrigger;
 199   Bool_t dirok = f->cd(direct);
 200   if (!dirok) { Error(__FUNCTION__,Form("No input drectory %s",direct.Data())); return -1;}
 201   TList *outputList = (TList*)gDirectory->Get(direct);
 202   if(!outputList){ Error(__FUNCTION__,"No input list! "); return -1;}
 203   outputList->SetOwner();
 204
 205   fout->mkdir(Form("%s/%s/%ld/%s/%s",period.Data(),pass.Data(),run,"RunLevelQA",fTrigger.Data()));
 206   fout->cd();
 207   fout->Cd(Form("%s/%s/%ld/%s/%s",period.Data(),pass.Data(),run,"RunLevelQA",fTrigger.Data()));
 208
 209   hNEvents =(TH1F *)outputList->FindObject("hNEvents");
 210   if(!hNEvents){ Error(__FUNCTION__,Form("hNEvent histogram not found for trigger %s ",fTrigger.Data())); return -2;}
 211   Events = (Int_t)hNEvents->GetEntries();
 212   if(Events==0){ Error(__FUNCTION__,Form("No event in trigger %s",fTrigger.Data())); return -3;}
 213
 214   Double_t Eth=1;
 215   if(system=="PbPb"){
 216     Eth = 5.;
 217     if (fTrigger.Contains("EMC")) Eth=20.;
 218   }
 219   if(system=="pp"){
 220     Eth = 1.;
 221     if (fTrigger.Contains("EMC")) Eth=5.;
 222   }
 223
 224   hCellAmplitude =(TH2F *)outputList->FindObject("EMCAL_hAmpId");
 225
 226   for(Int_t i = 0; i < geom->GetNCells() ; i++){
 227     Double_t Esum = 0;
 228     Double_t Nsum = 0;
 229
 230     for (Int_t j = 1; j <= hCellAmplitude->GetNbinsX(); j++)
 231       {
 232         Double_t E = hCellAmplitude->GetXaxis()->GetBinCenter(j);
 233         Double_t N = hCellAmplitude->GetBinContent(j, i+1);
 234
 235         if (E < 0.3) continue;
 236
 237         if (E <= Eth) {
 238           Esum += E*N;
 239           Nsum += N;
 240         }
 241       }
 242
 243     Int_t absId = i;
 244     if(n!=0) {if(mask[n]<=mask[n-1]) Warning(__FUNCTION__,"The list of bad cells is not sorted !!");}
 245     if(i==mask[n]){n++ ; continue; } // skip bad cells
 246
 247     geom->GetCellIndex(absId,  nSupMod, nModule, nIphi, nIeta);
 248     geom->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta);
 249
 250     realbinphi = 120-(nSupMod/2)*24 -iphi -1; //
 251     if (nSupMod%2==0) realbineta= 48-ieta -1;
 252     if (nSupMod%2==1) realbineta= -ieta -1;
 253
 254     hEnergyMapReal->Fill(realbineta,realbinphi,Esum/(Double_t)Events);
 255     hOccupancyMapReal->Fill(realbineta,realbinphi,Nsum/(Double_t)Events);
 256   }
 257
 258   cout <<" Run: " << run << " trigger: " << fTrigger << " N_events: "<<Events<<endl;
 259
 260   TPMERegexp r("_\\w+");
 261   TString Energy;   Energy   = QAPATH + "MapEnergy"  + fTrigger(r) + ".pdf";
 262   TString Energy2;  Energy2  = QAPATH + "MapEnergy"  + fTrigger(r) + ".png";
 263   TString Entries;  Entries  = QAPATH + "MapEntries" + fTrigger(r) + ".pdf";
 264   TString Entries2; Entries2 = QAPATH + "MapEntries" + fTrigger(r) + ".png";
 265
 266   TCanvas *c1 = new TCanvas("Energymap","Energy Map",600,600);
 267   c1->SetFillColor(0);
 268   c1->SetGrid();
 269   c1->SetRightMargin(0.14);
 270   TString title = "run ";
 271   title += run ;
 272   if(fTrigger.Contains("EMC")) { title += " EMC ";} else {title += " MB ";}
 273   title += " Summed energy map";
 274   hEnergyMapReal->SetTitle(title);
 275   hEnergyMapReal->DrawCopy("colz");
 276   if(SavePlots==2) c1->SaveAs(Energy);
 277   if(SavePlots) c1->SaveAs(Energy2);
 278   c1->Write();
 279   delete c1;
 280
 281
 282   TCanvas *c2 = new TCanvas("Occupancy","Occupancy Map",600,600);
 283   c2->SetFillColor(0);
 284   c2->SetGrid();
 285   c2->SetRightMargin(0.14);
 286   TString title2 = "run ";
 287   title2 += run ;
 288   if(fTrigger.Contains("EMC")) { title2 += " EMC ";} else { title2 += " MB ";}
 289   title2 += " Occupancy map";
 290   hOccupancyMapReal->SetTitle(title2);
 291   hOccupancyMapReal->DrawCopy("colz");
 292   if(SavePlots==2) c2->SaveAs(Entries);
 293   if(SavePlots) c2->SaveAs(Entries2);
 294   c2->Write();
 295   delete c2;
 296
 297   if (outputList) {outputList->Delete();}
 298
 299   delete hEnergyMapReal;
 300   delete hOccupancyMapReal;
 301
 302   return 0;
 303 }
 304
 305 //-----------------------------------------------------------------------------------------------------------------------
 306 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)
 307 {
 308
 309   TString direct;
 310   if(!fTrigger.Contains("QA")) {
 311     direct = "CaloQA_";
 312   }
 313   direct += fTrigger;
 314
 315   f->cd(direct);
 316   if(!direct) { Error(__FUNCTION__,Form("No input directory %s",direct.Data())); return -1;}
 317   TList *outputList = (TList*)gDirectory->Get(direct);
 318   if(!outputList){ Error(__FUNCTION__,Form("No input list! %s",direct.Data())); return -2;}
 319   outputList->SetOwner();
 320   if (kFilter)
 321   {
 322      fout->mkdir(Form("%s/%s/%ld/%s",period.Data(),pass.Data(),run,fTrigger.Data()));
 323      fout->cd();
 324      fout->Cd(Form("%s/%s/%ld/%s",period.Data(),pass.Data(),run,fTrigger.Data()));
 325      outputList->Write();
 326   }
 327   fout->cd();
 328   fout->Cd(Form("%s/%s/%ld/%s/%s",period.Data(),pass.Data(),run,"RunLevelQA",fTrigger.Data()));
 329
 330
 331   set_plot_style();
 332   gStyle->SetPalette(1);
 333   TH1::AddDirectory(kFALSE);
 334   TString outfilename;
 335   TString outfilename2;
 336   const char* legend="";
 337   TPMERegexp r("_\\w+");
 338
 339   if (fTrigger.Contains("EMC")){ legend = Form(" Run %d EMC ",(int)run);}
 340   else legend = Form(" Run %d MB ",(int)run);
 341
 342   TH1F* hNEvents =(TH1F *)outputList->FindObject("hNEvents");
 343   if(!hNEvents){ Error(__FUNCTION__,Form("hNEvent histogram not found for trigger %s ",fTrigger.Data())); return -3;}
 344   Int_t Events = (Int_t)hNEvents->GetEntries();
 345   if(Events==0){ Error(__FUNCTION__,Form("No event in trigger %s",fTrigger.Data())); return -4 ;}
 346
 347   TCanvas* c1 = new TCanvas("TimeVsE", "Cluster Time Vs Energy", 600, 600);
 348   c1->SetLogz();
 349   c1->SetFillColor(0);
 350   c1->SetBorderSize(0);
 351   c1->SetFrameBorderMode(0);
 352
 353   TH2F*  hClusterTimeEnergy =(TH2F *)outputList->FindObject("EMCAL_hClusterTimeEnergy");
 354   if(!hClusterTimeEnergy) { Error(__FUNCTION__,Form("EMCAL_hClusterTimeEnergy: Histogram for trigger %s not found!",fTrigger.Data())); return -5;}
 355   FormatRunHisto(hClusterTimeEnergy,Form("Time Vs Energy%s",legend),"EMCAL ToF(ns)");
 356
 357   AutoZoom(hClusterTimeEnergy,"maxx")->DrawCopy("colz");
 358   outfilename =  QAPATH + "TimeRun" + fTrigger(r) + ".pdf" ;
 359   outfilename2 = QAPATH + "TimeRun" + fTrigger(r) + ".png" ;
 360
 361   if(SavePlots==2) c1->SaveAs(outfilename);
 362   if(SavePlots) c1->SaveAs(outfilename2);
 363   c1->Write();
 364   delete c1;
 365
 366   TCanvas  * c2 = new TCanvas("ClusterVsTrack ","Correlation calo Mult Vs Track Multiplicity", 600, 600);
 367   c2->SetLogz();
 368   c2->SetFillColor(0);
 369   c2->SetBorderSize(0);
 370   c2->SetFrameBorderMode(0);
 371
 372   TH2F* hClusterVsTrack =(TH2F *)outputList->FindObject("EMCAL_hCaloTrackMNClusters");
 373   FormatRunHisto(hClusterVsTrack,Form("N cluster Vs N track%s",legend));
 374
 375   AutoZoom(hClusterVsTrack,"maxx,maxy",1)->DrawCopy("colz");
 376   outfilename = QAPATH + "CaloTrackMult" + fTrigger(r) + ".pdf";
 377   outfilename2 = QAPATH + "CaloTrackMult" + fTrigger(r) + ".png";
 378   if(SavePlots==2) c2->SaveAs(outfilename);
 379   if(SavePlots) c2->SaveAs(outfilename2);
 380   c2->Write();
 381   delete c2;
 382
 383   TCanvas* c3 = new TCanvas("ClusterEVsTrack ","Correlation E calo Vs Track Multiplicity", 600, 600);
 384   c3->SetLogz();
 385   c3->SetFillColor(0);
 386   c3->SetBorderSize(0);
 387   c3->SetFrameBorderMode(0);
 388
 389   TH2F* hClusterEVsTrack =(TH2F*)outputList->FindObject("EMCAL_hCaloTrackMEClusters");
 390   FormatRunHisto(hClusterEVsTrack,Form("Sum E cluster Vs N track%s",legend));
 391
 392   AutoZoom(hClusterEVsTrack,"maxx,maxy",1)->DrawCopy("colz");
 393   outfilename =  QAPATH + "ETrackMult" + fTrigger(r) + ".pdf";
 394   outfilename2 = QAPATH + "ETrackMult" + fTrigger(r) + ".png";
 395   if(SavePlots==2) c3->SaveAs(outfilename);
 396   if(SavePlots) c3->SaveAs(outfilename2);
 397   c3->Write();
 398   delete c3;
 399
 400   TCanvas* c4 = new TCanvas("ClusterEVsV0 ","Correlation E calo Vs V0 signal", 600, 600);
 401   c4->SetLogz();
 402   c4->SetFillColor(0);
 403   c4->SetBorderSize(0);
 404   c4->SetFrameBorderMode(0);
 405
 406   TH2F* hClusterEVsV0S =(TH2F*)outputList->FindObject("EMCAL_hCaloV0SEClusters");
 407   FormatRunHisto(hClusterEVsV0S,Form("Sum E cluster Vs V0 signal%s",legend));
 408
 409   AutoZoom(hClusterEVsV0S,"maxx,maxy",1)->DrawCopy("colz");
 410   outfilename = QAPATH +"EVsV0s" + fTrigger(r) + ".pdf";
 411   outfilename2 = QAPATH +"EVsV0s" + fTrigger(r) + ".png";
 412   if(SavePlots==2) c4->SaveAs(outfilename);
 413   if(SavePlots) c4->SaveAs(outfilename2);
 414   c4->Write();
 415   delete c4;
 416
 417   TCanvas* c5 = new TCanvas("CellsperCluster","Nb of cells per cluster for each SM", 600, 600);
 418   c5->SetLogz();
 419   c5->SetFillColor(0);
 420   c5->SetBorderSize(0);
 421   c5->SetFrameBorderMode(0);
 422   Bool_t mod3=0; if (nSM%3) mod3=1;
 423   c5->Divide(3,(nSM/3)+mod3);
 424
 425   for (int ism = 0; ism < nSM; ism++)
 426     {
 427       c5->cd(ism+1);
 428       gPad->SetLogz();
 429       if(TString(Form("Nb of cells per cluster%s Mod %d",legend,ism)).Length() > 60) { Error(__FUNCTION__,"Title too long!"); return -6;}
 430       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");
 431     }
 432
 433   outfilename =  QAPATH + "CellsperClusterSM" + fTrigger(r) + ".pdf";
 434   outfilename2 = QAPATH + "CellsperClusterSM" + fTrigger(r) + ".png";
 435   if(SavePlots==2) c5->SaveAs(outfilename);
 436   if(SavePlots) c5->SaveAs(outfilename2);
 437   c5->Write();
 438   delete c5;
 439
 440   if (outputList) outputList->Delete();
 441
 442   return 0;
 443 }
 444
 445 //----------------------------------------------------------------------------------------------------------------------------------
 446 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)
 447 {
 448
 449   TString fTrigger="";
 450   TString aCalorimeter;
 451   if (n<=12) {aCalorimeter = fCalorimeter;} else  {aCalorimeter = TString("EMCAL_and_DCAL");}
 452   TDatime now;
 453
 454   Double_t Nevent=0 ;
 455   Double_t xe=0.5;
 456
 457   Double_t CellMean=0;
 458   Double_t CellRMS=0;
 459   Double_t ClusterMean=0;
 460   Double_t ClusterRMS=0;
 461   Double_t EtotalMean=0;
 462   Double_t EtotalRMS=0;
 463
 464   Double_t CellPerClusterMean=0; //
 465   Double_t CellPerClusterRMS=0; //
 466
 467   Double_t mPDG = 134.9766;
 468   Double_t Npi0=0;
 469   Double_t Npi0Err=0;
 470   Double_t MeanPos=0;
 471   Double_t MeanPosErr=0;
 472   Double_t Width=0;
 473   Double_t WidthErr=0;
 474   Double_t Chi2NdfPi0=0;
 475   Double_t Ngg=0;
 476   Double_t NggErr=0;
 477   Double_t Signif=0; // !S/(S+B)
 478   Double_t SignifErr=0; // !S/(S+B)
 479
 480   TFile* ftree = new TFile(Form("%s/trending.root",QAPATH.Data()),"RECREATE");
 481
 482   TTree *tree = new TTree("trending","Trending QA Tree");
 483   tree->Branch("fDate",&now);
 484   tree->Branch("fCalorimeter",&aCalorimeter);
 485   tree->Branch("system",&system);
 486   tree->Branch("period",&period);
 487   tree->Branch("pass",&pass);
 488   tree->Branch("fTrigger",&fTrigger);
 489   tree->Branch("run",&RunId,"run/I");
 490   tree->Branch("xe",&xe,"xe/D");
 491
 492   tree->Branch("Nevent",&Nevent,"Nevent/D");
 493   tree->Branch("CellMean",&CellMean,"CellMean/D");
 494   tree->Branch("CellRMS",&CellRMS,"CellRMS/D");
 495   tree->Branch("ClusterMean",&ClusterMean,"ClusterMean/D");
 496   tree->Branch("ClusterRMS",&ClusterRMS,"ClusterRMS/D");
 497   tree->Branch("EtotalMean",&EtotalMean,"EtotalMean/D");
 498   tree->Branch("EtotalRMS",&EtotalRMS,"EtotalRMS/D");
 499
 500   tree->Branch("CellPerClusterMean",&CellPerClusterMean,"CellPerClusterMean/D"); //
 501   tree->Branch("CellPerClusterRMS",&CellPerClusterRMS,"CellPerClusterRMS/D");  //
 502
 503   tree->Branch("Npi0",&Npi0,"Npi0/D");
 504   tree->Branch("Npi0Err",&Npi0Err,"Npi0Err/D");
 505   tree->Branch("MeanPos",&MeanPos,"MeanPos/D");
 506   tree->Branch("MeanPosErr",&MeanPosErr,"MeanPosErr/D");
 507   tree->Branch("Width",&Width,"Width/D");
 508   tree->Branch("WidthErr",&WidthErr,"WidthErr/D");
 509   tree->Branch("Chi2NdfPi0",&Chi2NdfPi0,"Chi2NdfPi0/D");
 510   tree->Branch("Ngg",&Ngg,"Ngg/D");
 511   tree->Branch("NggErr",&NggErr,"NggErr/D");
 512   tree->Branch("Signif",&Signif,"Signif/D");
 513   tree->Branch("SignifErr",&SignifErr,"SignifErr/D");
 514
 515   tree->Branch("nSM",&n,"nSM/I");
 516
 517   int nMax = 22;
 518   Double_t CellMeanSM[nMax];
 519   Double_t CellRMSSM[nMax];
 520   Double_t ClusterMeanSM[nMax];
 521   Double_t ClusterRMSSM[nMax];
 522   Double_t EtotalMeanSM[nMax]; //mean total energy deposited per event
 523   Double_t EtotalRMSSM[nMax];
 524   Double_t CellPerClusterMeanSM[nMax];
 525   Double_t CellPerClusterRMSSM[nMax];
 526   Double_t ECell1MeanSM[nMax]; //total energy deposited per event without 1 cell clusters
 527   Double_t ECell1RMSSM[nMax];
 528
 529   Double_t MeanPosSM[nMax];
 530   Double_t MeanPosErrSM[nMax];
 531   Double_t WidthSM[nMax];
 532   Double_t WidthErrSM[nMax];
 533   Double_t Npi0SM[nMax];
 534   Double_t Npi0ErrSM[nMax];
 535
 536   tree->Branch("CellMeanSM",CellMeanSM,TString::Format("CellMeanSM[%i]/D",nMax));
 537   tree->Branch("CellRMSSM",CellRMSSM,TString::Format("CellRMSSM[%i]/D",nMax));
 538   tree->Branch("ClusterMeanSM",ClusterMeanSM,TString::Format("ClusterMeanSM[%i]/D",nMax));
 539   tree->Branch("ClusterRMSSM",ClusterRMSSM,TString::Format("ClusterRMSSM[%i]/D",nMax));
 540   tree->Branch("EtotalMeanSM",EtotalMeanSM,TString::Format("EtotalMeanSM[%i]/D",nMax));
 541   tree->Branch("EtotalRMSSM",EtotalRMSSM,TString::Format("EtotalRMSSM[%i]/D",nMax));
 542   tree->Branch("CellPerClusterMeanSM",CellPerClusterMeanSM,TString::Format("CellPerClusterMeanSM[%i]/D",nMax));
 543   tree->Branch("CellPerClusterRMSSM",CellPerClusterRMSSM,TString::Format("CellPerClusterRMSSM[%i]/D",nMax));
 544   tree->Branch("ECell1MeanSM",ECell1MeanSM,TString::Format("ECell1MeanSM[%i]/D",nMax));
 545   tree->Branch("ECell1RMSSM",ECell1RMSSM,TString::Format("ECell1RMSSM[%i]/D",nMax));
 546
 547   tree->Branch("MeanPosSM",MeanPosSM,TString::Format("MeanPosSM[%i]/D",nMax));
 548   tree->Branch("MeanPosErrSM",MeanPosErrSM,TString::Format("MeanPosErrSM[%i]/D",nMax));
 549   tree->Branch("WidthSM",WidthSM,TString::Format("WidthSM[%i]/D",nMax));
 550   tree->Branch("WidthErrSM",WidthErrSM,TString::Format("WidthErrSM[%i]/D",nMax));
 551   tree->Branch("Npi0SM",Npi0SM,TString::Format("Npi0SM[%i]/D",nMax));
 552   tree->Branch("Npi0ErrSM",Npi0ErrSM,TString::Format("Npi0ErrSM[%i]/D",nMax));
 553
 554   TF1* fitMass = new TF1("fitMass",pi0massP2,100,250,6);
 555   fitMass->SetParName(0,"A");
 556   fitMass->SetParName(1,"m_{0}");
 557   fitMass->SetParName(2,"sigma");
 558   fitMass->SetParName(3,"a_{0}");
 559   fitMass->SetParName(4,"a_{1}");
 560   fitMass->SetParName(5,"a_{2}");
 561   fitMass->SetParLimits(0,  1.e-5,1.e5);
 562   fitMass->SetParLimits(1, 0.11, 0.16); //
 563   fitMass->SetParLimits(2,  0.001,0.06);
 564
 565   TList* outputList;
 566
 567   TH1F* fhNEvents;
 568   TH1F* fhE;
 569   TH1F* fhNClusters = 0x0;
 570   TH1F* fhNCells = 0x0;
 571
 572   TH1F* NCells[n];
 573   TH1F* NClusters[n];
 574   TH2F* NCellsPerCluster[n];
 575   TH1F* E[n];
 576
 577   TH2F* fhIM ;
 578   TH1F* fhMgg;
 579   TH2F* IM[n];
 580   TH1F* MggSM[n];
 581
 582   TPMERegexp r("_\\w+");
 583   TIter next(TriggersList);
 584   int ret = 0;
 585   while (TObject *obj = next())
 586     {
 587       fTrigger= TString(obj->GetName());
 588       TString namefile = QAPATH + period + "_" +  pass + fTrigger(r).Data() + "_" + RunId + "_data.txt";
 589       ofstream QAData(namefile, ios::app); // write checks at the end
 590
 591       Npi0=0;
 592       Npi0Err=0;
 593       MeanPos=0;
 594       MeanPosErr=0;
 595       Width=0;
 596       WidthErr=0;
 597       Chi2NdfPi0=0;
 598       Ngg=0;
 599       NggErr=0;
 600       Signif=0;
 601       SignifErr=0;
 602
 603       memset (CellMeanSM, 0, sizeof (Double_t) * nMax);
 604       memset (CellRMSSM, 0, sizeof (Double_t) *  nMax);
 605       memset (ClusterMeanSM, 0, sizeof (Double_t) * nMax);
 606       memset (ClusterRMSSM, 0, sizeof (Double_t) * nMax);
 607       memset (EtotalMeanSM, 0, sizeof (Double_t) * nMax);
 608       memset (EtotalRMSSM, 0, sizeof (Double_t) * nMax);
 609       memset (CellPerClusterMeanSM, 0, sizeof (Double_t) * nMax);
 610       memset (CellPerClusterRMSSM, 0, sizeof (Double_t) * nMax);
 611       memset (ECell1MeanSM, 0, sizeof (Double_t) * nMax);
 612       memset (ECell1RMSSM, 0, sizeof (Double_t) * nMax);
 613
 614       memset (MeanPosSM, 0, sizeof (Double_t) * nMax);
 615       memset (MeanPosErrSM, 0, sizeof (Double_t) * nMax);
 616       memset (WidthSM, 0, sizeof (Double_t) * nMax);
 617       memset (WidthErrSM, 0, sizeof (Double_t) * nMax);
 618       memset (Npi0SM, 0, sizeof (Double_t) * nMax);
 619       memset (Npi0ErrSM, 0, sizeof (Double_t) * nMax);
 620
 621       TString dirname;
 622       if(!fTrigger.Contains("QA")) {
 623         dirname = "CaloQA_";
 624       }
 625       dirname += fTrigger;
 626
 627       Bool_t dirok = f->cd(dirname);
 628       if(!dirok) { Error(__FUNCTION__,Form("No input directory %s",dirname.Data())); tree->Fill(); ftree->cd(); tree->Write(); ret= -1; continue;}
 629       outputList = (TList*)gDirectory->Get(dirname);
 630       if(!outputList){ Error(__FUNCTION__,Form("No input list! %s",dirname.Data())); tree->Fill();  ftree->cd(); tree->Write(); ret=-2; continue;;}
 631       outputList->SetOwner();
 632
 633       // number of events
 634       fhNEvents =(TH1F *)outputList->FindObject("hNEvents");
 635       if(!fhNEvents){ Error(__FUNCTION__,Form("NEvent histogram not found for trigger %s",fTrigger.Data())); tree->Fill();  ftree->cd(); tree->Write();  ret=-3; continue;}
 636       Nevent=fhNEvents->GetEntries();
 637       if(Nevent==0) {Error(__FUNCTION__,Form("No event in trigger %s",fTrigger.Data())); tree->Fill();  ftree->cd(); tree->Write(); ret=-4; continue;}
 638       if(Nevent<20) {Error(__FUNCTION__,Form("Less than 20 events in trigger %s",fTrigger.Data())); tree->Fill();  ftree->cd(); tree->Write(); ret=-5; continue;}
 639
 640       // first do clusters trending
 641       fhE = (TH1F *)outputList->FindObject(fCalorimeter+"_hE");
 642       Double_t energy = 0. ;
 643
 644       for(Int_t ibin = fhE->FindBin(0.6) ; ibin <fhE->FindBin(50.) ; ibin++){
 645         energy+=fhE->GetBinCenter(ibin)*fhE->GetBinContent(ibin);
 646       }
 647       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;}
 648       EtotalMean=energy/fhE->Integral(fhE->FindBin(0.6), fhE->FindBin(50.)) ;
 649       EtotalRMS=fhE->GetMeanError();
 650
 651       TString nameNCell = Form("%s_hNCells_Mod",fCalorimeter.Data());
 652       TString nameNCluster = Form("%s_hNClusters_Mod",fCalorimeter.Data());
 653       TString nameE = Form("%s_hE_Mod",fCalorimeter.Data());
 654       TH2F* hNCellsMod= (TH2F*)outputList->FindObject(nameNCell.Data());
 655       TH2F* hNClusterMod=(TH2F*)outputList->FindObject(nameNCluster.Data());
 656       TH2F* hEMod=(TH2F*)outputList->FindObject(nameE.Data());
 657
 658       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;}
 659
 660       TCanvas* c1 = new TCanvas("Pi0InvMassSM","Pi0 Invariant Mass for each SM", 600, 600);
 661       c1->SetFillColor(0);
 662       c1->SetBorderSize(0);
 663       c1->SetFrameBorderMode(0);
 664       Bool_t mod3=0; if (n%3) mod3=1;
 665       c1->Divide(3,n/3+mod3);
 666
 667       //per sm trending
 668       TString nameNCellPerCluster;
 669       for(Int_t ism = 0 ; ism < n ; ism++){
 670         cout << "#########################"<< endl;
 671         cout      << " Super Module " << ism << " Run " << RunId << endl;
 672         // first do clusters trending
 673         nameNCellPerCluster = Form("%s_hNCellsPerCluster_Mod%d",fCalorimeter.Data(),ism);
 674         NCellsPerCluster[ism] = (TH2F*)outputList->FindObject(nameNCellPerCluster.Data());
 675  if(!  (TH2F*)outputList->FindObject(nameNCellPerCluster.Data()) ) { Error(__FUNCTION__,Form("NCellsPerCluster histogram not found for super module %i",ism));ret=-8; continue;}
 676    NCellsPerCluster[ism] = (TH2F*)outputList->FindObject(nameNCellPerCluster.Data());
 677
 678         NCells[ism] = (TH1F*)hNCellsMod->ProjectionX(Form("NCells%d",ism),ism+1,ism+2,"");
 679         NClusters[ism] = (TH1F*)hNClusterMod->ProjectionX(Form("NClusters%d",ism),ism+1,ism+2,"");
 680         E[ism] = (TH1F*)hEMod->ProjectionX(Form("E%d",ism),ism+1,ism+2,"");
 681         CellMeanSM[ism]=NCells[ism]->GetMean();
 682         CellRMSSM[ism]=NCells[ism]->GetMeanError();
 683         ClusterMeanSM[ism]=NClusters[ism]->GetMean();
 684         ClusterRMSSM[ism]=NClusters[ism]->GetMeanError();
 685         CellPerClusterMeanSM[ism]=NCellsPerCluster[ism]->GetMean(2);
 686         CellPerClusterRMSSM[ism]=NCellsPerCluster[ism]->GetMeanError(2);
 687
 688         ECell1MeanSM[ism] =NCellsPerCluster[ism]->ProjectionX("",2,50,"")->Integral(5,50)/(Nevent);
 689         ECell1RMSSM[ism] =NCellsPerCluster[ism]->ProjectionX("",2,50,"")->GetMeanError();
 690         Double_t energySM = 0. ;
 691         for(Int_t ibin = E[ism]->FindBin(0.6) ; ibin <E[ism]->FindBin(50.) ; ibin++){
 692           energySM+=E[ism]->GetBinCenter(ibin)*(E[ism]->GetBinContent(ibin));
 693         }
 694         if(E[ism]->Integral(E[ism]->FindBin(0.6),E[ism]->FindBin(50.))==0){Error(__FUNCTION__,Form("Energy: Not enough events/SM")); continue;}
 695         EtotalMeanSM[ism]=energySM/(E[ism]->Integral(E[ism]->FindBin(0.6),E[ism]->FindBin(50.)));
 696
 697         EtotalRMSSM[ism]=E[ism]->GetMeanError();
 698
 699         if(ism==0) {
 700           fhNCells = (TH1F*)NCells[ism]->Clone("NCells");
 701           fhNClusters = (TH1F*)NClusters[ism]->Clone("NClusters");
 702         }
 703         else {
 704           fhNCells->Add(NCells[ism],1);
 705           fhNClusters->Add(NClusters[ism],1);
 706         }
 707
 708         //Pi0
 709         c1->cd(ism+1);
 710         TString namePair = Form("%s_hIM_Mod%d",fCalorimeter.Data(),ism);
 711         IM[ism] = (TH2F*)outputList->FindObject(namePair.Data());
 712         IM[ism]->Sumw2();
 713
 714         TString projname = Form("SM_%d",ism);
 715         MggSM[ism] = (TH1F *)IM[ism]->ProjectionY(projname.Data(), 2, 150, "") ;
 716
 717
 718         if(MggSM[ism]->GetEntries()>100) {
 719           fitMass->SetParameter(0, MggSM[ism]->GetBinContent(MggSM[ism]->GetMaximumBin()));
 720           fitMass->SetParameter(1, mPDG); //
 721           fitMass->SetParameter(2, 15.);  //
 722           fitMass->SetParameter(3,0.);
 723           fitMass->SetParameter(4,MggSM[ism]->GetBinContent(MggSM[ism]->FindBin(0.11)));
 724           fitMass->SetParameter(5,MggSM[ism]->GetBinContent(MggSM[ism]->FindBin(0.20)));
 725
 726           if(MggSM[ism]->GetEntries()<1000){ MggSM[ism]->Rebin(4);} else {MggSM[ism]->Rebin();}
 727           MggSM[ism]->Fit("fitMass", "WL R +","",0.05, 0.30);
 728           MggSM[ism]->SetTitle(Form("Pi0 Mass for super module %i",ism));
 729           MggSM[ism]->SetTitleSize(0.1);
 730           MggSM[ism]->SetXTitle("Pi0 Mass");
 731           MggSM[ism]->SetYTitle("Nb of entries");
 732           MggSM[ism]->GetXaxis()->SetLabelSize(0.05);
 733           MggSM[ism]->GetXaxis()->SetTitleSize(0.07);
 734           MggSM[ism]->GetXaxis()->SetTitleOffset(0.68);
 735           MggSM[ism]->GetYaxis()->SetLabelSize(0.05);
 736           MggSM[ism]->GetYaxis()->SetTitleSize(0.06);
 737           MggSM[ism]->GetYaxis()->SetTitleOffset(0.78);
 738
 739           MeanPosSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(1)*1000;
 740           MeanPosErrSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParError(1)*1000;
 741           WidthSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(2)*1000;
 742           WidthErrSM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParError(2)*1000;
 743           Npi0SM[ism] = MggSM[ism]->GetFunction("fitMass")->GetParameter(0)*(MggSM[ism]->GetFunction("fitMass")->GetParameter(2))*TMath::Sqrt(2*TMath::Pi())/(Nevent*MggSM[ism]->GetBinWidth(1));
 744           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))
 745                                        +(MggSM[ism]->GetFunction("fitMass")->GetParError(2)/MggSM[ism]->GetFunction("fitMass")->GetParameter(2))*(MggSM[ism]->GetFunction("fitMass")->GetParError(2)/MggSM[ism]->GetFunction("fitMass")->GetParameter(2)));
 746           Npi0ErrSM[ism] = 0.; //
 747
 748         }// end if enough events for Pi0 fit and trending
 749         else { Info(__FUNCTION__,Form("Not enough events for Pi0 fit and trending for super module %i",ism));} ;
 750       } //per SM loop
 751
 752       // Now Pi0 global trending
 753       TCanvas* c2 = new TCanvas("Pi0InvMass","Pi0 Invariant Mass", 600, 600);
 754       c2->SetFillColor(0);
 755       c2->SetBorderSize(0);
 756       c2->SetFrameBorderMode(0);
 757
 758       fhIM = (TH2F *)outputList->FindObject(fCalorimeter+"_hIM");
 759       fhIM->Sumw2();
 760       fhMgg = (TH1F *)fhIM->ProjectionY("Mgg", 2, 150, "") ;
 761       if(fhMgg->GetEntries()==0) {Error(__FUNCTION__,"The Pi0 histogram is empty !");  tree->Fill(); ret=-8; continue;}
 762       fitMass->SetParameter(0, 4500);
 763       fitMass->SetParameter(1, mPDG);
 764       fitMass->SetParameter(2, 0.01);
 765       fitMass->SetParameter(3,0.);
 766       fitMass->SetParameter(4,fhMgg->GetBinContent(fhMgg->FindBin(0.11)));
 767       fitMass->SetParameter(5,fhMgg->GetBinContent(fhMgg->FindBin(0.20)));
 768
 769       if(fhMgg->GetEntries()<5000){
 770         fhMgg->Rebin(4);}
 771       else   fhMgg->Rebin();
 772
 773       fhMgg->Fit("fitMass", "L R +", "", 0.05, 0.20);
 774
 775       fhMgg->SetTitle("Pi0 Mass");
 776       fhMgg->SetTitleSize(0.1);
 777       fhMgg->SetXTitle("Pi0 Mass");
 778       fhMgg->SetYTitle("Nb of entries");
 779       fhMgg->GetXaxis()->SetLabelSize(0.03);
 780       fhMgg->GetXaxis()->SetTitleSize(0.03);
 781       fhMgg->GetXaxis()->SetTitleOffset(1.3);
 782       fhMgg->GetYaxis()->SetLabelSize(0.03);
 783       fhMgg->GetYaxis()->SetTitleSize(0.03);
 784       fhMgg->GetYaxis()->SetTitleOffset(1.3);
 785
 786       MeanPos = fhMgg->GetFunction("fitMass")->GetParameter(1)*1000;
 787       MeanPosErr = fhMgg->GetFunction("fitMass")->GetParError(1)*1000;
 788       Width = fhMgg->GetFunction("fitMass")->GetParameter(2)*1000;
 789       WidthErr = fhMgg->GetFunction("fitMass")->GetParError(2)*1000;
 790       Chi2NdfPi0 = fhMgg->GetFunction("fitMass")->GetChisquare()/fhMgg->GetFunction("fitMass")->GetNDF();
 791       Npi0 = fhMgg->GetFunction("fitMass")->GetParameter(0)*fhMgg->GetFunction("fitMass")->GetParameter(2)*TMath::Sqrt(2*TMath::Pi())/(Nevent*fhMgg->GetBinWidth(10));
 792       Npi0Err = TMath::Sqrt((fhMgg->GetFunction("fitMass")->GetParError(0)/fhMgg->GetFunction("fitMass")->GetParameter(0))*(fhMgg->GetFunction("fitMass")->GetParError(0)/fhMgg->GetFunction("fitMass")->GetParameter(0))+(WidthErr/Width)*(WidthErr/Width));
 793       Npi0Err = 0.; //
 794       Ngg = fhMgg->GetFunction("fitMass")->Integral(0.11, 0.16)/(Nevent*fhMgg->GetBinWidth(10));
 795       NggErr = fhMgg->GetFunction("fitMass")->IntegralError(0.11, 0.16)/(fhMgg->Integral()*fhMgg->GetBinWidth(10));
 796       Signif = Npi0/Ngg;
 797       SignifErr = TMath::Sqrt((Npi0Err/Npi0*(Npi0Err/Npi0)+(NggErr/Ngg*(NggErr/Ngg))));
 798       SignifErr = Signif*SignifErr;
 799
 800       cout<<"******************"<<endl;
 801       //end of global trending
 802
 803       ClusterMean=fhNClusters->GetMean();
 804       ClusterRMS=fhNClusters->GetMeanError();
 805       CellMean=fhNCells->GetMean();
 806       CellRMS=fhNCells->GetMeanError();
 807       tree->Fill();
 808
 809       TString outfilename = QAPATH +  "Pi0InvMass" + fTrigger(r) + ".pdf";
 810       TString outfilename2 = QAPATH + "Pi0InvMass" + fTrigger(r) + ".png";
 811       if(SavePlots==2) c2->SaveAs(outfilename);
 812       if(SavePlots) c2->SaveAs(outfilename2);
 813
 814       outfilename = QAPATH + "Pi0InvMassSM" + fTrigger(r) + ".pdf";
 815       outfilename2 = QAPATH + "Pi0InvMassSM" + fTrigger(r) + ".png";
 816       if(SavePlots==2) c1->SaveAs(outfilename);
 817       if(SavePlots) c1->SaveAs(outfilename2);
 818
 819       fout->cd();
 820       fout->Cd(Form("%s/%s/%ld/%s/%s",period.Data(),pass.Data(),RunId,"RunLevelQA",fTrigger.Data()));
 821       c2->Write();
 822       c1->Write();
 823       delete c1;
 824       delete c2;
 825       if (outputList) outputList->Delete() ;
 826
 827       QAData  << RunId<<"    "<< Nevent
 828               <<"\n";
 829
 830       QAData.close();
 831
 832     }
 833
 834   ftree->cd();
 835   tree->Write();
 836   ftree->Close();
 837
 838   return ret;
 839
 840 }
 841
 842 //-------------------------------------------------------------------------
 843 TH2F* FormatRunHisto(TH2F* aHisto,const char* title,const char* YTitle)
 844 {
 845
 846   if(!aHisto) {Error(__FUNCTION__,Form("The histogram with title \"%s\" was not found!",title)); return new TH2F();}
 847   aHisto->SetStats(kFALSE);
 848   aHisto->SetTitle(title);
 849   aHisto->SetStats(kFALSE);
 850   aHisto->SetYTitle(YTitle);
 851   aHisto->GetYaxis()->SetTitleOffset(1.2);
 852   aHisto->GetYaxis()->SetLabelSize(0.03);
 853   aHisto->GetZaxis()->SetLabelSize(0.02);
 854
 855   return aHisto;
 856
 857 }
 858
 859 //--------------------------------------------------------------------------------------------------------------
 860 TH2F* HistoPerMod(TH2F* hTmpPerMod,const char* title)
 861 {
 862
 863   if(!hTmpPerMod) {Error(__FUNCTION__,Form("The histogram with title \"%s\" was not found!",title)); return new TH2F();}
 864   hTmpPerMod->SetStats(kFALSE);
 865   hTmpPerMod->SetTitle(title);
 866   hTmpPerMod->SetTitleSize(0.1);
 867   hTmpPerMod->GetXaxis()->SetTitleOffset(1.1);
 868   hTmpPerMod->GetXaxis()->SetTitleSize(0.05);
 869   hTmpPerMod->GetXaxis()->SetLabelSize(0.06);
 870   hTmpPerMod->GetYaxis()->SetTitleOffset(1.1);
 871   hTmpPerMod->GetYaxis()->SetTitleSize(0.05);
 872   hTmpPerMod->GetYaxis()->SetLabelSize(0.06);
 873   hTmpPerMod->GetZaxis()->SetLabelSize(0.04);
 874
 875   return hTmpPerMod;
 876
 877 }
 878
 879 //---------------------------------------------------------------------------------------------------
 880 TH2F* AutoZoom(TH2F* H,Option_t* aType, Int_t EntryMin)
 881 {
 882
 883   Int_t shiftX = (Int_t)(H->GetNbinsX()/30.);
 884   Int_t shiftY = (Int_t)(H->GetNbinsY()/30.);
 885
 886   TString opt = aType;
 887   opt.ToLower();
 888
 889   int minX = 0;
 890   int maxX = H->GetNbinsX();
 891   int New_minX = minX;
 892   int New_maxX = maxX;
 893
 894   int minY = 0;
 895   int maxY = H->GetNbinsY();
 896   int New_minY = minY;
 897   int New_maxY = maxY;
 898
 899   if (opt.Contains("all")) opt = TString("minx,maxx,miny,maxy");
 900
 901   if (opt.Contains("maxx"))
 902     {
 903
 904       for  (New_maxX = maxX;New_maxX >=minX; New_maxX--)
 905         {  Stat_t c = 0;
 906           for  (int i_y = maxY; i_y >= minY;i_y--)
 907             { c = H->GetBinContent(New_maxX,i_y);  if (c>EntryMin) break;}
 908           if (c>EntryMin) break;
 909         }
 910     }
 911
 912   if (opt.Contains("maxy"))
 913     {
 914
 915       for  (New_maxY = maxY;New_maxY >=minY;New_maxY--)
 916         {  Stat_t c = 0;
 917           for  (int i_x=maxX; i_x>=minX;i_x--)
 918             { c = H->GetBinContent(i_x, New_maxY );  if (c>EntryMin) break;}
 919           if (c>EntryMin) break;
 920         }
 921
 922     }
 923
 924   if (opt.Contains("minx"))
 925     {
 926
 927       for  (New_minX = minX;New_minX <=maxX; New_minX++)
 928         {  Stat_t c = 0;
 929           for  (int i_y = minY; i_y <= maxY;i_y++)
 930             { c = H->GetBinContent(New_minX,i_y);  if (c>EntryMin) break;}
 931           if (c>EntryMin) break;
 932         }
 933     }
 934
 935   if (opt.Contains("miny"))
 936     {
 937       for  (New_minY = minY;New_minY <=maxY;New_minY++)
 938         {  Stat_t c = 0;
 939           for  (int i_x=minX; i_x<=maxX;i_x++)
 940             { c = H->GetBinContent(i_x, New_minY );  if (c>EntryMin) break;}
 941           if (c>EntryMin) break;
 942         }
 943     }
 944
 945  if (New_maxX!=-1 && New_maxY!=-1) H->GetXaxis()->SetRange(New_minX - shiftX  , New_maxX + shiftX);
 946  if (New_maxX!=-1 && New_maxY!=-1) H->GetYaxis()->SetRange(New_minY - shiftY  , New_maxY + shiftY);
 947
 948   return H;
 949
 950 }
 951
 952 //----------------------------------------------------------------------------------------------------
 953 int FindNumberOfSM(TFile* f, TString fTrigger, TString period)
 954 {
 955
 956   TString direct;
 957   if(!fTrigger.Contains("QA")) {
 958     direct = "CaloQA_";
 959   }
 960   direct += fTrigger;
 961
 962   Int_t nSMt=-1;
 963   Int_t year = 2000 + TString(period(3,2)).Atoi();
 964   if ( year == 2010 ) { nSMt=6; }
 965   else if ( year == 2011 || year == 2012 ) { nSMt=10; }
 966   else if ( year == 2013 || year == 2014 ) { nSMt=12; }
 967   else { nSMt=20; }
 968
 969   TList *outputList;
 970   Bool_t dirok = f->cd(direct);
 971   if (!dirok) { Error(__FUNCTION__,Form("No input directory %s, the number SMs will be returned based on the year!",direct.Data()));}
 972   else { outputList = (TList*)gDirectory->Get(direct);}
 973   if(!outputList) { Error(__FUNCTION__,"No input list, the number SMs will be returned based on the year! ");}
 974   else {
 975     outputList->SetOwner();
 976     TH2F* hNSM =(TH2F *)outputList->FindObject("EMCAL_hE_Mod");
 977     if (!hNSM || (!hNSM->GetEntries())) { Error(__FUNCTION__,"hNSM Histogram not found or it is empty, the number SMs will be returned based on the year!");}
 978     else {
 979       nSMt = hNSM->GetYaxis()->GetBinUpEdge(hNSM->FindLastBinAbove(0,2));
 980     }
 981   }
 982   if (outputList) {outputList->Delete();}
 983
 984   return nSMt;
 985
 986 }