PWGLF/SPECTRA/PiKaPr/TestAOD/MainAnalysis.C

   1 #include <TString.h>
   2 #include <TStyle.h>
   3 #include <TSystem.h>
   4 #include <TFile.h>
   5 #include <iostream>
   6 #include "TFile.h"
   7 #include "TH1I.h"
   8 #include "TH2.h"
   9 #include "TCanvas.h"
  10 #include "TLegend.h"
  11 #include "TROOT.h"
  12
  13 TString Charge[]={"Pos","Neg"};
  14 TString Sign[]={"Plus","Minus"};
  15 TString Particle[]={"Pion","Kaon","Proton"};
  16 Int_t Color[3]={1,2,4};
  17 Int_t Marker[6]={20,21,22,24,25,26};
  18 Double_t projl[3]={0,20,80};
  19 Double_t proju[3]={10,40,90};
  20
  21
  22 void MainAnalysis() {
  23   //TString fold="5SigmaPIDFilterBit6";
  24   TString fold="5SigmaPID";
  25   Int_t ibinToCompare=-1;
  26
  27   //TString sname="Cent0to100_QVec0.0to100.0";
  28
  29   TString sname="Cent0to5_QVec0.0to100.0";ibinToCompare=0;
  30   //TString sname="Cent5to10_QVec0.0to100.0";ibinToCompare=1;
  31   //TString sname="Cent10to20_QVec0.0to100.0";ibinToCompare=2;
  32   //TString sname="Cent20to30_QVec0.0to100.0";ibinToCompare=3;
  33   //TString sname="Cent30to40_QVec0.0to100.0";ibinToCompare=4;
  34
  35   //TString sname="Cent20to40_QVec0.0to2.0";
  36   //TString sname="Cent20to40_QVec3.0to100.0";
  37
  38   TString dataFile = Form("outputAOD%s/Pt.AOD.1._data_ptcut_%s.root",fold.Data(),sname.Data());
  39   TString mcFile =Form("outputAOD%s/Pt.AOD.1._MC_%s.root",fold.Data(),sname.Data());
  40
  41   gSystem->Load("libCore.so");
  42   gSystem->Load("libGeom.so");
  43   gSystem->Load("libPhysics.so");
  44   gSystem->Load("libVMC");
  45   gSystem->Load("libTree");
  46   gSystem->Load("libProof");
  47   gSystem->Load("libMatrix");
  48   gSystem->Load("libSTEERBase");
  49   gSystem->Load("libESD");
  50   gSystem->Load("libAOD");
  51   gSystem->Load("libANALYSIS");
  52   gSystem->Load("libOADB");
  53   gSystem->Load("libANALYSISalice");
  54   gSystem->Load("libTENDER");
  55   gSystem->Load("libCORRFW");
  56   //gSystem->Load("libPWG0base");
  57   gSystem->Load("libMinuit");
  58   gSystem->Load("libPWGTools");
  59   gSystem->Load("libPWGLFSPECTRA");
  60   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  61   gStyle->SetPalette(1);
  62
  63   // Open root MC file and get classes
  64   cout << "Analysis Macro" << endl;
  65   cout << "  > Reading MC data" << endl;
  66   TFile *_mc = TFile::Open(mcFile.Data());
  67   AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) _mc->Get("SpectraHistos");
  68   AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) _mc->Get("Event Cuts");
  69   AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) _mc->Get("Track Cuts");
  70   // print info about mc track and Event cuts
  71   cout << " -- Info about MC -- "<< endl;
  72   ecuts_mc->PrintCuts();
  73   // tcuts_mc->PrintCuts();
  74   // proceed likewise for data
  75   TFile *_data = TFile::Open(dataFile.Data());
  76   AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) _data->Get("SpectraHistos");
  77   AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
  78   AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) _data->Get("Track Cuts");
  79   // print info about track and Event cuts
  80   cout << " -- Info about data -- " << endl;
  81   ecuts_data->PrintCuts();
  82   tcuts_data->PrintCuts();
  83
  84   //dedx in data and MC
  85   TCanvas *cPIDSig=new TCanvas("cPIDSig","cPIDSig",700,500);
  86   cPIDSig->Divide(2,2);
  87   cPIDSig->cd(1);
  88   TH2F *PIDSig_data = (TH2F*)((TH2F*)hman_data->GetPIDHistogram("histPIDTPC"))->Clone();
  89   gPad->SetLogz();
  90   gPad->SetGridy();
  91   gPad->SetGridx();
  92   PIDSig_data->DrawClone("colz");
  93   cPIDSig->cd(2);
  94   TH2F *PIDSig_mc = (TH2F*)((TH2F*)hman_mc->GetPIDHistogram("histPIDTPC"))->Clone();
  95   gPad->SetLogz();
  96   gPad->SetGridy();
  97   gPad->SetGridx();
  98   PIDSig_mc->DrawClone("colz");
  99   cPIDSig->cd(3);
 100   TH2F *PIDSig_data = (TH2F*)((TH2F*)hman_data->GetPIDHistogram("histPIDTOF"))->Clone();
 101   gPad->SetLogz();
 102   gPad->SetGridy();
 103   gPad->SetGridx();
 104   PIDSig_data->DrawClone("colz");
 105   cPIDSig->cd(4);
 106   TH2F *PIDSig_mc = (TH2F*)((TH2F*)hman_mc->GetPIDHistogram("histPIDTOF"))->Clone();
 107   gPad->SetLogz();
 108   gPad->SetGridy();
 109   gPad->SetGridx();
 110   PIDSig_mc->DrawClone("colz");
 111
 112   //nsig in data and MC
 113   for(Int_t ipart=0;ipart<3;ipart++){
 114     TCanvas *cnsig=new TCanvas(Form("cnsig%s",Particle[ipart].Data()),Form("cnsig%s",Particle[ipart].Data()),700,500);
 115     cnsig->Divide(2,3);
 116     cnsig->cd(1);
 117     TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTPC",Particle[ipart].Data())))->Clone();
 118     nsig_data->SetXTitle("Pt (GeV/c)");
 119     gPad->SetLogz();
 120     gPad->SetGridy();
 121     gPad->SetGridx();
 122     nsig_data->DrawClone("colz");
 123     cnsig->cd(2);
 124     TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTPC",Particle[ipart].Data())))->Clone();
 125     nsig_mc->SetXTitle("Pt (GeV/c)");
 126     gPad->SetLogz();
 127     gPad->SetGridy();
 128     gPad->SetGridx();
 129     nsig_mc->DrawClone("colz");
 130     cnsig->cd(3);
 131     TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTOF",Particle[ipart].Data())))->Clone();
 132     nsig_data->SetXTitle("Pt (GeV/c)");
 133     gPad->SetLogz();
 134     gPad->SetGridy();
 135     gPad->SetGridx();
 136     nsig_data->DrawClone("colz");
 137     cnsig->cd(4);
 138     TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTOF",Particle[ipart].Data())))->Clone();
 139     nsig_mc->SetXTitle("Pt (GeV/c)");
 140     gPad->SetLogz();
 141     gPad->SetGridy();
 142     gPad->SetGridx();
 143     nsig_mc->DrawClone("colz");
 144     cnsig->cd(5);
 145     TH2F *nsig_data = (TH2F*)((TH2F*)hman_data->GetNSigHistogram(Form("histNSig%sPtTPCTOF",Particle[ipart].Data())))->Clone();
 146     nsig_data->SetXTitle("Pt (GeV/c)");
 147     gPad->SetLogz();
 148     gPad->SetGridy();
 149     gPad->SetGridx();
 150     nsig_data->DrawClone("colz");
 151     cnsig->cd(6);
 152     TH2F *nsig_mc = (TH2F*)((TH2F*)hman_mc->GetNSigHistogram(Form("histNSig%sPtTPCTOF",Particle[ipart].Data())))->Clone();
 153     nsig_mc->SetXTitle("Pt (GeV/c)");
 154     gPad->SetLogz();
 155     gPad->SetGridy();
 156     gPad->SetGridx();
 157     nsig_mc->DrawClone("colz");
 158   }
 159   //qvec in data and MC
 160   for (Int_t icharge=0;icharge<2;icharge++){
 161     TCanvas *cqvec=new TCanvas(Form("cqvec%s",Charge[icharge].Data()),Form("cqvec%s",Charge[icharge].Data()),700,500);
 162     cqvec->Divide(2,1);
 163     cqvec->cd(1);
 164     TString hname=Form("histq%s",Charge[icharge].Data());
 165     TH2F *qvec_data = (TH2F*)((TH2F*)hman_data->GetqVecHistogram(hname.Data()))->Clone();
 166     gPad->SetLogz();
 167     gPad->SetGridy();
 168     gPad->SetGridx();
 169     qvec_data->DrawClone("colz");
 170     cqvec->cd(2);
 171     TH2F *qvec_mc = (TH2F*)((TH2F*)hman_mc->GetqVecHistogram(hname.Data()))->Clone();
 172     gPad->SetLogz();
 173     gPad->SetGridy();
 174     gPad->SetGridx();
 175     qvec_mc->DrawClone("colz");
 176
 177     TCanvas *cProjqvec=new TCanvas(Form("cProjqvec%s",Charge[icharge].Data()),Form("cProjqvec%s",Charge[icharge].Data()),700,500);
 178     cProjqvec->Divide(3,2);
 179     for(Int_t iproj=0;iproj<3;iproj++){
 180       TH1F *proj_data=(TH1F*)qvec_data->ProjectionX("data",qvec_data->GetYaxis()->FindBin(projl[iproj]),qvec_data->GetYaxis()->FindBin(proju[iproj]));
 181       if(proj_data->GetEntries()==0)continue;
 182       proj_data->Scale(1/proj_data->GetEntries());
 183       proj_data->SetTitle(Form("data q%s [%.0f,%.0f]",Charge[icharge].Data(),projl[iproj],proju[iproj]));
 184       TH1F *proj_mc=(TH1F*)qvec_mc->ProjectionX("mc",qvec_mc->GetYaxis()->FindBin(projl[iproj]),qvec_mc->GetYaxis()->FindBin(proju[iproj]));
 185       proj_mc->Scale(1/proj_mc->GetEntries());
 186       proj_mc->SetTitle(Form("mc q%s [%.0f,%.0f]",Charge[icharge].Data(),projl[iproj],proju[iproj]));
 187       proj_mc->SetLineColor(2);
 188       cProjqvec->cd(iproj+1);
 189       gPad->SetGridy();
 190       gPad->SetGridx();
 191       proj_data->DrawClone();
 192       proj_mc->DrawClone("same");
 193       gPad->BuildLegend();
 194       proj_data->Divide(proj_mc);
 195       cProjqvec->cd(iproj+4);
 196       proj_data->DrawClone();
 197     }
 198   }
 199
 200   //efficiencies
 201   Printf("\n\n-> Calculating MC Correction Factors");
 202   TH1F *CorrFact[6];
 203   TCanvas *ceff=new TCanvas("ceff","ceff",700,500);
 204   for(Int_t icharge=0;icharge<2;icharge++){
 205     for(Int_t ipart=0;ipart<3;ipart++){
 206       Int_t index=ipart+3*icharge;
 207       //TString hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 208       //TString hname=Form("histPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 209       //TString hname=Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 210       TString hname=Form("histPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 211       Printf("Getting %s",hname.Data());
 212       CorrFact[index]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 213       CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 214       CorrFact[index]->SetTitle(Form("CorrFact %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 215       CorrFact[index]->SetMarkerStyle(Marker[index]);
 216       CorrFact[index]->SetMarkerColor(Color[ipart]);
 217       CorrFact[index]->SetLineColor(Color[ipart]);
 218       hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 219       Printf("... and divide it by %s",hname.Data());
 220       CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
 221       if(index==0)CorrFact[index]->DrawClone();
 222       else CorrFact[index]->DrawClone("same");
 223     }
 224   }
 225   TFile *fESD=new TFile("EffAlex/pionEffPbPb.root");
 226   TH1F *hEffESD=(TH1F*)fESD->Get("effMapPionTpcOnlyNeg0");
 227   hEffESD->DrawClone("same");
 228   gPad->BuildLegend();
 229
 230   //Normalization
 231   printf("\n\n-> Spectra Normalization");
 232   AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
 233   Double_t events_data =  ecuts_data->NumberOfEvents();
 234   Printf(": accepted events: %.0f",events_data);
 235
 236   //divide RAW for Correction Factor
 237   Printf("\n\n-> Using MC correction factor to correct RAW spectra");
 238   TH1F *Spectra[6];
 239   for(Int_t icharge=0;icharge<2;icharge++){
 240     for(Int_t ipart=0;ipart<3;ipart++){
 241       Int_t index=ipart+3*icharge;
 242       TString hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 243       printf("Getting %s",hname.Data());
 244       Spectra[index] =(TH1F*)((TH1F*) hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 245       Spectra[index]->SetName(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 246       Spectra[index]->SetTitle(Form("Spectra %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 247       Spectra[index]->SetMarkerStyle(Marker[index]);
 248       Spectra[index]->SetMarkerColor(Color[ipart]);
 249       Spectra[index]->SetLineColor(Color[ipart]);
 250       Printf("... and divide it by %s",hname.Data());
 251       Spectra[index]->Divide(CorrFact[index]);//////////////////////////////////////////////////////////////////////////////////////////FIXME
 252       // if(index!=3)Spectra[index]->Divide(CorrFact[index]);
 253       // else{
 254       //        Spectra[index]=AliPWGHistoTools::MyDivideHistosDifferentBins(Spectra[index],hEffESD);
 255       // }
 256       Spectra[index]->Scale(1./events_data,"width");//NORMALIZATION
 257     }
 258   }
 259
 260
 261   //Geant/Fluka Correction
 262   Printf("\nGF correction for Kaons");
 263   TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
 264   TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
 265   TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
 266   TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
 267   for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){
 268     Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
 269     Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
 270     Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
 271     Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
 272     Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
 273     Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
 274   }
 275
 276   //Geant Fluka for P
 277   Printf("\nGF correction for Protons");
 278   const Int_t kNCharge=2;
 279   Int_t kPos=0;
 280   Int_t kNeg=1;
 281   TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
 282   TH2D * hCorrFluka[kNCharge];
 283   TH2D * hCorrFluka[2];
 284   hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
 285   hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
 286   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 287     Int_t nbins = Spectra[2]->GetNbinsX();
 288     Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
 289     for(Int_t ibin = 0; ibin < nbins; ibin++){
 290       Float_t pt = Spectra[2]->GetBinCenter(ibin);
 291       Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
 292       Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
 293       if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 294       if (pt > maxPtCorrection) pt = maxPtCorrection;
 295       Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
 296       //cout<<correction<<"     charge "<<icharge<<endl;
 297       if(icharge==0){
 298         if (correction != 0) {// If the bin is empty this is a  0
 299           Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction);
 300           Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError  (ibin)*correction);
 301         }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 302           cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries," << endl;
 303           cout << " Bin content: " << Spectra[2]->GetBinContent(ibin)  << endl;
 304         }
 305       }
 306       if(icharge==1){
 307         if (correction != 0) {// If the bin is empty this is a  0
 308           Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction);
 309           Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError  (ibin)*correction);
 310         }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 311           cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries," << endl;
 312           cout << " Bin content: " << Spectra[5]->GetBinContent(ibin)  << endl;
 313         }
 314       }
 315     }
 316   }
 317
 318
 319   DCACorrection(Spectra,hman_data,hman_mc);
 320
 321   TH1F *MCTruth[6];
 322   for(Int_t icharge=0;icharge<2;icharge++){
 323     for(Int_t ipart=0;ipart<3;ipart++){
 324       Int_t index=ipart+3*icharge;
 325       hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 326       MCTruth[index]=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),0.9,1.1))->Clone();
 327       MCTruth[index]->Scale(1./events_data,"width");//NORMALIZATION
 328     }
 329   }
 330
 331   //Contamination
 332   Printf("\n\n-> Contamination from MC");
 333   TH1F *Cont[6];
 334   TCanvas *ccont=new TCanvas("ccont","ccont",700,500);
 335   for(Int_t icharge=0;icharge<2;icharge++){
 336     for(Int_t ipart=0;ipart<3;ipart++){
 337       Int_t index=ipart+3*icharge;
 338       TString hname=Form("histPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 339       Printf("Getting %s",hname.Data());
 340       Cont[index] =(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 341       Cont[index]->SetName(Form("Cont_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 342       Cont[index]->SetTitle(Form("RecTrue/RecSigma %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 343       Cont[index]->SetMarkerStyle(Marker[index]);
 344       Cont[index]->SetMarkerColor(Color[ipart]);
 345       Cont[index]->SetLineColor(Color[ipart]);
 346       hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 347       Printf("... and divide it by %s",hname.Data());
 348       Cont[index]->Divide((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1));
 349       if(index==0)Cont[index]->DrawClone();
 350       else Cont[index]->DrawClone("same");
 351       //Spectra[index]->Multiply(Cont[index]);
 352     }
 353   }
 354   gPad->BuildLegend();
 355
 356
 357   //Drawing Final Spectra
 358   TCanvas *cspectra=new TCanvas("cspectra","cspectra",700,500);
 359   for(Int_t icharge=0;icharge<2;icharge++){
 360     for(Int_t ipart=0;ipart<3;ipart++){
 361       Int_t index=ipart+3*icharge;
 362       if(index==0)Spectra[index]->DrawClone();
 363       else Spectra[index]->DrawClone("same");
 364       //MCTruth[index]->DrawClone("same");
 365     }
 366   }
 367   gPad->BuildLegend();
 368
 369   //saving spectra
 370   Printf("\n\n-> Saving spectra in Out%s_%s.root",sname.Data(),fold.Data());
 371   TFile * fout=new TFile(Form("results/Out_%s_%s.root",sname.Data(),fold.Data()),"RECREATE");
 372   for(Int_t icharge=0;icharge<2;icharge++){
 373     for(Int_t ipart=0;ipart<3;ipart++){
 374       Int_t index=ipart+3*icharge;
 375       Spectra[index]->Write();
 376       CorrFact[index]->Write();
 377       MCTruth[index]->Write();
 378     }
 379   }
 380
 381
 382   //if Bin 0-5% with no cut ratio with combined analysis
 383   if(ibinToCompare!=-1){
 384     TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",700,500);
 385     CratioComb->Divide(3,2);
 386     TString nameComb[6]={Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare),
 387                          Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)};
 388     TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120323.root");
 389     for(Int_t icharge=0;icharge<2;icharge++){
 390       for(Int_t ipart=0;ipart<3;ipart++){
 391         Int_t index=ipart+3*icharge;
 392         Spectra[index]->GetXaxis()->SetRangeUser(0,2.5);
 393         TH1F *hcomb=fComb->Get(nameComb[index].Data())->Clone();
 394         CratioComb->cd(ipart+1);
 395         gPad->SetGridy();
 396         gPad->SetGridx();
 397         if(icharge==0)Spectra[index]->DrawClone();
 398         else Spectra[index]->DrawClone("same");
 399         hcomb->DrawClone("same");
 400         Spectra[index]->Divide(hcomb);
 401         Spectra[index]->SetMaximum(1.3);
 402         Spectra[index]->SetMinimum(0.7);
 403         CratioComb->cd(ipart+4);
 404         gPad->SetGridy();
 405         gPad->SetGridx();
 406         if(icharge==0)Spectra[index]->DrawClone();
 407         else Spectra[index]->DrawClone("same");
 408       }
 409     }
 410   }
 411
 412   //comparison with charged hadron
 413   Printf("\n\n-> ChargedHadron comparison");
 414   TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500);
 415   cAllCh->Divide(1,4);
 416   TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("histPtRec",-1,-1))->Clone();
 417   hChHad_data->Scale(1./events_data,"width");//NORMALIZATION
 418   //fraction of sec in MC
 419   TH1F *hSecAllMC=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtGen",0,0))->Clone();
 420   TH1F *hAllMC=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtGen",0,1))->Clone();
 421   hSecAllMC->Divide(hAllMC);
 422   cAllCh->cd(1);
 423   hSecAllMC->DrawClone();
 424   cAllCh->cd(2);
 425   hChHad_data->DrawClone();
 426   for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){
 427     Double_t en=hChHad_data->GetBinContent(ibin);
 428     Double_t sec=hSecAllMC->GetBinContent(ibin);
 429     hChHad_data->SetBinContent(ibin,en-(en*sec*0.2));
 430     //Printf("%f %f %d",en,sec,ibin);
 431     //Printf("%f",hChHad_data->GetBinContent(ibin));
 432   }
 433   hChHad_data->DrawClone("lhistsame");
 434   //efficiency for primaries
 435   TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtRec",-1,-1))->Clone();
 436   hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("histPtGen",0.5,1.5))->Clone());
 437   cAllCh->cd(3);
 438   hEff_mc->DrawClone();
 439   hChHad_data->Divide(hEff_mc);
 440   cAllCh->cd(4);
 441   hChHad_data->Draw();
 442
 443   //Comparison of efficiency with TPCTOF ESD analysis
 444   Printf("\n\n-> Calculating Efficiency to be compared with ESD analysis");
 445   TH1F *EffTRUEPions;
 446   TH1F *EffSIGMAPions;
 447   TCanvas *cEffESD=new TCanvas("cEffESD","cEffESD",700,500);
 448   cEffESD->Divide(1,2);
 449   Int_t icharge=1;
 450   Int_t ipart=0;
 451   //using MC truth
 452   //TString hname=Form("histPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 453   TString hname=Form("histPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 454   Printf("Getting %s",hname.Data());
 455   EffTRUEPions=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 456   EffTRUEPions->SetName(Form("Eff TRUE_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 457   EffTRUEPions->SetTitle(Form("Eff TRUE %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 458   EffTRUEPions->SetMarkerStyle(Marker[icharge]);
 459   EffTRUEPions->SetMarkerColor(Color[ipart]);
 460   EffTRUEPions->SetLineColor(Color[ipart]);
 461   hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 462   Printf("... and divide it by %s",hname.Data());
 463   EffTRUEPions->Divide(EffTRUEPions,(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
 464   //using NSigma
 465   //hname=Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 466   hname=Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 467   Printf("Getting %s",hname.Data());
 468   EffSIGMAPions=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 469   EffSIGMAPions->SetName(Form("Eff SIGMA_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 470   EffSIGMAPions->SetTitle(Form("Eff SIGMA %s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 471   EffSIGMAPions->SetMarkerStyle(Marker[icharge]);
 472   EffSIGMAPions->SetMarkerColor(Color[ipart+1]);
 473   EffSIGMAPions->SetLineColor(Color[ipart+1]);
 474   hname=Form("histPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 475   Printf("... and divide it by %s",hname.Data());
 476   EffSIGMAPions->Divide(EffSIGMAPions,(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
 477   cEffESD->cd(1);
 478   //if(icharge==0)EffTRUEPions->DrawClone();
 479   //else EffTRUEPions->DrawClone("same");
 480   EffTRUEPions->DrawClone("lhist");
 481   EffSIGMAPions->DrawClone("lhistsame");
 482   hEffESD->Draw("lhistsame");
 483   gPad->BuildLegend();
 484   cEffESD->cd(2);
 485   TH1F *hRatioTRUE=AliPWGHistoTools::MyDivideHistosDifferentBins(EffTRUEPions,hEffESD);
 486   hRatioTRUE->Draw("lhist");
 487   TH1F *hRatioSIGMA=AliPWGHistoTools::MyDivideHistosDifferentBins(EffSIGMAPions,hEffESD);
 488   hRatioSIGMA->Draw("lhistsame");
 489
 490
 491   //Muon over Pion Ratio
 492   TCanvas *cMu=new TCanvas("cMu","cMu");
 493   TH1F *hMuOverPi[2];
 494   for(Int_t icharge=0;icharge<2;icharge++){
 495     TString hname=Form("histPtRecTruePrimaryMuon%s",Sign[icharge].Data());
 496     hMuOverPi[icharge]=(TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 497     hname=Form("histPtRecTruePrimaryPion%s",Sign[icharge].Data());
 498     hMuOverPi[icharge]->Divide((TH1F*)((TH1F*) hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone());
 499     if(icharge==0)hMuOverPi[icharge]->DrawClone();
 500     else hMuOverPi[icharge]->DrawClone("same");
 501   }
 502
 503 }
 504
 505
 506
 507 void DCACorrection(TH1F **Spectra, AliSpectraAODHistoManager* hman_data, AliSpectraAODHistoManager* hman_mc){
 508
 509   Double_t FitRange[2]={-3,3};
 510   Int_t nrebin=20;
 511   Printf("\DCACorr");
 512   TH1F *hcorrection[2];
 513   TCanvas *ccorrection=new TCanvas("DCAcorrection","DCAcorrection",700,500);
 514   ccorrection->Divide(2,1);
 515   TString sample[2]={"data","mc"};
 516   for(Int_t icharge=0;icharge<2;icharge++){
 517     for(Int_t ipart=0;ipart<3;ipart++){
 518       Int_t index=ipart+3*icharge;
 519       for(Int_t isample=0;isample<2;isample++){
 520         TCanvas *cDCA=new TCanvas(Form("cDCA%d%s",index,sample[isample].Data()),Form("cDCA%d%s",index,sample[isample].Data()),700,500);
 521         hcorrection[isample]=(TH1F*)Spectra[index]->Clone();
 522         hcorrection[isample]->Reset("all");
 523         cDCA->Divide(8,4);
 524         for(Int_t ibin_data=6;ibin_data<38;ibin_data++){
 525           Double_t lowedge=Spectra[index]->GetBinLowEdge(ibin_data);
 526           Double_t binwidth=Spectra[index]->GetBinWidth(ibin_data);
 527           if(isample==0)TH1F *hToFit =(TH1F*) ((TH1F*)hman_data->GetDCAHistogram1D(Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 528           if(isample==1)TH1F *hToFit =(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 529           TH1F *hmc1=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaPrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 530           TH1F *hmc2=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaSecondaryWeakDecay%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 531           TH1F *hmc3=(TH1F*) ((TH1F*)hman_mc->GetDCAHistogram1D(Form("histPtRecSigmaSecondaryMaterial%s%s",Particle[ipart].Data(),Sign[icharge].Data()),lowedge,lowedge+binwidth))->Clone();
 532           Double_t minentries=5;
 533           if(hToFit->GetEntries()<=minentries || hmc1->GetEntries()<=minentries || hmc2->GetEntries()<=minentries || hmc3->GetEntries()<=minentries)continue;
 534           hmc1->Rebin(nrebin);
 535           hmc2->Rebin(nrebin);
 536           hmc3->Rebin(nrebin);
 537           hToFit->Rebin(nrebin);
 538           //Data and MC can have different stat
 539           hToFit->Sumw2();
 540           hmc1->Sumw2();
 541           hmc2->Sumw2();
 542           hmc3->Sumw2();
 543           hToFit->Scale(1./hToFit->GetEntries());
 544           Double_t normMC=hmc1->GetEntries()+hmc2->GetEntries()+hmc3->GetEntries();
 545           hmc1->Scale(1./normMC);
 546           hmc2->Scale(1./normMC);
 547           hmc3->Scale(1./normMC);
 548           cDCA->cd(ibin_data-5);
 549           gPad->SetGridy();
 550           gPad->SetGridx();
 551           gPad->SetLogy();
 552           hToFit->DrawClone("lhist");
 553           hmc3->DrawClone("lhist");
 554           TObjArray *mc = new TObjArray(3);        // MC histograms are put in this array
 555           mc->Add(hmc1);
 556           mc->Add(hmc2);
 557           mc->Add(hmc3);
 558           TFractionFitter* fit = new TFractionFitter(hToFit,mc); // initialise
 559           fit->Constrain(0,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 560           fit->Constrain(1,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 561           fit->Constrain(2,0.0,1.0);               // constrain fraction 1 to be between 0 and 1
 562           fit->SetRangeX(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
 563           hToFit->GetXaxis()->SetRange(hToFit->GetXaxis()->FindBin(FitRange[0]),hToFit->GetXaxis()->FindBin(FitRange[1]));
 564           hToFit->SetTitle(Form("DCA distr - %s",sample[isample].Data()));
 565           Int_t status = fit->Fit();               // perform the fit
 566           cout << "fit status: " << status << endl;
 567           if (status == 0) {                       // check on fit status
 568             TH1F* result = (TH1F*) fit->GetPlot();
 569             TH1F* PrimMCPred=(TH1F*)fit->GetMCPrediction(0);
 570             TH1F* secStMCPred=(TH1F*)fit->GetMCPrediction(1);
 571             TH1F* secMCPred=(TH1F*)fit->GetMCPrediction(2);
 572             //Drawing section
 573             PrimMCPred->SetLineColor(2);
 574             secStMCPred->SetLineColor(6);
 575             secMCPred->SetLineColor(4);
 576             hToFit->DrawClone("lhist");
 577             result->SetTitle("Fit result");
 578             result->DrawClone("lhistsame");
 579             PrimMCPred->DrawClone("lhistsame");
 580             secStMCPred->DrawClone("lhistsame");
 581             secMCPred->DrawClone("lhistsame");
 582             Double_t v1=0,v2=0,v3=0;
 583             Double_t ev1=0,ev2=0,ev3=0;
 584             fit->GetResult(0,v1,ev1);
 585             fit->GetResult(1,v2,ev2);
 586             fit->GetResult(2,v3,ev3);
 587             Printf("\n\n\n\n\nv1:%f  v2:%f  v3:%f   ev1:%f  ev2:%f  ev3:%f   ",v1,v2,v3,ev1,ev2,ev3);
 588             hcorrection[isample]->SetBinContent(ibin_data,v1/(v1+v2+v3));
 589             //hcorrection[isample]->SetBinError(ibin_data,TMath::Sqrt(ev1*ev1+ev2*ev2+ev3*ev3));
 590             hcorrection[isample]->SetBinError(ibin_data,0);
 591           }
 592           else{
 593             hcorrection[isample]->SetBinContent(ibin_data,1);
 594             hcorrection[isample]->SetBinError(ibin_data,0);
 595           }
 596           Printf("deleting");
 597           delete hToFit;
 598         }
 599
 600         ccorrection->cd(icharge+1);
 601         gPad->SetGridy();
 602         gPad->SetGridx();
 603         hcorrection[isample]->SetTitle(Form("DCA corr %s %s %s",Particle[ipart].Data(),Charge[icharge].Data(),sample[isample].Data()));
 604         hcorrection[isample]->SetLineWidth(2);
 605         hcorrection[isample]->SetLineColor(Color[ipart]);
 606         hcorrection[isample]->SetLineStyle(isample+1);
 607         hcorrection[isample]->SetMarkerColor(Color[ipart]);
 608         hcorrection[isample]->GetXaxis()->SetRangeUser(0.2,2.5);
 609         if(ipart==0 && isample==0)hcorrection[isample]->DrawClone("lhist");
 610         else hcorrection[isample]->DrawClone("lhistsame");
 611         // smooth the DCA correction
 612         // TF1 *fitfun = new TF1("fitfun","[0]+[1]/x^[2]",0.2,0.8);
 613         // hcorrection[isample]->Fit(fitfun,"WRVMN","",0.35,0.8);
 614         // fitfun->SetLineWidth(1.5);
 615         // fitfun->SetLineColor(ipart+1);
 616         // fitfun->SetLineStyle(ipart+1);
 617         // fitfun->DrawClone("same");
 618         // for(Int_t ibin=1;ibin<30;ibin++){
 619         //      hcorrection[isample]->SetBinContent(ibin,fitfun->Eval(hcorrection[isample]->GetBinCenter(ibin)));
 620         // }
 621
 622       }
 623       Spectra[index]->Multiply(hcorrection[0]);//multiply for data
 624       //Spectra[index]->Divide(hcorrection[1]); //divide by Monte Carlo
 625     }
 626   }
 627   ccorrection->cd(1);
 628   gPad->BuildLegend();
 629   ccorrection->cd(2);
 630   gPad->BuildLegend();
 631
 632 }
 633
 634