PWGHF/vertexingHF/macros/DrawQAoutput.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include <TFile.h>
   4 #include <TString.h>
   5 #include <TH2F.h>
   6 #include <TH1F.h>
   7 #include <TF1.h>
   8 #include <TGraph.h>
   9 #include <TDirectoryFile.h>
  10 #include <TList.h>
  11 #include <TCanvas.h>
  12 #include <TLegend.h>
  13 #include <TPaveText.h>
  14 #include <TPaveStats.h>
  15 #include <TStyle.h>
  16 #include <TClass.h>
  17 #include <TDatabasePDG.h>
  18 #include <TParameter.h>
  19 #include <AliCounterCollection.h>
  20 #include <AliRDHFCuts.h>
  21 #endif
  22
  23 /* $Id$ */
  24
  25 TString *periodsname;
  26
  27 //read the file and take list and stat
  28
  29 Bool_t ReadFile(TList* &list,TH1F* &hstat, TString listname,TString partname,TString path="./",TString filename=/*"AnalysisResults.root"*/"PWG3histograms.root", TString dirname="PWG3_D2H_QA");
  30 Bool_t ReadFileMore(TList* &list,TH1F* &hstat, AliRDHFCuts* &cutobj, TString listname,TString partname,TString path="./",TString filename=/*"AnalysisResults.root"*/"PWG3histograms.root", TString dirname="PWG3_D2H_QA");
  31 void SuperimposeBBToTPCSignal(Int_t period /*0=LHC10bc, 1=LHC10d, 2=LHC10h*/,TCanvas* cpid, Int_t set);
  32 void TPCBetheBloch(Int_t set);
  33
  34 Bool_t ReadFile(TList* &list,TH1F* &hstat, TString listname,TString partname,TString path,TString filename, TString dirname){
  35
  36   TString hstatname="nEntriesQA", cutobjname="";
  37   filename.Prepend(path);
  38   listname+=partname;
  39   hstatname+=partname;
  40
  41   TFile* f=new TFile(filename.Data());
  42   if(!f->IsOpen()){
  43     cout<<filename.Data()<<" not found"<<endl;
  44     return kFALSE;
  45   }
  46   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname);
  47   if(!dir){
  48     cout<<dirname.Data()<<" not found in "<<filename.Data()<<endl;
  49     f->ls();
  50     return kFALSE;
  51   }
  52
  53   list=(TList*)dir->Get(listname);
  54   if(!list){
  55     cout<<"List "<<listname.Data()<<" not found"<<endl;
  56     dir->ls();
  57     return kFALSE;
  58   }
  59
  60   hstat=(TH1F*)dir->Get(hstatname);
  61   if(!hstat){
  62     cout<<hstatname.Data()<<" not found"<<endl;
  63     return kFALSE;
  64   }
  65
  66   return kTRUE;
  67 }
  68
  69 Bool_t ReadFileMore(TList* &list,TH1F* &hstat, AliRDHFCuts* &cutobj, TString listname,TString partname,TString path,TString filename,TString dirname){
  70
  71   TString hstatname="nEntriesQA", cutobjname="";
  72   filename.Prepend(path);
  73   listname+=partname;
  74   hstatname+=partname;
  75
  76   if(partname.Contains("Dplus")) cutobjname="AnalysisCuts";//"DplustoKpipiCutsStandard";
  77   else{
  78     if(partname.Contains("D0")) cutobjname="D0toKpiCutsStandard";//"D0toKpiCuts"
  79     else{
  80       if(partname.Contains("Dstar")) cutobjname="DStartoKpipiCuts";
  81       else{
  82         if(partname.Contains("Ds")) cutobjname="DstoKKpiCuts";
  83         else{
  84           if(partname.Contains("D04")) cutobjname="D0toKpipipiCuts";
  85           else{
  86             if(partname.Contains("Lc")) cutobjname="LctopKpiAnalysisCuts";
  87           }
  88         }
  89       }
  90     }
  91   }
  92
  93   TFile* f=new TFile(filename.Data());
  94   if(!f->IsOpen()){
  95     cout<<filename.Data()<<" not found"<<endl;
  96     return kFALSE;
  97   }
  98   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname);
  99   if(!dir){
 100     cout<<dirname.Data()<<" not found  in "<<filename.Data()<<endl;
 101     return kFALSE;
 102   }
 103
 104   list=(TList*)dir->Get(listname);
 105   if(!list){
 106     cout<<"List "<<listname.Data()<<" not found"<<endl;
 107     dir->ls();
 108     return kFALSE;
 109   }
 110
 111   hstat=(TH1F*)dir->Get(hstatname);
 112   if(!hstat){
 113     cout<<hstatname.Data()<<" not found"<<endl;
 114     return kFALSE;
 115   }
 116
 117   cutobj=(AliRDHFCuts*)dir->Get(cutobjname);
 118   if(!cutobj){
 119     cout<<cutobjname.Data()<<" not found"<<endl;
 120     return kFALSE;
 121   }
 122
 123   return kTRUE;
 124 }
 125
 126 //draw "track related" histograms (list "outputTrack")
 127 void DrawOutputTrack(TString partname="D0",TString textleg="",TString path="./", Bool_t superimpose=kFALSE, TString suffixdir="",TString filename=/*"AnalysisResults.root"*/"PWG3histograms.root", Bool_t normint=kTRUE /*normalize at integral, or at nevents*/){
 128   gStyle->SetCanvasColor(0);
 129   gStyle->SetTitleFillColor(0);
 130   gStyle->SetStatColor(0);
 131   gStyle->SetPalette(1);
 132
 133   TString listname="outputTrack",name1="",name2="",path2="",filename2="PWG3histograms.root",legtext1="",legtext2="",suffixdir2="";
 134   TString tmp="y";
 135
 136   if(superimpose){
 137     cout<<"Enter the names:\n>";
 138     cin>>name1;
 139     cout<<">";
 140     cin>>name2;
 141     cout<<"Are they in the same output file? (y/n)"<<endl;
 142     cin>>tmp;
 143     if(tmp=="n"){
 144       cout<<"Path: \n";
 145       cout<<">";
 146       cin>>path2;
 147       cout<<"Filename: "<<endl;
 148       cout<<">";
 149       cin>>filename2;
 150       cout<<"Suffix in the dir name, if any (otherwhise write no)\n>";
 151       cin>>suffixdir2;
 152       if(suffixdir2=="no") suffixdir2="";
 153       cout<<"Text in the legend:\n 1)";
 154       cin>>legtext1;
 155       cout<<"2)";
 156       cin>>legtext2;
 157     }
 158
 159   }
 160
 161   Float_t nevents,nevents2;
 162   TList* list;
 163   TH1F * hstat;
 164   TString dirname="PWG3_D2H_QA";
 165   //dirname+=suffixdir;
 166   Bool_t isRead=ReadFile(list,hstat,listname,Form("%s%s",partname.Data(),name1.Data()),path,filename,Form("%s%s",dirname.Data(),suffixdir.Data()));
 167   if(!isRead) return;
 168   if(!list || !hstat){
 169     cout<<":-( null pointers..."<<endl;
 170     return;
 171   }
 172   nevents=hstat->GetBinContent(1);
 173   TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 174   pvtxt->SetBorderSize(0);
 175   pvtxt->SetFillStyle(0);
 176   pvtxt->AddText(legtext1);
 177
 178   TList* llist;
 179   TH1F* hhstat;
 180   if(superimpose){
 181     isRead=ReadFile(llist,hhstat,listname,Form("%s%s",partname.Data(),name2.Data()),path2,filename2,Form("%s%s",dirname.Data(),suffixdir2.Data()));
 182     if(!isRead) return;
 183     if(!llist || !hhstat){
 184       cout<<":-( null pointers..."<<endl;
 185       return;
 186     }
 187     nevents2=hhstat->GetBinContent(1);
 188     TText *redtext=pvtxt->AddText(legtext2);
 189     redtext->SetTextColor(kRed);
 190     hhstat->Scale(hstat->Integral()/hhstat->Integral());
 191
 192   }
 193
 194   for(Int_t i=0;i<list->GetEntries();i++){
 195     TH1F* h=(TH1F*)list->At(i);
 196     TH1F* hh=0x0;
 197     TH1F* hr=0x0;
 198     if(superimpose){
 199       hh=(TH1F*)llist->At(i);
 200       hr=(TH1F*)hh->Clone(Form("%s_ratio",hh->GetName()));
 201       if(hh && TMath::Abs(hh->Integral()) > 1e-6)hh->Scale(h->Integral()/hh->Integral());
 202     }
 203     if(!h || (superimpose && !hh)){
 204       cout<<"Histogram "<<i<<" not found"<<endl;
 205       continue;
 206     }
 207     if(superimpose){
 208       if(normint) hh->Scale(h->Integral()/hh->Integral());
 209       else hh->Scale(nevents/nevents2);
 210       hhstat->SetLineColor(kRed);
 211       hh->SetLineColor(kRed);
 212       hr->Divide(h);
 213     }
 214
 215     TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
 216     c->cd();
 217     c->SetGrid();
 218     TString hname=h->GetName();
 219     if(!hname.Contains("nCls")){
 220       c->SetLogy();
 221       if(hname.Contains("Layer")){
 222         for(Int_t ibin=1;ibin<=h->GetNbinsX();ibin++){
 223           h->GetXaxis()->SetBinLabel(ibin+1,Form("%d",ibin));
 224         }
 225         h->GetXaxis()->SetLabelSize(0.06);
 226         h->GetXaxis()->SetRangeUser(0,6); //comment to see ntracks!
 227       }
 228       //h->SetMinimum(1);
 229       h->Draw();
 230       if(superimpose)
 231         {
 232           hh->Draw("sames");
 233           TCanvas* c2=new TCanvas(Form("c2%s",h->GetName()),h->GetName());
 234           c2->cd();
 235           c2->SetGrid();
 236           hr->Draw();
 237           c2->SaveAs(Form("%s%s%s%sRatio.png",c->GetName(),name1.Data(),name2.Data(),textleg.Data()));
 238
 239         }
 240     } else {
 241       h->Draw("htext0");
 242       if(superimpose)hh->Draw("htext0sames");
 243     }
 244     c->cd();
 245     pvtxt->Draw();
 246     c->SaveAs(Form("%s%s%s%s.png",c->GetName(),name1.Data(),name2.Data(),textleg.Data()));
 247     c->SaveAs(Form("%s%s%s%s.eps",c->GetName(),name1.Data(),name2.Data(),textleg.Data()));
 248   }
 249
 250   TCanvas* cst=new TCanvas("cst","Stat");
 251   cst->SetGridy();
 252   cst->cd();
 253   hstat->Draw("htext0");
 254   if(superimpose) {
 255     hhstat->Draw("htext0sames");
 256     pvtxt->Draw();
 257   }
 258   cst->SaveAs(Form("%s%s.png",hstat->GetName(),textleg.Data()));
 259   cst->SaveAs(Form("%s%s.eps",hstat->GetName(),textleg.Data()));
 260
 261   TH1F* hd0fb4=(TH1F*)list->FindObject("hd0TracksFilterBit4");
 262   TH1F* hd0SPD1=(TH1F*)list->FindObject("hd0TracksSPDin");
 263   TH1F* hd0SPDany=(TH1F*)list->FindObject("hd0TracksSPDany");
 264   TH1F* hd0TPCITScuts=(TH1F*)list->FindObject("hd0TracksTPCITSSPDany");
 265   TH1F* hhd0fb4=0x0; TH1F* hhd0SPD1=0x0; TH1F* hhd0SPDany=0x0; TH1F* hhd0TPCITScuts=0x0;
 266   if(superimpose){
 267     hhd0fb4=(TH1F*)llist->FindObject("hd0TracksFilterBit4");
 268     hhd0SPD1=(TH1F*)llist->FindObject("hd0TracksSPDin");
 269     hhd0SPDany=(TH1F*)llist->FindObject("hd0TracksSPDany");
 270     hhd0TPCITScuts=(TH1F*)llist->FindObject("hd0TracksTPCITSSPDany");
 271
 272   }
 273   if(hd0fb4 && hd0SPD1 && hd0SPDany && hd0TPCITScuts){
 274     TCanvas* ctrsel=new TCanvas("ctrsel","Track Sel");
 275     ctrsel->SetLogy();
 276     hd0SPD1->SetLineColor(kCyan+3);
 277     hd0SPD1->Draw();
 278     ctrsel->Update();
 279     TPaveStats *st1=(TPaveStats*)hd0SPD1->GetListOfFunctions()->FindObject("stats");
 280     st1->SetTextColor(kCyan+3);
 281     st1->SetY1NDC(0.71);
 282     st1->SetY2NDC(0.9);
 283     hd0SPDany->SetLineColor(4);
 284     hd0SPDany->Draw("sames");
 285     ctrsel->Update();
 286     TPaveStats *st2=(TPaveStats*)hd0SPDany->GetListOfFunctions()->FindObject("stats");
 287     st2->SetY1NDC(0.51);
 288     st2->SetY2NDC(0.7);
 289     st2->SetTextColor(4);
 290     hd0fb4->SetLineColor(2);
 291     hd0fb4->Draw("sames");
 292     ctrsel->Update();
 293     TPaveStats *st3=(TPaveStats*)hd0fb4->GetListOfFunctions()->FindObject("stats");
 294     st3->SetY1NDC(0.31);
 295     st3->SetY2NDC(0.5);
 296     st3->SetTextColor(2);
 297     hd0TPCITScuts->SetLineColor(kGreen+1);
 298     hd0TPCITScuts->Draw("sames");
 299     ctrsel->Update();
 300     TPaveStats *st4=(TPaveStats*)hd0TPCITScuts->GetListOfFunctions()->FindObject("stats");
 301     st4->SetY1NDC(0.71);
 302     st4->SetY2NDC(0.9);
 303     st4->SetX1NDC(0.55);
 304     st4->SetX2NDC(0.75);
 305     st4->SetTextColor(kGreen+1);
 306     ctrsel->Modified();
 307     TLegend* leg=new TLegend(0.15,0.5,0.45,0.78);
 308     leg->SetFillStyle(0);
 309     leg->SetBorderSize(0);
 310     leg->AddEntry(hd0SPD1,"kITSrefit+SPD inner","L");
 311     leg->AddEntry(hd0SPDany,"kITSrefit+SPD any","L");
 312     leg->AddEntry(hd0TPCITScuts,"TPC+ITS cuts+SPD any","L");
 313     leg->AddEntry(hd0fb4,"Filter Bit 4","L");
 314     leg->Draw();
 315
 316     if(superimpose && hhd0fb4 && hhd0SPD1 && hhd0SPDany && hhd0TPCITScuts){
 317       hhd0SPD1->SetLineColor(kCyan+3);
 318       hhd0SPD1->SetStats(0);
 319       hhd0SPD1->SetLineStyle(3);
 320       hhd0SPDany->SetLineColor(4);
 321       hhd0SPDany->SetStats(0);
 322       hhd0SPDany->SetLineStyle(3);
 323       hhd0TPCITScuts->SetLineColor(kGreen+1);
 324       hhd0TPCITScuts->SetStats(0);
 325       hhd0TPCITScuts->SetLineStyle(3);
 326       hhd0fb4->SetLineColor(2);
 327       hhd0fb4->SetStats(0);
 328       hhd0fb4->SetLineStyle(3);
 329       ctrsel->cd();
 330       hhd0SPD1->Draw("sames");
 331       hhd0SPDany->Draw("sames");
 332       hhd0TPCITScuts->Draw("sames");
 333       hhd0fb4->Draw("sames");
 334
 335     }
 336     ctrsel->SaveAs("ImpactParameterTrackSel.eps");
 337     ctrsel->SaveAs("ImpactParameterTrackSel.png");
 338
 339   }
 340 }
 341
 342 //draw "pid related" histograms (list "outputPID")
 343 //period=-999 to draw the pull instead of the cut
 344 void DrawOutputPID(TString partname="D0", Int_t mode=0/*0=with pull, 1=with nsigma*/,TString textleg="",TString path="./",TString suffixdir="", TString filename="AnalysisResults.root"){
 345   gStyle->SetCanvasColor(0);
 346   gStyle->SetTitleFillColor(0);
 347   gStyle->SetOptStat(0);
 348   gStyle->SetPalette(1);
 349
 350   Int_t period=2 ,set=0;
 351   if(mode==1){
 352     cout<<"Choose period: \n-LHC10h -> 2;\n-LHC10de -> 1;\n-LHC10bc -> 0"<<endl;
 353     cin>>period;
 354     if(period>0){
 355       cout<<"Choose set: "<<endl;
 356       if(period==2) cout<<"-pass1 -> 0;\n-pass2 -> 1"<<endl;
 357       cin>>set;
 358     }
 359   }
 360
 361   TString listname="outputPid";
 362   TString dirname="PWG3_D2H_QA";
 363   dirname+=suffixdir;
 364
 365   TList* list;
 366   TH1F * hstat;
 367   //needed only for mode 1
 368   AliRDHFCuts* cutobj;
 369   AliAODPidHF* aodpid;
 370   Double_t nsigmaTOF=0;
 371   Double_t nsigmaTPC[3]={},plimTPC[2]={};
 372
 373   if(mode==1){
 374     Bool_t isRead=ReadFileMore(list,hstat,cutobj,listname,partname,path,filename,dirname);
 375     if(!isRead) return;
 376     if(!list || !hstat){
 377       cout<<":-( null pointers..."<<endl;
 378       return;
 379     }
 380     aodpid=(AliAODPidHF*)cutobj->GetPidHF();
 381     if(!aodpid){
 382       cout<<"PidHF object not found! cannot get the nsigma values"<<endl;
 383       return;
 384     }
 385     nsigmaTOF=aodpid->GetSigma(3);
 386
 387     nsigmaTPC[0]=aodpid->GetSigma(0);
 388     nsigmaTPC[1]=aodpid->GetSigma(1);
 389     nsigmaTPC[2]=aodpid->GetSigma(2);
 390     aodpid->GetPLimit(plimTPC);
 391
 392   }else{
 393     Bool_t isRead=ReadFile(list,hstat,listname,partname,path,filename,dirname);
 394     if(!isRead) return;
 395     if(!list || !hstat){
 396       cout<<":-( null pointers..."<<endl;
 397       return;
 398     }
 399   }
 400
 401
 402   TPaveText *txtsigmaTOF=new TPaveText(0.1,0.65,0.5,0.9,"NDC");
 403   txtsigmaTOF->SetBorderSize(0);
 404   txtsigmaTOF->SetFillStyle(0);
 405   txtsigmaTOF->AddText(Form("nsigmacut from cutobj = %.1f",nsigmaTOF));
 406   TLine lTOF;
 407   lTOF.SetLineColor(kMagenta+1);
 408   lTOF.SetLineStyle(2);
 409   lTOF.SetLineWidth(3);
 410
 411   TPaveText *txtsigmaTPC=new TPaveText(0.3,0.6,0.6,0.9,"NDC");
 412   txtsigmaTPC->SetBorderSize(0);
 413   txtsigmaTPC->SetFillStyle(0);
 414   txtsigmaTPC->AddText("nsigmacut from cutobj \n");
 415   txtsigmaTPC->AddText(Form("p < %.1f : %.1f \n",plimTPC[0],nsigmaTPC[0]));
 416   txtsigmaTPC->AddText(Form("%.1f < p < %.1f : %.1f \n",plimTPC[0],plimTPC[1],nsigmaTPC[1]));
 417   txtsigmaTPC->AddText(Form("p > %.1f : %.1f \n",plimTPC[1],nsigmaTPC[2]));
 418   TLine lTPC;
 419   lTPC.SetLineColor(kMagenta+1);
 420   lTPC.SetLineStyle(2);
 421   lTPC.SetLineWidth(3);
 422
 423   // TCanvas *ctest=new TCanvas("text","Test text");
 424   // ctest->cd();
 425   // txtsigmaTPC->Draw();
 426   // txtsigmaTOF->Draw();
 427
 428
 429   for(Int_t i=0;i<list->GetEntries();i++){
 430     TClass* objtype=list->At(i)->IsA();
 431     TString tpname=objtype->GetName();
 432
 433     if(tpname=="TH1F"){
 434       TH1F* h=(TH1F*)list->At(i);
 435
 436       if(!h){
 437         cout<<"Histogram "<<i<<" not found"<<endl;
 438         continue;
 439       }
 440       //h->Scale(1./h->Integral("width"));
 441       TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
 442       c->SetLogz();
 443       c->cd();
 444       h->Draw();
 445
 446       //write
 447       c->SaveAs(Form("%s%s.png",h->GetName(),textleg.Data()));
 448       c->SaveAs(Form("%s%s.eps",h->GetName(),textleg.Data()));
 449       TFile* fout=new TFile(Form("%s.root",h->GetName()),"recreate");
 450       fout->cd();
 451       c->Write();
 452     }
 453
 454     if(tpname=="TH2F"){
 455       TH2F* h=(TH2F*)list->At(i);
 456
 457       if(!h){
 458         cout<<"Histogram "<<i<<" not found"<<endl;
 459         continue;
 460       }
 461       TString hname=h->GetName();
 462       h->Sumw2();
 463       if(h->Integral("width")==0) {cout<<"Empty histogram, skip\n"; continue;}
 464       h->Scale(1./h->Integral("width"));
 465
 466       Double_t maxzaxis=h->GetBinContent(h->GetMaximumBin());
 467       Double_t minzaxis=h->GetBinContent(h->GetMinimumBin());
 468       printf("Minimum = %f, maximum = %f\n",minzaxis,maxzaxis);
 469       TH2F* hallzrange=(TH2F*)h->Clone(Form("%swholez",hname.Data()));
 470       hallzrange->SetAxisRange(1e-07,maxzaxis,"Z");
 471       //hallzrange->SetAxisRange(minzaxis,maxzaxis,"Z");
 472
 473       TCanvas* cwholez=new TCanvas(Form("c%swholez",hname.Data()),Form("%s down to lowest z",hname.Data()));
 474       cwholez->SetLogz();
 475       hallzrange->Draw("colz");
 476       cwholez->SaveAs(Form("%swholez.png",h->GetName()));
 477       cwholez->SaveAs(Form("%swholez.eps",h->GetName()));
 478
 479       if(hname.Contains("hTOFtimeKaonHyptime")){
 480         TCanvas* cz=new TCanvas(Form("c%szoom",hname.Data()),Form("%szoom",hname.Data()));
 481         cz->SetLogz();
 482         TH2F* hz=(TH2F*)h->Clone(Form("%sz",hname.Data()));
 483         hz->Draw("colz");
 484         hz->SetAxisRange(-1500,1500,"Y");
 485         hz->SetAxisRange(0.,5.,"X");
 486         //write
 487         cz->SaveAs(Form("%szoom.png",h->GetName()));
 488         cz->SaveAs(Form("%szoom.eps",h->GetName()));
 489       }
 490
 491       TCanvas* c=new TCanvas(Form("c%s",hname.Data()),hname.Data());
 492       c->SetLogz();
 493       //c->SetLogx();
 494       TCanvas* c2=new TCanvas(Form("c2%s",hname.Data()),hname.Data());
 495       c2->SetLogz();
 496
 497       c->cd();
 498       h->DrawClone("colz");
 499
 500       if (hname.Contains("Sig") || hname.Contains("sigma"))h->SetAxisRange(-5,5,"Y");
 501       c2->cd();
 502       //if (hname.Contains("TOFtime"))h->SetAxisRange(-1500,1500,"Y");
 503       h->SetAxisRange(0.,5.,"X");
 504
 505       h->Draw("colz");
 506
 507       //TCanvas *test=new TCanvas("test","test");
 508       if(mode==0){
 509         //mean and pull, code from Jens Wiechula
 510         TF1 fg("fg","gaus",-2.,2.); // fit range +- 2 sigma
 511         TLine l;
 512         TObjArray arr;
 513
 514         //h->Draw("colz");
 515         fg.SetParameters(1,0,1);
 516         h->FitSlicesY(&fg,0,-1,0,"NQR",&arr);
 517
 518         TH1 *hM=(TH1*)arr.At(1);
 519         hM->SetMarkerStyle(20);
 520         hM->SetMarkerSize(.5);
 521         hM->DrawClone("sames");
 522
 523         TH1 *hS=(TH1*)arr.At(2);
 524         hS->SetMarkerStyle(20);
 525         hS->SetMarkerSize(.5);
 526         hS->SetMarkerColor(kRed);
 527         hS->SetLineColor(kRed);
 528         hS->DrawClone("same");
 529
 530         l.SetLineColor(kBlack);
 531         l.DrawLine(.2,0,20,0);
 532         l.SetLineColor(kRed);
 533         l.DrawLine(.2,1,20,1);
 534
 535       }else{ //mode 1
 536
 537         if(hname.Contains("TOFsigma")) {
 538
 539           c->cd();
 540           txtsigmaTOF->Draw();
 541           lTOF.DrawLine(.2,nsigmaTOF,20,nsigmaTOF);
 542           lTOF.DrawLine(.2,-1*nsigmaTOF,4.,-1*nsigmaTOF);
 543
 544         }
 545
 546
 547         if(hname.Contains("TPCsigma")){
 548
 549           c->cd();
 550           txtsigmaTPC->Draw();
 551           lTPC.DrawLine(0.,nsigmaTPC[0],plimTPC[0],nsigmaTPC[0]);
 552           lTPC.DrawLine(plimTPC[0],nsigmaTPC[1],plimTPC[1],nsigmaTPC[1]);
 553           lTPC.DrawLine(plimTPC[1],nsigmaTPC[2],4,nsigmaTPC[2]);
 554           lTPC.DrawLine(0.,-1*nsigmaTPC[0],plimTPC[0],-1*nsigmaTPC[0]);
 555           lTPC.DrawLine(plimTPC[0],-1*nsigmaTPC[1],plimTPC[1],-1*nsigmaTPC[1]);
 556           lTPC.DrawLine(plimTPC[1],-1*nsigmaTPC[2],4,-1*nsigmaTPC[2]);
 557         }
 558
 559         if(hname.Contains("TPCsigvsp")){
 560           SuperimposeBBToTPCSignal(period,c,set);
 561         }
 562       }
 563
 564       //write
 565       c->SaveAs(Form("%s%d.png",h->GetName(),mode));
 566       c->SaveAs(Form("%s%d.eps",h->GetName(),mode));
 567       c2->SaveAs(Form("%s2%d.png",h->GetName(),mode));
 568       c2->SaveAs(Form("%s2%d.eps",h->GetName(),mode));
 569
 570       TFile* fout=new TFile(Form("%s%d.root",h->GetName(),mode),"recreate");
 571       fout->cd();
 572       c->Write();
 573       c2->Write();
 574     }
 575   }
 576 }
 577
 578 void SuperimposeBBToTPCSignal(Int_t period /*0=LHC10bc, 1=LHC10d, 2=LHC10h*/,TCanvas* cpid,Int_t set /*see below*/){
 579
 580   TFile* fBethe=new TFile("BetheBlochTPC.root");
 581   if(!fBethe->IsOpen()){
 582     TPCBetheBloch(set);
 583     fBethe=new TFile("BetheBlochTPC.root");
 584   }
 585   const Int_t npart=4;
 586   TString partnames[npart]={"Kaon","Pion","Electron","Proton"};
 587   for(Int_t ipart=0;ipart<npart;ipart++){
 588     TString grname=Form("%sP%d",partnames[ipart].Data(),period);
 589     TGraph* gr=(TGraph*)fBethe->Get(grname);
 590     cpid->cd();
 591     gr->SetLineColor(1);
 592     gr->SetLineWidth(2);
 593     gr->Draw("L");
 594   }
 595
 596   //cpid->SaveAs(Form("%sBB.png",hname.Data()));
 597 }
 598
 599 //draw and save Bethe Bloch from TPC in different periods
 600 void TPCBetheBloch(Int_t set){
 601   gStyle->SetOptTitle(0);
 602   gStyle->SetCanvasColor(0);
 603
 604   AliTPCPIDResponse *tpcResp=new AliTPCPIDResponse();
 605
 606   const Int_t npart=4;
 607   //Double_t masses[npart]={TDatabasePDG::Instance()->GetParticle(321)->Mass()/*Kaon*/,TDatabasePDG::Instance()->GetParticle(211)->Mass()/*Pion*/,TDatabasePDG::Instance()->GetParticle(11)->Mass()/*Electron*/,TDatabasePDG::Instance()->GetParticle(2212)->Mass()/*Proton*/};
 608   TString partnames[npart]={"Kaon","Pion","Electron","Proton"};
 609   //printf("%s = %.4f,%s = %.4f,%s = %.4f\n",partnames[0].Data(),masses[0],partnames[1].Data(),masses[1],partnames[2].Data(),masses[2]);
 610   TCanvas *cBethe=new TCanvas("cBethe","Bethe Bloch K pi e p");
 611   Int_t nperiods=4; //LHC10b+c, LHC10d, LHC10h, MC
 612   Double_t alephParameters[5]={};
 613   Int_t nsets=1/*LHC10bc*/+2/*LHC10de*/+2/*LHC10h*/+3/*MC*/;
 614
 615   periodsname=new TString[nsets];
 616   cout<<"Creating the file of the Bethe Bloch"<<endl;
 617   TFile* fout=new TFile("BetheBlochTPC.root","recreate");
 618
 619   for(Int_t iperiod=0;iperiod<nperiods;iperiod++){
 620     cout<<"Period "<<iperiod<<" : ";
 621     if(iperiod==0){ //LHC10bc
 622
 623       alephParameters[0] = 0.0283086/0.97;
 624       alephParameters[1] = 2.63394e+01;
 625       alephParameters[2] = 5.04114e-11;
 626       alephParameters[3] = 2.12543e+00;
 627       alephParameters[4] = 4.88663e+00;
 628       periodsname[0]="dataLHC10bc";
 629     }
 630     if(iperiod==1){ //LHC10de,low energy
 631       if(set==0){
 632         alephParameters[0] = 1.63246/50.;
 633         alephParameters[1] = 2.20028e+01;
 634         alephParameters[2] = TMath::Exp(-2.48879e+01);
 635         alephParameters[3] = 2.39804e+00;
 636         alephParameters[4] = 5.12090e+00;
 637         periodsname[1]="dataLHC10deold";
 638       }
 639       if(set==1){
 640         alephParameters[0] = 1.34490e+00/50.;
 641         alephParameters[1] =  2.69455e+01;
 642         alephParameters[2] =  TMath::Exp(-2.97552e+01);
 643         alephParameters[3] = 2.35339e+00;
 644         alephParameters[4] = 5.98079e+00;
 645         periodsname[2]="dataLHC10denew";
 646       }
 647     }
 648
 649     if(iperiod==2){//LHC10h
 650       if(set==0){//pass1
 651         alephParameters[0]=1.25202/50.;
 652         alephParameters[1]=2.74992e+01;
 653         alephParameters[2]=TMath::Exp(-3.31517e+01);
 654         alephParameters[3]=2.46246;
 655         alephParameters[4]=6.78938;
 656         periodsname[3]="dataLHC10hpass1";
 657       }
 658       if (set==1){//pass2 (AOD044)
 659         alephParameters[0]=1.25202/50.;
 660         alephParameters[1]=2.74992e+01;
 661         alephParameters[2]=TMath::Exp(-3.31517e+01);
 662         alephParameters[3]=2.46246;
 663         alephParameters[4]=6.78938;
 664         periodsname[4]="dataLHC10hpass2";
 665       }
 666     }
 667     if(iperiod==3){ //MC
 668       if(set==0){
 669         alephParameters[0] = 2.15898e+00/50.;
 670         alephParameters[1] = 1.75295e+01;
 671         alephParameters[2] = 3.40030e-09;
 672         alephParameters[3] = 1.96178e+00;
 673         alephParameters[4] = 3.91720e+00;
 674         periodsname[5]="MCold";
 675       }
 676       if(set==1){ //new
 677         alephParameters[0] = 1.44405/50;
 678         alephParameters[1] = 2.35409e+01;
 679         alephParameters[2] = TMath::Exp(-2.90330e+01);
 680         alephParameters[3] = 2.10681;
 681         alephParameters[4] = 4.62254;
 682         periodsname[6]="MCnew";
 683       }
 684
 685       if(set==2){ //new BB from Francesco
 686         alephParameters[0] = 0.029021;
 687         alephParameters[1] = 25.4181;
 688         alephParameters[2] = 4.66596e-08;
 689         alephParameters[3] = 1.90008;
 690         alephParameters[4] = 4.63783;
 691         periodsname[7]="MCBBFrancesco";
 692       }
 693
 694       if(set==3){ //low energy 2011
 695         alephParameters[0] = 0.0207667;
 696         alephParameters[1] = 29.9936;
 697         alephParameters[2] = 3.87866e-11;
 698         alephParameters[3] = 2.17291;
 699         alephParameters[4] = 7.1623;
 700         //periodsname[8]="MClowen2011";
 701       }
 702
 703
 704     }
 705     //cout<<periodsname[iperiod]<<endl;
 706     tpcResp->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
 707     cout<<"here"<<endl;
 708     for(Int_t ipart=0;ipart<npart;ipart++){
 709
 710       const Int_t n=1000;
 711       Double_t p[n],bethe[n];
 712
 713       for(Int_t k=0;k<n;k++){ //loop on the momentum steps
 714         p[k]=0.0001+k*4./n; //limits 0.-4. GeV/c
 715         //cout<<p[k]<<"\t";
 716         //bethe[k]=-tpcResp->Bethe(p[k]/masses[ipart]);
 717         AliPID::EParticleType ptype=AliPID::kKaon;
 718         if(ipart==1) ptype=AliPID::kPion;
 719         if(ipart==2) ptype=AliPID::kElectron;
 720         if(ipart==3) ptype=AliPID::kProton;
 721         bethe[k]=tpcResp->GetExpectedSignal(p[k],ptype);
 722       }
 723       //cout<<endl;
 724       TGraph *gr=new TGraph(n,p,bethe);
 725       gr->SetName(Form("%sP%d",partnames[ipart].Data(),iperiod));
 726       gr->SetTitle(Form("%sP%d;p (GeV/c);",partnames[ipart].Data(),iperiod));
 727       gr->SetLineColor(ipart+1);
 728       gr->SetMarkerColor(ipart+1);
 729       gr->GetYaxis()->SetRangeUser(35,100);
 730       cBethe->cd();
 731       if(iperiod==0 && ipart==0)gr->DrawClone("AL");
 732       else gr->DrawClone("L");
 733
 734       fout->cd();
 735       gr->Write();
 736     }
 737
 738   }
 739   TParameter<int> sett;
 740   sett.SetVal(set);
 741   fout->cd();
 742   sett.Write();
 743
 744   fout->Close();
 745 }
 746
 747 void DrawOutputCentrality(TString partname="D0",TString textleg="",TString path="./", Bool_t superimpose=kFALSE,TString suffixdir="",TString filename=/*"AnalysisResults.root"*/"PWG3histograms.root"){
 748   gStyle->SetCanvasColor(0);
 749   gStyle->SetTitleFillColor(0);
 750   gStyle->SetStatColor(0);
 751   gStyle->SetPalette(1);
 752
 753   TString listname="outputCentrCheck",partname2="",path2="",suffixdir2="",filename2="PWG3histograms.root";
 754
 755   // if(superimpose){
 756   //   cout<<"Enter the names:\n>";
 757   //   cin>>name1;
 758   //   cout<<">";
 759   //   cin>>name2;
 760   // }
 761   // TString listname="outputTrack",name1="",name2="";
 762   TString tmp="y";
 763
 764   if(superimpose){
 765     cout<<"##Second file\n";
 766     cout<<"Enter the name:\n";
 767     cout<<">";
 768     cin>>partname2;
 769     cout<<"Are they in the same output file? (y/n)"<<endl;
 770     cin>>tmp;
 771     if(tmp=="n"){
 772       cout<<"Path: \n";
 773       cout<<">";
 774       cin>>path2;
 775       cout<<"Dir name:\n";
 776       cout<<">";
 777       cin>>suffixdir2;
 778       cout<<"Filename: "<<endl;
 779       cout<<">";
 780       cin>>filename2;
 781     }
 782
 783   }
 784   // Int_t nhist=1;
 785   // TString *name=0x0;
 786   // if(superimpose){
 787   //   cout<<"Number of histogram to superimpose: ";
 788   //   cin>>nhist;
 789   //   name=new TString[nhist];
 790   //   for (Int_t j=0;j<nhist;j++){
 791   //     cout<<">";
 792   //     cin>>name[j];
 793   //   }
 794   // }
 795
 796   TList* list;
 797   TH1F * hstat;
 798
 799   TString dirname="PWG3_D2H_QA",dirname2=dirname;
 800   dirname+=suffixdir;
 801   dirname2+=suffixdir2;
 802   Bool_t isRead=ReadFile(list,hstat,listname,partname.Data(),path,filename,dirname);
 803   if(!isRead) return;
 804   if(!list || !hstat){
 805     cout<<":-( null pointers..."<<endl;
 806     return;
 807   }
 808
 809   TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 810   pvtxt->SetBorderSize(0);
 811   pvtxt->SetFillStyle(0);
 812   pvtxt->AddText(partname);
 813
 814   TList* llist;
 815   TH1F* hhstat;
 816   if(superimpose){
 817     isRead=ReadFile(llist,hhstat,listname,partname2.Data(),path2,filename2,dirname2);
 818     if(!isRead) return;
 819     if(!llist || !hhstat){
 820       cout<<":-( null pointers..."<<endl;
 821       return;
 822     }
 823     TText *redtext=pvtxt->AddText(partname2);
 824     redtext->SetTextColor(kRed);
 825
 826   }
 827
 828
 829   TCanvas* cst=new TCanvas("cst","Stat");
 830   cst->SetGridy();
 831   cst->cd();
 832   Int_t nevents=hstat->GetBinContent(1);
 833   hstat->Draw("htext0");
 834   cst->SaveAs(Form("%s%s.png",hstat->GetName(),textleg.Data()));
 835   cst->SaveAs(Form("%s%s.eps",hstat->GetName(),textleg.Data()));
 836   Int_t nevents080=1,nnevents080=1;
 837
 838   //TCanvas *spare=new TCanvas("sparecv","Spare");
 839
 840   for(Int_t i=0;i<list->GetEntries();i++){
 841
 842     TClass* objtype=list->At(i)->IsA();
 843     TString tpname=objtype->GetName();
 844
 845     if(tpname=="TH1F"){
 846
 847       TH1F* h=(TH1F*)list->At(i);
 848       TH1F* hh=0x0;
 849       if(superimpose){
 850         hh=(TH1F*)llist->At(i);
 851       }
 852       if(!h || (superimpose && !hh)){
 853         cout<<"Histogram "<<i<<" not found"<<endl;
 854         continue;
 855       }
 856       if(superimpose){
 857         hhstat->SetLineColor(kRed);
 858         hh->SetLineColor(kRed);
 859       }
 860
 861       TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
 862       TPaveText *pvtxt2=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 863       pvtxt2->SetBorderSize(0);
 864       pvtxt2->SetFillStyle(0);
 865
 866       c->cd();
 867       c->SetGrid();
 868       c->SetLogy();
 869       Int_t entries=h->Integral();
 870       pvtxt2->AddText(Form("%.1f %s of the events",(Double_t)entries/(Double_t)nevents*100,"%"));
 871       h->Draw();
 872       if(superimpose) {
 873         hh->Draw("sames");
 874         pvtxt->Draw();
 875       }
 876       pvtxt2->Draw();
 877       c->SaveAs(Form("%s%s.pdf",c->GetName(),textleg.Data()));
 878       c->SaveAs(Form("%s%s.eps",c->GetName(),textleg.Data()));
 879     }
 880     if(tpname=="TH2F"){
 881       TH2F* h=(TH2F*)list->At(i);
 882       if(!h){
 883         cout<<"Histogram "<<i<<" not found"<<endl;
 884         continue;
 885       }
 886       TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
 887       TPaveText *pvtxt3=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 888       pvtxt3->SetBorderSize(0);
 889       pvtxt3->SetFillStyle(0);
 890
 891       c->cd();
 892       c->SetGrid();
 893       Int_t entries=h->Integral();
 894       pvtxt3->AddText(Form("%.1f %s of the events",(Double_t)entries/(Double_t)nevents*100,"%"));
 895       h->Draw("colz");
 896       c->SetLogz();
 897       pvtxt3->Draw();
 898       c->SaveAs(Form("%s%s.pdf",c->GetName(),textleg.Data()));
 899       c->SaveAs(Form("%s%s.eps",c->GetName(),textleg.Data()));
 900     }
 901   }
 902
 903
 904   listname="countersCentrality";
 905
 906   isRead=ReadFile(list,hstat,listname,partname.Data(),path,filename,dirname);
 907   if(!isRead) return;
 908   if(!list || !hstat){
 909     cout<<":-( null pointers..."<<endl;
 910     return;
 911   }
 912
 913
 914   if(superimpose){
 915     isRead=ReadFile(llist,hhstat,listname,partname2.Data(),path2,filename2,dirname2);
 916     if(!isRead) return;
 917     if(!llist || !hhstat){
 918       cout<<":-( null pointers..."<<endl;
 919       return;
 920     }
 921     TText *redtext=pvtxt->AddText(partname2);
 922     redtext->SetTextColor(kRed);
 923
 924   }
 925
 926   TH1F* hallcntr=0x0;
 927   TH1F* hhallcntr=0x0;
 928   cout<<"normalizing to 0-80% as a check"<<endl;
 929   Int_t ncentr=10;//check this
 930   TH1F* h020=0x0;
 931   TH1F* h2080=0x0;
 932   TH1F* hh020=0x0;
 933   TH1F* hh2080=0x0;
 934
 935   TCanvas *cvnocnt=new TCanvas("cvnocnt","No Centrality estimation",800,400);
 936   cvnocnt->Divide(2,1);
 937   TCanvas *ccent=0x0;
 938
 939   for(Int_t i=0;i<list->GetEntries();i++){
 940     AliCounterCollection* coll=(AliCounterCollection*)list->At(i);
 941     AliCounterCollection* colle=0x0;
 942     if(superimpose) colle=(AliCounterCollection*)llist->At(i);
 943     coll->SortRubric("run");//sort by run number
 944
 945     h020=0x0;
 946     h2080=0x0;
 947     hh020=0x0;
 948     hh2080=0x0;
 949
 950     hallcntr=0x0;
 951     hhallcntr=0x0;
 952
 953     TH1F* hbad=(TH1F*)coll->Get("run",Form("centralityclass:-990_-980"));
 954     cvnocnt->cd(i+1);
 955     if(hbad) hbad->Draw();
 956
 957     ccent=new TCanvas(Form("ccent%s",coll->GetName()),Form("Centrality vs Run (%s)",coll->GetName()),1400,800);
 958     ccent->SetTicky();
 959     ccent->Divide(4,2);
 960
 961     TH1F* hh=0x0;
 962
 963     for(Int_t ic=0;ic<8/*ncentr*/;ic++){ //normalizing to 0-80% as a check
 964
 965       TH1F* h=(TH1F*)coll->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
 966       h->SetName(Form("h%d%d",i,ic));
 967       if(!hallcntr) {
 968         hallcntr=(TH1F*)h->Clone("hallcntr");
 969         hallcntr->Sumw2();
 970       } else {
 971         hallcntr->Add(h);
 972       }
 973
 974       nevents080+=h->Integral();
 975
 976       if(superimpose){
 977         hh=(TH1F*)colle->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
 978         hh->SetName(Form("hh%d%d",i,ic));
 979         if(!hhallcntr) {
 980           hhallcntr=(TH1F*)hh->Clone("hhallcntr");
 981           hhallcntr->Sumw2();
 982         }else hhallcntr->Add(hh);
 983
 984         nnevents080+=hh->Integral();
 985
 986       }
 987     }
 988
 989     for(Int_t ic=0;ic<ncentr;ic++){
 990
 991       TH1F* h=(TH1F*)coll->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
 992       h->SetName(Form("h%d%d",i,ic));
 993       h->Sumw2();
 994
 995       if(ic>=0 && ic<=1){ //0-20
 996         if(!h020) {
 997           h020=(TH1F*)h->Clone(Form("h020%s",coll->GetName()));
 998           h020->SetTitle(Form("Centrality 0-20 %s",coll->GetName()));
 999           if(superimpose){
1000             hh020=(TH1F*)hh->Clone(Form("hh020%s",coll->GetName()));
1001             hh020->SetTitle(Form("Centrality 0-20 %s",coll->GetName()));
1002           }
1003         }
1004         else {
1005           h020->Add(h);
1006           if(superimpose)hh020->Add(hh);
1007         }
1008       }
1009       if(ic>=2 && ic<=7){ //20-80
1010         if(!h2080) {
1011           h2080=(TH1F*)h->Clone(Form("h2080%s",coll->GetName()));
1012           h2080->SetTitle(Form("Centrality 20-80 %s",coll->GetName()));
1013           if(superimpose){
1014             hh2080=(TH1F*)hh->Clone(Form("hh2080%s",coll->GetName()));
1015             hh2080->SetTitle(Form("Centrality 20-80 %s",coll->GetName()));
1016           }
1017         }
1018         else {
1019           h2080->Add(h);
1020           if(superimpose)hh2080->Add(hh);
1021         }
1022
1023       }
1024
1025       h->Divide(hallcntr);
1026
1027       if(ic<8){
1028         ccent->cd(ic+1);
1029         h->GetYaxis()->SetLabelSize(0.05);
1030         h->GetYaxis()->SetTitleOffset(1.5);
1031         h->SetMinimum(0);
1032         //h->GetYaxis()->SetRangeUser(0.,0.15);
1033         h->DrawClone();
1034       }
1035       /*
1036         if(ic==0&&i==0){
1037         spare->cd();
1038         h->Draw();
1039         }
1040       */
1041       // ccent->cd(1);
1042       // h->SetLineColor(ic+1);
1043       // if(ic==0)h->DrawClone();
1044       // else h->DrawClone("sames");
1045     }
1046     h020->Divide(hallcntr);
1047     if(superimpose){
1048       hh020->Divide(hhallcntr);
1049       hh020->SetLineColor(2);
1050       hh020->SetMarkerColor(2);
1051     }
1052
1053     /*//draw 0-20 and 20-80 in the multi pad canvas (increase divisions before uncommenting)
1054     ccent->cd(ncentr+1);
1055     h020->DrawClone();
1056     if(superimpose){
1057       hh020->DrawClone("sames");
1058     }
1059     */
1060     TCanvas* cv020=new TCanvas(Form("cv020-%d",i),"0-20% vs run number",1400,600);
1061     cv020->cd();
1062     cv020->SetTicky();
1063     h020->GetYaxis()->SetRangeUser(0.,1.);
1064     h020->DrawClone();
1065     if(superimpose)hh020->DrawClone("sames");
1066     cv020->SaveAs(Form("cv020-%d.pdf",i));
1067     cv020->SaveAs(Form("cv020-%d.eps",i));
1068
1069     h2080->Divide(hallcntr);
1070     if(superimpose) {
1071       hh2080->Divide(hhallcntr);
1072       hh2080->SetLineColor(2);
1073       hh2080->SetMarkerColor(2);
1074     }
1075
1076     /*
1077     ccent->cd(ncentr+2);
1078     h2080->DrawClone();
1079
1080     if(superimpose){
1081       hh2080->DrawClone("sames");
1082     }
1083     */
1084     TCanvas* cv2080=new TCanvas(Form("cv2080-%d",i),"20-80% vs run number",1400,600);
1085     cv2080->cd();
1086     cv2080->SetTicky();
1087     h2080->GetYaxis()->SetRangeUser(0.,1.);
1088     h2080->DrawClone();
1089     if(superimpose)hh2080->DrawClone("sames");
1090     cv2080->SaveAs(Form("cv2080-%d.pdf",i));
1091     cv2080->SaveAs(Form("cv2080-%d.eps",i));
1092
1093     ccent->SaveAs(Form("%s%s.pdf",ccent->GetName(),textleg.Data()));
1094     ccent->SaveAs(Form("%s%s.eps",ccent->GetName(),textleg.Data()));
1095   }
1096
1097 }
1098
1099 void DrawProjections(TString partname="D0",TString h2dname="hMultvsPercentile",Int_t groupnbins=5,Float_t fitmin=15,Float_t fitmax=50,TString direction="X",TString path="./",TString suffixdir="", TString filename="AnalysisResults.root", TString fitfunc="pol0"/*option "nofit" does not fit*/){
1100   gStyle->SetCanvasColor(0);
1101   gStyle->SetTitleFillColor(0);
1102   gStyle->SetStatColor(0);
1103   gStyle->SetPalette(1);
1104
1105   TString listname="outputCentrCheck";
1106   TString dirname="PWG3_D2H_QA";
1107   dirname+=suffixdir;
1108
1109   TList* list;
1110   TH1F * hstat;
1111
1112   Bool_t isRead=ReadFile(list,hstat,listname,partname,path,filename,dirname);
1113   if(!isRead) return;
1114   if(!list || !hstat){
1115     cout<<":-( null pointers..."<<endl;
1116     return;
1117   }
1118   Double_t nevents=hstat->GetBinContent(5); //ev good vertex
1119
1120   TH2F* h2=(TH2F*)list->FindObject(h2dname);
1121   if(!h2){
1122     cout<<h2dname.Data()<<" not found"<<endl;
1123     list->ls();
1124     return;
1125   }
1126   TCanvas* cv2d=new TCanvas("cv2d",h2->GetName());
1127   cv2d->cd();
1128   cv2d->SetLogz();
1129   cv2d->SetGrid();
1130   h2->Draw("colz");
1131   TPaveText *pvst=new TPaveText(0.6,0.2,0.9,0.7,"NDC");
1132   pvst->SetBorderSize(0);
1133   pvst->SetFillStyle(0);
1134   pvst->AddText("Bin -> Cont/nEvVtx");
1135
1136   //nsteps=group bins in the Y(X) direction if projecting on the X(Y) direction
1137   Int_t nsteps=0;
1138
1139   if(direction=="X") nsteps=h2->GetNbinsY()/groupnbins;
1140   if(direction=="Y") nsteps=h2->GetNbinsX()/groupnbins;
1141   cout<<"Grouping bins by " <<groupnbins<<" I obtaine "<<nsteps<<" projections"<<endl;
1142
1143   TCanvas *cvpj=new TCanvas(Form("cvpj%s%s",direction.Data(),h2dname.Data()),Form("cvpj%s",direction.Data()),1200,800);
1144   cvpj->Divide((Int_t)(nsteps/3)+1,3);
1145   TFile* fout=new TFile(Form("proj%s%s.root",direction.Data(),h2dname.Data()), "recreate");
1146   //Float_t maxx[nsteps];
1147   //Float_t maxx[12]={9000,9000,6000,4000,2000,1400,800,500,200,100,40,25};
1148   Double_t integralpernev[nsteps];
1149
1150   Double_t minx=0,maxx=0;
1151   if(direction=="X"){
1152     minx=h2->GetYaxis()->GetXmin();
1153     maxx=h2->GetYaxis()->GetXmax();
1154   }
1155   if(direction=="Y"){
1156     minx=h2->GetXaxis()->GetXmin();
1157     maxx=h2->GetXaxis()->GetXmax();
1158   }
1159   printf("Plotting from %.1f to %.1f\n",minx,maxx);
1160   TCanvas *cintegral=new TCanvas("cintegral","Integral of each projection");
1161   TH1F* hint=new TH1F("hint","Integral of each projection;Centrality (%);Entries",nsteps,minx,maxx);
1162   Double_t minint=999999999,maxint=0;
1163
1164   for(Int_t i=0;i<nsteps;i++){
1165     TH1F* h=0x0;
1166     // if(direction=="X")h=(TH1F*)h2->ProjectionX(Form("px%d",i),i+kbins,i+2*kbins);
1167     // if(direction=="Y")h=(TH1F*)h2->ProjectionY(Form("py%d",i),i+kbins,i+2*kbins);
1168     if(direction=="X")h=(TH1F*)h2->ProjectionX(Form("px%d",i),groupnbins*i+1,groupnbins*(i+1));
1169     if(direction=="Y")h=(TH1F*)h2->ProjectionY(Form("py%d",i),groupnbins*i+1,groupnbins*(i+1));
1170     Double_t projint=h->Integral();
1171     cout<<"Integral of projection "<<i<<" = "<<projint<<endl;
1172     hint->SetBinContent(i+1,projint);
1173     hint->SetBinError(i+1,TMath::Sqrt(projint));
1174
1175     if(projint<1e-7) continue;
1176     if(minint>projint) minint=projint;
1177     if(projint>maxint) maxint=projint;
1178     integralpernev[i]=h->Integral()/nevents;
1179
1180     TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
1181     pvtxt->SetBorderSize(0);
1182     pvtxt->SetFillStyle(0);
1183     pvtxt->AddText(Form("%.0f - %.0f",h2->GetYaxis()->GetBinLowEdge(groupnbins*i+1),h2->GetYaxis()->GetBinLowEdge(groupnbins*(i+1))));
1184     pvst->AddText(Form("%.0f - %.0f -> %.2f",h2->GetYaxis()->GetBinLowEdge(groupnbins*i+1),h2->GetYaxis()->GetBinLowEdge((groupnbins*(i+1))),integralpernev[i]));
1185
1186     cvpj->cd(i+1);
1187     //h->GetXaxis()->SetRangeUser(0,maxx[i]);
1188     h->Draw();
1189     pvtxt->Draw();
1190     fout->cd();
1191     h->Write();
1192   }
1193   cvpj->SaveAs(Form("cvpj%s%s.pdf",direction.Data(),h2dname.Data()));
1194   cvpj->SaveAs(Form("cvpj%s%s.eps",direction.Data(),h2dname.Data()));
1195
1196   cv2d->cd();
1197   pvst->Draw();
1198   cv2d->SaveAs(Form("%s.pdf",h2->GetName()));
1199   cv2d->SaveAs(Form("%s.eps",h2->GetName()));
1200
1201   cintegral->cd();
1202   hint->SetMarkerStyle(20);
1203   hint->Draw("PE");
1204   if(!fitfunc.Contains("nofit")){
1205     hint->Fit(fitfunc.Data(),"RL","PE",fitmin,fitmax);
1206     TF1* fpolfit=hint->GetFunction(fitfunc.Data());
1207     TPaveText *txtvar=new TPaveText(0.3,0.1,0.9,0.4,"NDC");
1208     txtvar->SetBorderSize(0);
1209     txtvar->SetFillStyle(0);
1210     //txtvar->AddText(Form("Full spread %.0f- %.0f",maxint,minint));
1211     txtvar->AddText(Form("Fit in %.1f-%.1f; ",fitmin,fitmax));
1212     for(Int_t ipar=0;ipar<fpolfit->GetNpar();ipar++){
1213       txtvar->AddText(Form("par%d = %.0f, ",ipar, fpolfit->GetParameter(ipar)));
1214     }
1215     txtvar->AddText(Form("#tilde{#chi}^{2} = %.2f",fpolfit->GetChisquare()/fpolfit->GetNDF()));
1216     txtvar->AddText(Form("bin width = %.1f %s",hint->GetBinWidth(3),"%"));
1217     txtvar->Draw();
1218   }
1219   fout->cd();
1220   hint->Write();
1221   cintegral->SaveAs(Form("%s.pdf",hint->GetName()));
1222   cintegral->SaveAs(Form("%s.eps",hint->GetName()));
1223 }
1224
1225 void DrawEventSelection(TString partname="D0", TString path="./",TString suffixdir="",TString filename="AnalysisResults.root"){
1226   gStyle->SetCanvasColor(0);
1227   gStyle->SetTitleFillColor(0);
1228   gStyle->SetStatColor(0);
1229   gStyle->SetPalette(1);
1230   gStyle->SetOptStat(0);
1231
1232   TString listname="outputEvSel";
1233   TString dirname="PWG3_D2H_QA";
1234   dirname+=suffixdir;
1235
1236   TList* list;
1237   TH1F * hstat;
1238
1239   Bool_t isRead=ReadFile(list,hstat,listname,partname,path,filename,dirname);
1240   if(!isRead) return;
1241   if(!list || !hstat){
1242     cout<<":-( null pointers..."<<endl;
1243     return;
1244   }
1245   //Double_t neventsgv=hstat->Integral(5,5); //ev good vertex
1246
1247   for(Int_t i=0;i<list->GetEntries();i++){
1248
1249     TClass* objtype=list->At(i)->IsA();
1250     TString tpname=objtype->GetName();
1251
1252     if(tpname=="TH1F"){
1253       TH1F* htmp=(TH1F*)list->At(i);
1254       TCanvas* c=new TCanvas(Form("c%s",htmp->GetName()),Form("c%s",htmp->GetName()));
1255       c->cd();
1256       htmp->Draw();
1257       c->SaveAs(Form("%s.pdf",htmp->GetName()));
1258       c->SaveAs(Form("%s.eps",htmp->GetName()));
1259     }
1260
1261     if(tpname=="TH2F"){
1262       TH2F* htmp=(TH2F*)list->At(i);
1263       TCanvas* c=new TCanvas(Form("c%s",htmp->GetName()),Form("c%s",htmp->GetName()),1200,800);
1264       c->cd();
1265       htmp->SetMarkerSize(1.3);
1266       htmp->Draw("colzhtext45");
1267       c->SaveAs(Form("%s.pdf",htmp->GetName()));
1268       c->SaveAs(Form("%s.eps",htmp->GetName()));
1269     }
1270   }
1271
1272   AliCounterCollection* coll=(AliCounterCollection*)list->FindObject("trigCounter");
1273   if(!coll) {
1274     cout<<"Trigger counter not found"<<endl;
1275     return;
1276   }
1277
1278   coll->SortRubric("run");//sort by run number
1279
1280   TString collname=coll->GetName();
1281
1282   TString keywords=coll->GetKeyWords("triggertype");
1283
1284   Int_t nkeys=keywords.CountChar(',')+1;
1285
1286   TString *words = new TString[nkeys];
1287   for(Int_t k=0;k<nkeys;k++) words[k]="";
1288   printf("Keywords: ");
1289   Int_t count=0;
1290   for(Int_t l=0;l<keywords.Length();l++){
1291     if(keywords[l] != ',') words[count]+=keywords[l];
1292     else {
1293       printf("%s ",words[count].Data());
1294       count++;
1295     }
1296   }
1297   cout<<endl;
1298
1299   TH1D** htrig=new TH1D*[nkeys]; //each trigger type in one histogram of counts vs run
1300   TH1D** htrignorm=new TH1D*[nkeys]; //normalized to the counts in kAny
1301   TCanvas* ctrigfraction=new TCanvas("cvtrigfrac","Fraction of given trigger type vs run",1400,800);
1302   TLegend* legtr=new TLegend(0.15,0.5,0.35,0.8);
1303   legtr->SetBorderSize(0);
1304   legtr->SetFillStyle(0);
1305   for(Int_t k=0;k<nkeys;k++){
1306     htrig[k]=coll->Get("run",Form("triggertype:%s",words[k].Data()));
1307     htrig[k]->SetName(Form("h%s",words[k].Data()));
1308     htrig[k]->SetTitle("Trigger type;RUN; counts");
1309     htrig[k]->SetMarkerColor(k+1);
1310     htrig[k]->SetMarkerStyle(k+20);
1311     htrig[k]->Sumw2();
1312     legtr->AddEntry(htrig[k],Form("%s",words[k].Data()),"P");
1313     //drawings
1314     //1) counts of a given trigger over counts in kAny
1315     htrignorm[k]=(TH1D*)htrig[k]->Clone(Form("h%snormAny",words[k].Data()));
1316     htrignorm[k]->SetTitle("Trigger type over ANY trigger;RUN; counts/countsANY");
1317     htrignorm[k]->Divide(htrig[k],htrig[0],1.,1.,"B");
1318     htrignorm[k]->GetXaxis()->SetRangeUser(0,1.1);
1319
1320     ctrigfraction->cd();
1321     if(k>0)htrignorm[k]->Draw("PEsames");
1322     else htrignorm[k]->Draw("PE");
1323   }
1324
1325   ctrigfraction->cd();
1326   legtr->Draw();
1327   ctrigfraction->SaveAs("TrgFractionOverANY.pdf");
1328   ctrigfraction->SaveAs("TrgFractionOverANY.eps");
1329
1330   delete [] words;
1331
1332 }