PWG3/vertexingHF/macros/DrawQAoutput.C

   1 #include <Riostream.h>
   2 #include <TFile.h>
   3 #include <TString.h>
   4 #include <TH2F.h>
   5 #include <TH1F.h>
   6 #include <TF1.h>
   7 #include <TGraph.h>
   8 #include <TDirectoryFile.h>
   9 #include <TList.h>
  10 #include <TCanvas.h>
  11 #include <TLegend.h>
  12 #include <TPaveText.h>
  13 #include <TStyle.h>
  14 #include <TClass.h>
  15 #include <TDatabasePDG.h>
  16 #include <TParameter.h>
  17 #include <AliCounterCollection.h>
  18 #include <AliRDHFCuts.h>
  19
  20 TString *periodsname;
  21
  22 //read the file and take list and stat
  23
  24 Bool_t ReadFile(TList* &list,TH1F* &hstat, TString listname,TString partname,TString path="./",TString filename=/*"AnalysisResults.root"*/"PWG3histograms.root");
  25 Bool_t ReadFileMore(TList* &list,TH1F* &hstat, AliRDHFCuts* &cutobj, TString listname,TString partname,TString path="./",TString filename=/*"AnalysisResults.root"*/"PWG3histograms.root");
  26 void SuperimposeBBToTPCSignal(Int_t period /*0=LHC10bc, 1=LHC10d, 2=LHC10h*/,TCanvas* cpid, Int_t set);
  27 void TPCBetheBloch(Int_t set);
  28
  29 Bool_t ReadFile(TList* &list,TH1F* &hstat, TString listname,TString partname,TString path,TString filename){
  30
  31   TString hstatname="nEntriesQA",dirname="PWG3_D2H_QA", cutobjname="";
  32   filename.Prepend(path);
  33   listname+=partname;
  34   hstatname+=partname;
  35
  36   TFile* f=new TFile(filename.Data());
  37   if(!f){
  38     cout<<filename.Data()<<" not found"<<endl;
  39     return kFALSE;
  40   }
  41   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname);
  42   if(!f){
  43     cout<<dirname.Data()<<" not found  in "<<filename.Data()<<endl;
  44     return kFALSE;
  45   }
  46
  47   list=(TList*)dir->Get(listname);
  48   if(!list){
  49     cout<<"List "<<listname.Data()<<" not found"<<endl;
  50     dir->ls();
  51     return kFALSE;
  52   }
  53
  54   hstat=(TH1F*)dir->Get(hstatname);
  55   if(!hstat){
  56     cout<<hstatname.Data()<<" not found"<<endl;
  57     return kFALSE;
  58   }
  59
  60   return kTRUE;
  61 }
  62
  63 Bool_t ReadFileMore(TList* &list,TH1F* &hstat, AliRDHFCuts* &cutobj, TString listname,TString partname,TString path,TString filename){
  64
  65   TString hstatname="nEntriesQA",dirname="PWG3_D2H_QA", cutobjname="";
  66   filename.Prepend(path);
  67   listname+=partname;
  68   hstatname+=partname;
  69
  70   if(partname.Contains("Dplus")) cutobjname="AnalysisCuts";//"DplustoKpipiCutsStandard";
  71   else{
  72     if(partname.Contains("D0")) cutobjname="D0toKpiCuts";//"D0toKpiCutsStandard"
  73     else{
  74       if(partname.Contains("Dstar")) cutobjname="DStartoKpipiCuts";
  75       else{
  76         if(partname.Contains("Ds")) cutobjname="DstoKKpiCuts";
  77         else{
  78           if(partname.Contains("D04")) cutobjname="D0toKpipipiCuts";
  79           else{
  80             if(partname.Contains("Lc")) cutobjname="LctopKpiAnalysisCuts";
  81           }
  82         }
  83       }
  84     }
  85   }
  86
  87   TFile* f=new TFile(filename.Data());
  88   if(!f){
  89     cout<<filename.Data()<<" not found"<<endl;
  90     return kFALSE;
  91   }
  92   TDirectoryFile* dir=(TDirectoryFile*)f->Get(dirname);
  93   if(!f){
  94     cout<<dirname.Data()<<" not found  in "<<filename.Data()<<endl;
  95     return kFALSE;
  96   }
  97
  98   list=(TList*)dir->Get(listname);
  99   if(!list){
 100     cout<<"List "<<listname.Data()<<" not found"<<endl;
 101     dir->ls();
 102     return kFALSE;
 103   }
 104
 105   hstat=(TH1F*)dir->Get(hstatname);
 106   if(!hstat){
 107     cout<<hstatname.Data()<<" not found"<<endl;
 108     return kFALSE;
 109   }
 110
 111   cutobj=(AliRDHFCuts*)dir->Get(cutobjname);
 112   if(!cutobj){
 113     cout<<cutobjname.Data()<<" not found"<<endl;
 114     return kFALSE;
 115   }
 116
 117   return kTRUE;
 118 }
 119
 120 //draw "track related" histograms (list "outputTrack")
 121 void DrawOutputTrack(TString partname="D0",TString textleg="",TString path="./", Bool_t superimpose=kFALSE){
 122   gStyle->SetCanvasColor(0);
 123   gStyle->SetTitleFillColor(0);
 124   gStyle->SetStatColor(0);
 125   gStyle->SetPalette(1);
 126
 127   TString listname="outputTrack",name1="",name2="",path2="",filename=/*"AnalysisResults.root"*/"PWG3histograms.root",filename2="PWG3histograms.root";
 128   TString tmp="y";
 129
 130   if(superimpose){
 131     cout<<"Enter the names:\n>";
 132     cin>>name1;
 133     cout<<">";
 134     cin>>name2;
 135     cout<<"Are they in the same output file? (y/n)"<<endl;
 136     cin>>tmp;
 137     if(tmp=="n"){
 138       cout<<"Path: \n";
 139       cout<<">";
 140       cin>>path2;
 141       cout<<"Filename: "<<endl;
 142       cout<<">";
 143       cin>>filename2;
 144     }
 145
 146   }
 147
 148   TList* list;
 149   TH1F * hstat;
 150
 151   Bool_t isRead=ReadFile(list,hstat,listname,Form("%s%s",partname.Data(),name1.Data()),path,filename);
 152   if(!isRead) return;
 153   if(!list || !hstat){
 154     cout<<":-( null pointers..."<<endl;
 155     return;
 156   }
 157   TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 158   pvtxt->SetBorderSize(0);
 159   pvtxt->SetFillStyle(0);
 160   pvtxt->AddText(name1);
 161
 162   TList* llist;
 163   TH1F* hhstat;
 164   if(superimpose){
 165     isRead=ReadFile(llist,hhstat,listname,Form("%s%s",partname.Data(),name2.Data()),path2,filename2);
 166     if(!isRead) return;
 167     if(!llist || !hhstat){
 168       cout<<":-( null pointers..."<<endl;
 169       return;
 170     }
 171     TText *redtext=pvtxt->AddText(name2);
 172     redtext->SetTextColor(kRed);
 173     hhstat->Scale(hstat->Integral()/hhstat->Integral());
 174
 175   }
 176
 177   for(Int_t i=0;i<list->GetEntries();i++){
 178     TH1F* h=(TH1F*)list->At(i);
 179     TH1F* hh=0x0;
 180     if(superimpose){
 181       hh=(TH1F*)llist->At(i);
 182       hh->Scale(h->Integral()/hh->Integral());
 183     }
 184     if(!h || (superimpose && !hh)){
 185       cout<<"Histogram "<<i<<" not found"<<endl;
 186       continue;
 187     }
 188     if(superimpose){
 189       hhstat->SetLineColor(kRed);
 190       hh->SetLineColor(kRed);
 191     }
 192
 193     TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
 194     c->cd();
 195     c->SetGrid();
 196     TString hname=h->GetName();
 197     if(!hname.Contains("nCls")){
 198       c->SetLogy();
 199       //h->SetMinimum(1);
 200       h->Draw();
 201       if(superimpose) hh->Draw("sames");
 202     } else {
 203       h->Draw("htext0");
 204       if(superimpose)hh->Draw("htext0sames");
 205     }
 206     pvtxt->Draw();
 207     c->SaveAs(Form("%s%s%s%s.png",c->GetName(),name1.Data(),name2.Data(),textleg.Data()));
 208   }
 209
 210   TCanvas* cst=new TCanvas("cst","Stat");
 211   cst->SetGridy();
 212   cst->cd();
 213   hstat->Draw("htext0");
 214   if(superimpose) {
 215     hhstat->Draw("htext0sames");
 216     pvtxt->Draw();
 217   }
 218   cst->SaveAs(Form("%s%s.png",hstat->GetName(),textleg.Data()));
 219
 220
 221 }
 222
 223 //draw "pid related" histograms (list "outputPID")
 224 //period=-999 to draw the pull instead of the cut
 225 void DrawOutputPID(TString partname="D0", Int_t mode=0/*0=with pull, 1=with nsigma*/,TString textleg="",TString path="./"){
 226   gStyle->SetCanvasColor(0);
 227   gStyle->SetTitleFillColor(0);
 228   gStyle->SetStatColor(0);
 229   gStyle->SetPalette(1);
 230
 231   Int_t period=2 ,set=0;
 232   if(mode==1){
 233     cout<<"Choose period: \n-LHC10h -> 2;\n-LHC10de -> 1;\n-LHC10bc -> 0"<<endl;
 234     cin>>period;
 235     if(period>0){
 236       cout<<"Choose set: "<<endl;
 237       if(period==2) cout<<"-pass1 -> 0;\n-pass2 -> 1"<<endl;
 238       cin>>set;
 239     }
 240   }
 241
 242   TString listname="outputPid";
 243
 244   TList* list;
 245   TH1F * hstat;
 246   //needed only for mode 1
 247   AliRDHFCuts* cutobj;
 248   AliAODPidHF* aodpid;
 249   Double_t nsigmaTOF=0;
 250   Double_t nsigmaTPC[3]={},plimTPC[2]={};
 251
 252   if(mode==1){
 253     Bool_t isRead=ReadFileMore(list,hstat,cutobj,listname,partname,path);
 254     if(!isRead) return;
 255     if(!list || !hstat){
 256       cout<<":-( null pointers..."<<endl;
 257       return;
 258     }
 259     aodpid=(AliAODPidHF*)cutobj->GetPidHF();
 260     if(!aodpid){
 261       cout<<"PidHF object not found! cannot get the nsigma values"<<endl;
 262       return;
 263     }
 264     nsigmaTOF=aodpid->GetSigma(3);
 265
 266     nsigmaTPC[0]=aodpid->GetSigma(0);
 267     nsigmaTPC[1]=aodpid->GetSigma(1);
 268     nsigmaTPC[2]=aodpid->GetSigma(2);
 269     aodpid->GetPLimit(plimTPC);
 270
 271   }else{
 272
 273     Bool_t isRead=ReadFile(list,hstat,listname,partname,path);
 274     if(!isRead) return;
 275     if(!list || !hstat){
 276       cout<<":-( null pointers..."<<endl;
 277       return;
 278     }
 279   }
 280
 281
 282   TPaveText *txtsigmaTOF=new TPaveText(0.1,0.65,0.5,0.9,"NDC");
 283   txtsigmaTOF->SetBorderSize(0);
 284   txtsigmaTOF->SetFillStyle(0);
 285   txtsigmaTOF->AddText(Form("nsigmacut from cutobj = %.1f",nsigmaTOF));
 286   TLine lTOF;
 287   lTOF.SetLineColor(kMagenta+1);
 288   lTOF.SetLineStyle(2);
 289   lTOF.SetLineWidth(3);
 290
 291   TPaveText *txtsigmaTPC=new TPaveText(0.3,0.6,0.6,0.9,"NDC");
 292   txtsigmaTPC->SetBorderSize(0);
 293   txtsigmaTPC->SetFillStyle(0);
 294   txtsigmaTPC->AddText("nsigmacut from cutobj \n");
 295   txtsigmaTPC->AddText(Form("p < %.1f : %.1f \n",plimTPC[0],nsigmaTPC[0]));
 296   txtsigmaTPC->AddText(Form("%.1f < p < %.1f : %.1f \n",plimTPC[0],plimTPC[1],nsigmaTPC[1]));
 297   txtsigmaTPC->AddText(Form("p > %.1f : %.1f \n",plimTPC[1],nsigmaTPC[2]));
 298   TLine lTPC;
 299   lTPC.SetLineColor(kMagenta+1);
 300   lTPC.SetLineStyle(2);
 301   lTPC.SetLineWidth(3);
 302
 303   // TCanvas *ctest=new TCanvas("text","Test text");
 304   // ctest->cd();
 305   // txtsigmaTPC->Draw();
 306   // txtsigmaTOF->Draw();
 307
 308
 309   for(Int_t i=0;i<list->GetEntries();i++){
 310     TClass* objtype=list->At(i)->IsA();
 311     TString tpname=objtype->GetName();
 312
 313     if(tpname=="TH1F"){
 314       TH1F* h=(TH1F*)list->At(i);
 315
 316       if(!h){
 317         cout<<"Histogram "<<i<<" not found"<<endl;
 318         continue;
 319       }
 320       //h->Scale(1./h->Integral("width"));
 321       TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
 322       c->SetLogz();
 323       c->cd();
 324       h->Draw();
 325
 326       //write
 327       c->SaveAs(Form("%s%s.png",h->GetName(),textleg.Data()));
 328       TFile* fout=new TFile(Form("%s.root",h->GetName()),"recreate");
 329       fout->cd();
 330       c->Write();
 331     }
 332
 333     if(tpname=="TH2F"){
 334       TH2F* h=(TH2F*)list->At(i);
 335
 336       if(!h){
 337         cout<<"Histogram "<<i<<" not found"<<endl;
 338         continue;
 339       }
 340       h->Sumw2();
 341       h->Scale(1./h->Integral("width"));
 342       TString hname=h->GetName();
 343
 344       if(hname.Contains("hTOFtimeKaonHyptime")){
 345         TCanvas* cz=new TCanvas(Form("c%szoom",hname.Data()),Form("%szoom",hname.Data()));
 346         cz->SetLogz();
 347         TH2F* hz=(TH2F*)h->Clone(Form("%sz",hname.Data()));
 348         hz->Draw("colz");
 349         hz->SetAxisRange(-1500,1500,"Y");
 350         //write
 351         cz->SaveAs(Form("%szoom.png",h->GetName()));
 352       }
 353
 354       TCanvas* c=new TCanvas(Form("c%s",hname.Data()),hname.Data());
 355       c->SetLogz();
 356       //c->SetLogx();
 357       c->cd();
 358
 359       h->Draw("colz");
 360
 361       //TCanvas *test=new TCanvas("test","test");
 362       if(mode==0){
 363         //mean and pull, code from Jens Wiechula
 364         TF1 fg("fg","gaus",-2.,2.); // fit range +- 2 sigma
 365         TLine l;
 366         TObjArray arr;
 367
 368         //h->Draw("colz");
 369         fg.SetParameters(1,0,1);
 370         h->FitSlicesY(&fg,0,-1,0,"NQR",&arr);
 371
 372         TH1 *hM=(TH1*)arr.At(1);
 373         hM->SetMarkerStyle(20);
 374         hM->SetMarkerSize(.5);
 375         hM->DrawClone("sames");
 376
 377         TH1 *hS=(TH1*)arr.At(2);
 378         hS->SetMarkerStyle(20);
 379         hS->SetMarkerSize(.5);
 380         hS->SetMarkerColor(kRed);
 381         hS->SetLineColor(kRed);
 382         hS->DrawClone("same");
 383
 384         l.SetLineColor(kBlack);
 385         l.DrawLine(.2,0,20,0);
 386         l.SetLineColor(kRed);
 387         l.DrawLine(.2,1,20,1);
 388
 389       }else{ //mode 1
 390
 391         if(hname.Contains("TOFsigma")) {
 392
 393           c->cd();
 394           txtsigmaTOF->Draw();
 395           lTOF.DrawLine(.2,nsigmaTOF,20,nsigmaTOF);
 396           lTOF.DrawLine(.2,-1*nsigmaTOF,4.,-1*nsigmaTOF);
 397
 398         }
 399
 400
 401         if(hname.Contains("TPCsigma")){
 402
 403           c->cd();
 404           txtsigmaTPC->Draw();
 405           lTPC.DrawLine(0.,nsigmaTPC[0],plimTPC[0],nsigmaTPC[0]);
 406           lTPC.DrawLine(plimTPC[0],nsigmaTPC[1],plimTPC[1],nsigmaTPC[1]);
 407           lTPC.DrawLine(plimTPC[1],nsigmaTPC[2],4,nsigmaTPC[2]);
 408           lTPC.DrawLine(0.,-1*nsigmaTPC[0],plimTPC[0],-1*nsigmaTPC[0]);
 409           lTPC.DrawLine(plimTPC[0],-1*nsigmaTPC[1],plimTPC[1],-1*nsigmaTPC[1]);
 410           lTPC.DrawLine(plimTPC[1],-1*nsigmaTPC[2],4,-1*nsigmaTPC[2]);
 411         }
 412
 413         if(hname.Contains("TPCsigvsp")){
 414           SuperimposeBBToTPCSignal(period,c,set);
 415         }
 416       }
 417
 418       //write
 419       c->SaveAs(Form("%s%d.png",h->GetName(),mode));
 420       TFile* fout=new TFile(Form("%s%d.root",h->GetName(),mode),"recreate");
 421       fout->cd();
 422       c->Write();
 423     }
 424   }
 425 }
 426
 427 void SuperimposeBBToTPCSignal(Int_t period /*0=LHC10bc, 1=LHC10d, 2=LHC10h*/,TCanvas* cpid,Int_t set /*see below*/){
 428
 429   TFile* fBethe=new TFile("BetheBlochTPC.root");
 430   if(!fBethe->IsOpen()){
 431     TPCBetheBloch(set);
 432     fBethe=new TFile("BetheBlochTPC.root");
 433   }
 434   const Int_t npart=4;
 435   TString partnames[npart]={"Kaon","Pion","Electron","Proton"};
 436   for(Int_t ipart=0;ipart<npart;ipart++){
 437     TString grname=Form("%sP%d",partnames[ipart].Data(),period);
 438     TGraph* gr=(TGraph*)fBethe->Get(grname);
 439     cpid->cd();
 440     gr->SetLineColor(1);
 441     gr->SetLineWidth(2);
 442     gr->Draw("L");
 443   }
 444
 445   //cpid->SaveAs(Form("%sBB.png",hname.Data()));
 446 }
 447
 448 //draw and save Bethe Bloch from TPC in different periods
 449 void TPCBetheBloch(Int_t set){
 450   gStyle->SetOptTitle(0);
 451   gStyle->SetCanvasColor(0);
 452
 453   AliTPCPIDResponse *tpcResp=new AliTPCPIDResponse();
 454
 455   const Int_t npart=4;
 456   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*/};
 457   TString partnames[npart]={"Kaon","Pion","Electron","Proton"};
 458   //printf("%s = %.4f,%s = %.4f,%s = %.4f\n",partnames[0].Data(),masses[0],partnames[1].Data(),masses[1],partnames[2].Data(),masses[2]);
 459   TCanvas *cBethe=new TCanvas("cBethe","Bethe Bloch K pi e p");
 460   Int_t nperiods=4; //LHC10b+c, LHC10d, LHC10h, MC
 461   Double_t alephParameters[5]={};
 462   Int_t nsets=1/*LHC10bc*/+2/*LHC10de*/+2/*LHC10h*/+3/*MC*/;
 463
 464   periodsname=new TString[nsets];
 465   cout<<"Creating the file of the Bethe Bloch"<<endl;
 466   TFile* fout=new TFile("BetheBlochTPC.root","recreate");
 467
 468   for(Int_t iperiod=0;iperiod<nperiods;iperiod++){
 469     cout<<"Period "<<iperiod<<" : ";
 470     if(iperiod==0){ //LHC10bc
 471
 472       alephParameters[0] = 0.0283086/0.97;
 473       alephParameters[1] = 2.63394e+01;
 474       alephParameters[2] = 5.04114e-11;
 475       alephParameters[3] = 2.12543e+00;
 476       alephParameters[4] = 4.88663e+00;
 477       periodsname[0]="dataLHC10bc";
 478     }
 479     if(iperiod==1){ //LHC10de,low energy
 480       if(set==0){
 481         alephParameters[0] = 1.63246/50.;
 482         alephParameters[1] = 2.20028e+01;
 483         alephParameters[2] = TMath::Exp(-2.48879e+01);
 484         alephParameters[3] = 2.39804e+00;
 485         alephParameters[4] = 5.12090e+00;
 486         periodsname[1]="dataLHC10deold";
 487       }
 488       if(set==1){
 489         alephParameters[0] = 1.34490e+00/50.;
 490         alephParameters[1] =  2.69455e+01;
 491         alephParameters[2] =  TMath::Exp(-2.97552e+01);
 492         alephParameters[3] = 2.35339e+00;
 493         alephParameters[4] = 5.98079e+00;
 494         periodsname[2]="dataLHC10denew";
 495       }
 496     }
 497
 498     if(iperiod==2){//LHC10h
 499       if(set==0){//pass1
 500         alephParameters[0]=1.25202/50.;
 501         alephParameters[1]=2.74992e+01;
 502         alephParameters[2]=TMath::Exp(-3.31517e+01);
 503         alephParameters[3]=2.46246;
 504         alephParameters[4]=6.78938;
 505         periodsname[3]="dataLHC10hpass1";
 506       }
 507       if (set==1){//pass2 (AOD044)
 508         alephParameters[0]=1.25202/50.;
 509         alephParameters[1]=2.74992e+01;
 510         alephParameters[2]=TMath::Exp(-3.31517e+01);
 511         alephParameters[3]=2.46246;
 512         alephParameters[4]=6.78938;
 513         periodsname[4]="dataLHC10hpass2";
 514       }
 515     }
 516     if(iperiod==3){ //MC
 517       if(set==0){
 518         alephParameters[0] = 2.15898e+00/50.;
 519         alephParameters[1] = 1.75295e+01;
 520         alephParameters[2] = 3.40030e-09;
 521         alephParameters[3] = 1.96178e+00;
 522         alephParameters[4] = 3.91720e+00;
 523         periodsname[5]="MCold";
 524       }
 525       if(set==1){ //new
 526         alephParameters[0] = 1.44405/50;
 527         alephParameters[1] = 2.35409e+01;
 528         alephParameters[2] = TMath::Exp(-2.90330e+01);
 529         alephParameters[3] = 2.10681;
 530         alephParameters[4] = 4.62254;
 531         periodsname[6]="MCnew";
 532       }
 533
 534       if(set==2){ //new BB from Francesco
 535         alephParameters[0] = 0.029021;
 536         alephParameters[1] = 25.4181;
 537         alephParameters[2] = 4.66596e-08;
 538         alephParameters[3] = 1.90008;
 539         alephParameters[4] = 4.63783;
 540         periodsname[7]="MCBBFrancesco";
 541       }
 542
 543       if(set==3){ //low energy 2011
 544         alephParameters[0] = 0.0207667;
 545         alephParameters[1] = 29.9936;
 546         alephParameters[2] = 3.87866e-11;
 547         alephParameters[3] = 2.17291;
 548         alephParameters[4] = 7.1623;
 549         //periodsname[8]="MClowen2011";
 550       }
 551
 552
 553     }
 554     //cout<<periodsname[iperiod]<<endl;
 555     tpcResp->SetBetheBlochParameters(alephParameters[0],alephParameters[1],alephParameters[2],alephParameters[3],alephParameters[4]);
 556     cout<<"here"<<endl;
 557     for(Int_t ipart=0;ipart<npart;ipart++){
 558
 559       const Int_t n=1000;
 560       Double_t p[n],bethe[n];
 561
 562       for(Int_t k=0;k<n;k++){ //loop on the momentum steps
 563         p[k]=0.0001+k*4./n; //limits 0.-4. GeV/c
 564         //cout<<p[k]<<"\t";
 565         //bethe[k]=-tpcResp->Bethe(p[k]/masses[ipart]);
 566         AliPID::EParticleType ptype=AliPID::kKaon;
 567         if(ipart==1) ptype=AliPID::kPion;
 568         if(ipart==2) ptype=AliPID::kElectron;
 569         if(ipart==3) ptype=AliPID::kProton;
 570         bethe[k]=tpcResp->GetExpectedSignal(p[k],ptype);
 571       }
 572       //cout<<endl;
 573       TGraph *gr=new TGraph(n,p,bethe);
 574       gr->SetName(Form("%sP%d",partnames[ipart].Data(),iperiod));
 575       gr->SetTitle(Form("%sP%d;p (GeV/c);",partnames[ipart].Data(),iperiod));
 576       gr->SetLineColor(ipart+1);
 577       gr->SetMarkerColor(ipart+1);
 578       gr->GetYaxis()->SetRangeUser(35,100);
 579       cBethe->cd();
 580       if(iperiod==0 && ipart==0)gr->DrawClone("AL");
 581       else gr->DrawClone("L");
 582
 583       fout->cd();
 584       gr->Write();
 585     }
 586
 587   }
 588   TParameter<int> sett;
 589   sett.SetVal(set);
 590   fout->cd();
 591   sett.Write();
 592
 593   fout->Close();
 594 }
 595
 596 void DrawOutputCentrality(TString partname="D0",TString textleg="",TString path="./", Bool_t superimpose=kFALSE){
 597   gStyle->SetCanvasColor(0);
 598   gStyle->SetTitleFillColor(0);
 599   gStyle->SetStatColor(0);
 600   gStyle->SetPalette(1);
 601
 602   TString listname="outputCentrCheck",name1="",name2="",path2="",filename=/*"AnalysisResults.root"*/"PWG3histograms.root",filename2="PWG3histograms.root";
 603
 604   // if(superimpose){
 605   //   cout<<"Enter the names:\n>";
 606   //   cin>>name1;
 607   //   cout<<">";
 608   //   cin>>name2;
 609   // }
 610   // TString listname="outputTrack",name1="",name2="";
 611   TString tmp="y";
 612
 613   if(superimpose){
 614     cout<<"Enter the names:\n>";
 615     cin>>name1;
 616     cout<<">";
 617     cin>>name2;
 618     cout<<"Are they in the same output file? (y/n)"<<endl;
 619     cin>>tmp;
 620     if(tmp=="n"){
 621       cout<<"Path: \n";
 622       cout<<">";
 623       cin>>path2;
 624       cout<<"Filename: "<<endl;
 625       cout<<">";
 626       cin>>filename2;
 627     }
 628
 629   }
 630   // Int_t nhist=1;
 631   // TString *name=0x0;
 632   // if(superimpose){
 633   //   cout<<"Number of histogram to superimpose: ";
 634   //   cin>>nhist;
 635   //   name=new TString[nhist];
 636   //   for (Int_t j=0;j<nhist;j++){
 637   //     cout<<">";
 638   //     cin>>name[j];
 639   //   }
 640   // }
 641
 642   TList* list;
 643   TH1F * hstat;
 644
 645   Bool_t isRead=ReadFile(list,hstat,listname,Form("%s%s",partname.Data(),name1.Data()),path);
 646   if(!isRead) return;
 647   if(!list || !hstat){
 648     cout<<":-( null pointers..."<<endl;
 649     return;
 650   }
 651
 652   TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 653   pvtxt->SetBorderSize(0);
 654   pvtxt->SetFillStyle(0);
 655   pvtxt->AddText(name1);
 656
 657   TList* llist;
 658   TH1F* hhstat;
 659   if(superimpose){
 660     isRead=ReadFile(llist,hhstat,listname,Form("%s%s",partname.Data(),name2.Data()),path2);
 661     if(!isRead) return;
 662     if(!llist || !hhstat){
 663       cout<<":-( null pointers..."<<endl;
 664       return;
 665     }
 666     TText *redtext=pvtxt->AddText(name2);
 667     redtext->SetTextColor(kRed);
 668
 669   }
 670
 671
 672   TCanvas* cst=new TCanvas("cst","Stat");
 673   cst->SetGridy();
 674   cst->cd();
 675   Int_t nevents=hstat->Integral(1,2);
 676   hstat->Draw("htext0");
 677   cst->SaveAs(Form("%s%s.png",hstat->GetName(),textleg.Data()));
 678   Int_t nevents080=1,nnevents080=1;
 679
 680   //TCanvas *spare=new TCanvas("sparecv","Spare");
 681
 682   for(Int_t i=0;i<list->GetEntries();i++){
 683
 684     TClass* objtype=list->At(i)->IsA();
 685     TString tpname=objtype->GetName();
 686
 687     if(tpname=="TH1F"){
 688
 689       TH1F* h=(TH1F*)list->At(i);
 690       TH1F* hh=0x0;
 691       if(superimpose){
 692         hh=(TH1F*)llist->At(i);
 693       }
 694       if(!h || (superimpose && !hh)){
 695         cout<<"Histogram "<<i<<" not found"<<endl;
 696         continue;
 697       }
 698       if(superimpose){
 699         hhstat->SetLineColor(kRed);
 700         hh->SetLineColor(kRed);
 701       }
 702
 703       TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
 704       TPaveText *pvtxt2=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 705       pvtxt2->SetBorderSize(0);
 706       pvtxt2->SetFillStyle(0);
 707
 708       c->cd();
 709       c->SetGrid();
 710       c->SetLogy();
 711       Int_t entries=h->Integral();
 712       pvtxt2->AddText(Form("%.1f %s of the events",(Double_t)entries/(Double_t)nevents*100,"%"));
 713       h->Draw();
 714       if(superimpose) {
 715         hh->Draw("sames");
 716         pvtxt->Draw();
 717       }
 718       pvtxt2->Draw();
 719       c->SaveAs(Form("%s%s.png",c->GetName(),textleg.Data()));
 720     }
 721     if(tpname=="TH2F"){
 722       TH2F* h=(TH2F*)list->At(i);
 723       if(!h){
 724         cout<<"Histogram "<<i<<" not found"<<endl;
 725         continue;
 726       }
 727       TCanvas* c=new TCanvas(Form("c%s",h->GetName()),h->GetName());
 728       TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 729       pvtxt->SetBorderSize(0);
 730       pvtxt->SetFillStyle(0);
 731
 732       c->cd();
 733       c->SetGrid();
 734       Int_t entries=h->Integral();
 735       pvtxt->AddText(Form("%.1f %s of the events",(Double_t)entries/(Double_t)nevents*100,"%"));
 736       h->Draw("colz");
 737       c->SetLogz();
 738       pvtxt->Draw();
 739       c->SaveAs(Form("%s%s.png",c->GetName(),textleg.Data()));
 740     }
 741   }
 742
 743
 744   listname="countersCentrality";
 745
 746   isRead=ReadFile(list,hstat,listname,Form("%s%s",partname.Data(),name1.Data()),path);
 747   if(!isRead) return;
 748   if(!list || !hstat){
 749     cout<<":-( null pointers..."<<endl;
 750     return;
 751   }
 752
 753
 754   if(superimpose){
 755     isRead=ReadFile(llist,hhstat,listname,Form("%s%s",partname.Data(),name2.Data()),path2);
 756     if(!isRead) return;
 757     if(!llist || !hhstat){
 758       cout<<":-( null pointers..."<<endl;
 759       return;
 760     }
 761     TText *redtext=pvtxt->AddText(name2);
 762     redtext->SetTextColor(kRed);
 763
 764   }
 765
 766   TH1F* hallcntr=0x0;
 767   TH1F* hhallcntr=0x0;
 768   cout<<"normalizing to 0-80% as a check"<<endl;
 769
 770   TCanvas *cvnocnt=new TCanvas("cvnocnt","No Centrality estimation",800,400);
 771   cvnocnt->Divide(2,1);
 772
 773   for(Int_t i=0;i<list->GetEntries();i++){
 774     AliCounterCollection* coll=(AliCounterCollection*)list->At(i);
 775     AliCounterCollection* colle=0x0;
 776     if(superimpose) colle=(AliCounterCollection*)llist->At(i);
 777     coll->SortRubric("run");//sort by run number
 778     //coll->PrintKeyWords();
 779     Int_t ncentr=10;//check this
 780     TH1F* h020=0x0;
 781     TH1F* h2080=0x0;
 782     TH1F* hh020=0x0;
 783     TH1F* hh2080=0x0;
 784     hallcntr=0x0;
 785     hhallcntr=0x0;
 786
 787     TH1F* hbad=(TH1F*)coll->Get("run",Form("centralityclass:-990_-980"));
 788     cvnocnt->cd(i+1);
 789     if(hbad) hbad->Draw();
 790
 791     TCanvas *ccent=new TCanvas(Form("ccent%s",coll->GetName()),Form("Centrality vs Run (%s)",coll->GetName()),1400,800);
 792     ccent->Divide(5,3);
 793
 794     TH1F* hh=0x0;
 795
 796     for(Int_t ic=0;ic<8/*ncentr*/;ic++){ //normalizing to 0-80% as a check
 797
 798       TH1F* h=(TH1F*)coll->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
 799       h->SetName(Form("h%d%d",i,ic));
 800       if(!hallcntr) {
 801         hallcntr=(TH1F*)h->Clone("hallcntr");
 802         hallcntr->Sumw2();
 803       } else {
 804         hallcntr->Add(h);
 805       }
 806
 807       nevents080+=h->Integral();
 808
 809       if(superimpose){
 810         hh=(TH1F*)colle->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
 811         hh->SetName(Form("hh%d%d",i,ic));
 812         if(!hhallcntr) {
 813           hhallcntr=(TH1F*)hh->Clone("hhallcntr");
 814           hhallcntr->Sumw2();
 815         }else hhallcntr->Add(hh);
 816         nnevents080+=hh->Integral();
 817
 818       }
 819     }
 820
 821     for(Int_t ic=0;ic<ncentr;ic++){
 822
 823       TH1F* h=(TH1F*)coll->Get("run",Form("centralityclass:%d_%d",ic*10,ic*10+10));
 824       h->SetName(Form("h%d%d",i,ic));
 825       h->Sumw2();
 826
 827       if(ic>=0 && ic<=1){ //0-20
 828         if(!h020) {
 829           h020=(TH1F*)h->Clone(Form("h020%s",coll->GetName()));
 830           h020->SetTitle(Form("Centrality 0-20 %s",coll->GetName()));
 831           if(superimpose){
 832             hh020=(TH1F*)hh->Clone(Form("hh020%s",coll->GetName()));
 833             hh020->SetTitle(Form("Centrality 0-20 %s",coll->GetName()));
 834           }
 835         }
 836         else {
 837           h020->Add(h);
 838           if(superimpose)hh020->Add(hh);
 839         }
 840       }
 841       if(ic>=2 && ic<=7){ //20-80
 842         if(!h2080) {
 843           h2080=(TH1F*)h->Clone(Form("h2080%s",coll->GetName()));
 844           h2080->SetTitle(Form("Centrality 20-80 %s",coll->GetName()));
 845           if(superimpose){
 846             hh2080=(TH1F*)hh->Clone(Form("hh2080%s",coll->GetName()));
 847             hh2080->SetTitle(Form("Centrality 20-80 %s",coll->GetName()));
 848           }
 849         }
 850         else {
 851           h2080->Add(h);
 852           if(superimpose)hh2080->Add(hh);
 853         }
 854
 855       }
 856
 857       h->Divide(hallcntr);
 858
 859       ccent->cd(ic+1);
 860       h->GetYaxis()->SetRangeUser(0.,0.15);
 861       h->DrawClone();
 862       /*
 863       if(ic==0&&i==0){
 864         spare->cd();
 865         h->Draw();
 866       }
 867       */
 868       // ccent->cd(1);
 869       // h->SetLineColor(ic+1);
 870       // if(ic==0)h->DrawClone();
 871       // else h->DrawClone("sames");
 872     }
 873
 874     ccent->cd(ncentr+1);
 875     h020->Divide(hallcntr);
 876     h020->DrawClone();
 877     if(superimpose){
 878       hh020->Divide(hhallcntr);
 879       hh020->SetLineColor(2);
 880       hh020->SetMarkerColor(2);
 881       hh020->DrawClone("sames");
 882     }
 883     TCanvas* cv020=new TCanvas(Form("cv020-%d",i),"0-20% vs run number",1400,600);
 884     cv020->cd();
 885     h020->GetYaxis()->SetRangeUser(0.,1.);
 886     h020->DrawClone();
 887     if(superimpose)hh020->DrawClone("sames");
 888     cv020->SaveAs(Form("cv020-%d.png",i));
 889
 890     ccent->cd(ncentr+2);
 891     h2080->Divide(hallcntr);
 892     h2080->DrawClone();
 893     if(superimpose){
 894       hh2080->Divide(hhallcntr);
 895       hh2080->SetLineColor(2);
 896       hh2080->SetMarkerColor(2);
 897       hh2080->DrawClone("sames");
 898     }
 899
 900     TCanvas* cv2080=new TCanvas(Form("cv2080-%d",i),"20-80% vs run number",1400,600);
 901     cv2080->cd();
 902     h2080->GetYaxis()->SetRangeUser(0.,1.);
 903     h2080->DrawClone();
 904     if(superimpose)hh2080->DrawClone("sames");
 905     cv2080->SaveAs(Form("cv2080-%d.png",i));
 906
 907     ccent->SaveAs(Form("%s%s.png",ccent->GetName(),textleg.Data()));
 908   }
 909
 910 }
 911
 912 void DrawProjections(TString partname="D0",TString h2dname="hMultvsPercentile",Int_t nsteps=0,TString direction="X",TString path="./"){
 913   gStyle->SetCanvasColor(0);
 914   gStyle->SetTitleFillColor(0);
 915   gStyle->SetStatColor(0);
 916   gStyle->SetPalette(1);
 917
 918   TString listname="outputCentrCheck";
 919
 920   TList* list;
 921   TH1F * hstat;
 922
 923   Bool_t isRead=ReadFile(list,hstat,listname,partname,path);
 924   if(!isRead) return;
 925   if(!list || !hstat){
 926     cout<<":-( null pointers..."<<endl;
 927     return;
 928   }
 929   Double_t nevents=hstat->Integral(5,6); //ev good vertex
 930
 931   TH2F* h2=(TH2F*)list->FindObject(h2dname);
 932   if(!h2){
 933     cout<<h2dname.Data()<<" not found"<<endl;
 934     return;
 935   }
 936   TCanvas* cv2d=new TCanvas("cv2d",h2->GetName());
 937   cv2d->cd();
 938   cv2d->SetLogz();
 939   cv2d->SetGrid();
 940   h2->Draw("colz");
 941   TPaveText *pvst=new TPaveText(0.6,0.2,0.9,0.7,"NDC");
 942   pvst->SetBorderSize(0);
 943   pvst->SetFillStyle(0);
 944   pvst->AddText("Bin -> Cont/nEvVtx");
 945
 946
 947   Int_t kbins=1;
 948   if(nsteps==0){
 949     if(direction=="X") nsteps=h2->GetNbinsY();
 950     if(direction=="Y") nsteps=h2->GetNbinsX();
 951   }
 952   else{
 953     if(direction=="X") kbins=h2->GetNbinsY()/nsteps;
 954     if(direction=="Y") kbins=h2->GetNbinsX()/nsteps;
 955   }
 956
 957   TCanvas *cvpj=new TCanvas(Form("cvpj%s%s",direction.Data(),h2dname.Data()),Form("cvpj%s",direction.Data()),1000,800);
 958   cvpj->Divide((Int_t)(nsteps/3)+1,3);
 959   TFile* fout=new TFile(Form("proj%s%s.root",direction.Data(),h2dname.Data()), "recreate");
 960   //Float_t maxx[nsteps];
 961   Float_t maxx[12]={9000,9000,6000,4000,2000,1400,800,500,200,100,40,25};
 962   Double_t integralpernev[nsteps];
 963
 964   for(Int_t i=0;i<nsteps;i++){
 965     TH1F* h=0x0;
 966     if(direction=="X")h=(TH1F*)h2->ProjectionX(Form("px%d",i),i+kbins,i+2*kbins);
 967     if(direction=="Y")h=(TH1F*)h2->ProjectionY(Form("py%d",i),i+kbins,i+2*kbins);
 968     integralpernev[i]=h->Integral()/nevents;
 969
 970     TPaveText *pvtxt=new TPaveText(0.6,0.6,0.9,0.9,"NDC");
 971     pvtxt->SetBorderSize(0);
 972     pvtxt->SetFillStyle(0);
 973     pvtxt->AddText(Form("%.0f - %.0f",h2->GetYaxis()->GetBinLowEdge((i+kbins)),h2->GetYaxis()->GetBinLowEdge((i+2*kbins))));
 974     pvst->AddText(Form("%.0f - %.0f -> %.2f",h2->GetYaxis()->GetBinLowEdge((i+kbins)),h2->GetYaxis()->GetBinLowEdge((i+2*kbins)),integralpernev[i]));
 975     cvpj->cd(i+1);
 976     h->GetXaxis()->SetRangeUser(0,maxx[i]);
 977     h->Draw();
 978     pvtxt->Draw();
 979     fout->cd();
 980     h->Write();
 981   }
 982   cvpj->SaveAs(Form("cvpj%s%s.png",direction.Data(),h2dname.Data()));
 983
 984   cv2d->cd();
 985   pvst->Draw();
 986   cv2d->SaveAs(Form("%s.png",h2->GetName()));
 987
 988 }