PWGLF/SPECTRA/PiKaPr/TestAOD/HighLevelQA/JCGAnalysis.C

   1 /////////////////////////////////////////////////////
   2 // JCGAnalysis.C                                   //
   3 // based on MainAnalysis.C                         //
   4 // uses a new method of calculating the efficiency //
   5 /////////////////////////////////////////////////////
   6
   7 class AliPID;
   8
   9 // global constants
  10 const Int_t nCharge = 2;
  11 const Int_t nPart = 3;
  12 TString Charge[nCharge] = {"Pos","Neg"};
  13 TString Sign[nCharge] = {"Plus","Minus"};
  14 TString Particle[nPart] = {"Pion","Kaon","Proton"};
  15 Int_t Color[nPart] = {1,2,4};
  16 Int_t Marker[nPart*nCharge] = {20,21,22,24,25,26};
  17 Double_t Range[nPart] = {0.3,0.3,0.5}; // LowPt range for pi k p
  18 enum ECharge_t {kPositive,kNegative,kNCharges};
  19 TString Names[nCharge*nPart] = {"#pi^{+}",
  20                                 "K^{+}",
  21                                 "p",
  22                                 "#pi^{-}",
  23                                 "K^{-}",
  24                                 "#bar{p}"};
  25
  26 void JCGAnalysis()
  27 {
  28   // load libraries 'n such
  29   gSystem->Load("libCore.so");
  30   gSystem->Load("libPhysics.so");
  31   gSystem->Load("libTree");
  32   gSystem->Load("libMatrix");
  33   gSystem->Load("libSTEERBase");
  34   gSystem->Load("libESD");
  35   gSystem->Load("libAOD");
  36   gSystem->Load("libANALYSIS");
  37   gSystem->Load("libOADB");
  38   gSystem->Load("libANALYSISalice");
  39   gSystem->Load("libTENDER");
  40   gSystem->Load("libCORRFW");
  41   gSystem->Load("libPWGTools");
  42   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  43   gStyle->SetOptStat(000000);
  44   gStyle->SetPalette(1);
  45
  46   // Set masses
  47   Double_t mass[3];
  48   mass[0]   = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
  49   mass[1]   = TDatabasePDG::Instance()->GetParticle("K+")->Mass();
  50   mass[2] = TDatabasePDG::Instance()->GetParticle("proton")->Mass();
  51
  52   // cuts 'n such
  53   //                            0    1    2    3
  54   Double_t CentCutMin[4]= {     0,  30,  30,  30};
  55   Double_t CentCutMax[4]= {     5,  40,  40,  40};
  56   Double_t QvecCutMin[4]={      0,   0,   0, 1.5};
  57   Double_t QvecCutMax[4]={    100, 100, 0.4, 100};
  58   Double_t EtaMin[4]={       -0.8,-0.8,-0.8,-0.8};
  59   Double_t EtaMax[4]={        0.8, 0.8, 0.8, 0.8};
  60   Double_t Nsigmapid=3.;
  61   Double_t pt=10.;
  62   Double_t p=10.;
  63   Double_t y=.5;
  64   Double_t ptTofMatch=.6;
  65   UInt_t trkbit=1024;
  66   UInt_t trkbitQVector=1;
  67   Bool_t UseCentPatchAOD049=kFALSE;
  68   Double_t DCA=100000;
  69   UInt_t minNclsTPC=70;
  70   Int_t nrebin=0;
  71   TString opt="";
  72   Int_t icut=1; // chooses which cut to use
  73   if(icut==0)Int_t ibinToCompare=0;
  74   else Int_t ibinToCompare=4;
  75
  76   TString sname_data=Form("OutputAODSpectraTask_Data_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_%.1fSigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],Nsigmapid,trkbit);
  77   TString sname_data_1sig=Form("OutputAODSpectraTask_Data_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_1.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
  78   // for MC we do not use the cut on Qvector
  79   TString sname_mc=Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec0.0to100.0_Eta%.1fto%.1f_%.1fSigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],EtaMin[icut],EtaMax[icut],Nsigmapid,trkbit);
  80   TString sname_mc_5sig=Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec0.0to100.0_Eta%.1fto%.1f_5.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
  81   TString sname_mc_1sig=Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec0.0to100.0_Eta%.1fto%.1f_1.0SigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],EtaMin[icut],EtaMax[icut],trkbit);
  82
  83   // input and output files
  84   TString fold="No_Outliers";
  85   TString dataFile = Form("output/%s/AnResDATATrain12_NoOutliers.root",fold.Data());
  86   TString mcFile = Form("output/%s/AnResMCTrain27_NoOutliers.root",fold.Data());
  87   TFile * fout=new TFile(Form("results/%s/Res_%s_%s.root",fold.Data(),sname_data.Data(),fold.Data()),"RECREATE");
  88
  89   // Open root MC file and get classes
  90   TFile *f_mc = TFile::Open(mcFile.Data());
  91   TDirectoryFile *Dir_mc=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc.Data()));
  92   AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) Dir_mc->Get("SpectraHistos");
  93   AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) Dir_mc->Get("Event Cuts");
  94   AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) Dir_mc->Get("Track Cuts");
  95   TDirectoryFile *Dir_mc_5sig=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc_5sig.Data()));
  96   AliSpectraAODHistoManager* hman_mc_5sig = (AliSpectraAODHistoManager*) Dir_mc_5sig->Get("SpectraHistos");
  97   AliSpectraAODEventCuts* ecuts_mc_5sig = (AliSpectraAODEventCuts*) Dir_mc_5sig->Get("Event Cuts");
  98   AliSpectraAODTrackCuts* tcuts_mc_5sig = (AliSpectraAODTrackCuts*) Dir_mc_5sig->Get("Track Cuts");
  99   TDirectoryFile *Dir_mc_1sig=(TDirectoryFile*)f_mc->Get(Form("%s",sname_mc_1sig.Data()));
 100   AliSpectraAODHistoManager* hman_mc_1sig = (AliSpectraAODHistoManager*) Dir_mc_1sig->Get("SpectraHistos");
 101   AliSpectraAODEventCuts* ecuts_mc_1sig = (AliSpectraAODEventCuts*) Dir_mc_1sig->Get("Event Cuts");
 102   AliSpectraAODTrackCuts* tcuts_mc_1sig = (AliSpectraAODTrackCuts*) Dir_mc_1sig->Get("Track Cuts");
 103   // proceed likewise for data:
 104   // Open root DATA file and get classes
 105   TFile *f_data = TFile::Open(dataFile.Data());
 106   TDirectoryFile *Dir_data=(TDirectoryFile*)f_data->Get(Form("%s",sname_data.Data()));
 107   AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) Dir_data->Get("SpectraHistos");
 108   AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) Dir_data->Get("Event Cuts");
 109   AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) Dir_data->Get("Track Cuts");
 110   TDirectoryFile *Dir_data_1sig=(TDirectoryFile*)f_data->Get(Form("%s",sname_data_1sig.Data()));
 111   AliSpectraAODHistoManager* hman_data_1sig = (AliSpectraAODHistoManager*) Dir_data_1sig->Get("SpectraHistos");
 112   AliSpectraAODEventCuts* ecuts_data_1sig = (AliSpectraAODEventCuts*) Dir_data_1sig->Get("Event Cuts");
 113   AliSpectraAODTrackCuts* tcuts_data_1sig = (AliSpectraAODTrackCuts*) Dir_data_1sig->Get("Track Cuts");
 114
 115   //---------------------------------------------
 116   // plot nsigma distributions for presentation -
 117   //---------------------------------------------
 118   // TPC
 119   TH2F * hTPCnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTPC"))->Clone());
 120   hTPCnsigPion->SetTitle("TPC NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
 121   TCanvas * cTPCnsigPion = new TCanvas("cTPCnsigPion","cTPCnsigPion");
 122   gPad->SetLogz();
 123   hTPCnsigPion->GetXaxis()->SetRangeUser(0,5);
 124   hTPCnsigPion->GetYaxis()->SetRangeUser(-10,40);
 125   hTPCnsigPion->DrawCopy("COLZ");
 126   // TOF
 127   TH2F * hTOFnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTOF"))->Clone());
 128   hTOFnsigPion->SetTitle("TOF NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
 129   TCanvas * cTOFnsigPion = new TCanvas("cTOFnsigPion","cTOFnsigPion");
 130   gPad->SetLogz();
 131   hTOFnsigPion->GetYaxis()->SetRangeUser(-10,40);
 132   hTOFnsigPion->GetXaxis()->SetRangeUser(0,5);
 133   hTOFnsigPion->DrawCopy("COLZ");
 134   // TPCTOF
 135   TH2F * hTPCTOFnsigPion = new TH2F(*(TH2F*)((TH2F*)hman_data->GetNSigHistogram("hHistNSigPionTPCTOF"))->Clone());
 136   hTPCTOFnsigPion->SetTitle("TPC+TOF NSigma Distribution for Pions;p_{T} (GeV/c);n#sigma");
 137   TCanvas * cTPCTOFnsigPion = new TCanvas("cTPCTOFnsigPion","cTPCTOFnsigPion");
 138   gPad->SetLogz();
 139   hTPCTOFnsigPion->GetYaxis()->SetRangeUser(-5,40);
 140   hTPCTOFnsigPion->GetXaxis()->SetRangeUser(0,5);
 141   hTPCTOFnsigPion->DrawCopy("COLZ");
 142
 143
 144   //------------------------------------------------------
 145   // plot the MC True Generated Primaries for comparison -
 146   //------------------------------------------------------
 147   TH1F * MCTruthSpectra[nCharge*nPart];
 148   Double_t events_mc = ecuts_mc->NumberOfEvents();
 149   for (Int_t icharge=0; icharge<nCharge; icharge++)
 150     {
 151       for (Int_t ipart=0; ipart<nPart; ipart++)
 152         {
 153           Int_t index = ipart + 3*icharge;
 154           TString hname = Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 155           MCTruthSpectra[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),1,1)->Clone());
 156           MCTruthSpectra[index]->SetTitle("Generated Primaries Spectra (MC Truth);p_{T} (GeV/c);");
 157           MCTruthSpectra[index]->SetMarkerStyle(Marker[index]);
 158           MCTruthSpectra[index]->SetMarkerColor(Color[ipart]);
 159           MCTruthSpectra[index]->SetLineColor(Color[ipart]);
 160           MCTruthSpectra[index]->Scale(1./events_mc,"width");
 161         }
 162     }
 163   // neg/pos ratios
 164   TH1F * NegPosRatio[nPart];
 165   for (Int_t ipart=0; ipart<nPart; ipart++)
 166     {
 167       NegPosRatio[ipart] = (TH1F*)MCTruthSpectra[ipart+nPart]->Clone();
 168       NegPosRatio[ipart]->Divide(MCTruthSpectra[ipart]);
 169       NegPosRatio[ipart]->SetTitle("Neg/Pos;p_{T} (GeV/c);");
 170     }
 171   // drawing
 172   TCanvas * cMCTruthSpectra = new TCanvas("cMCTruthSpectra","cMCTruthSpectra");
 173   cMCTruthSpectra->Divide(1,2);
 174   TLegend * lMCTruthSpectra = new TLegend(.69,.69,.89,.99);
 175   lMCTruthSpectra->SetFillColor(0);
 176   TLegend * lPosNegRatios = new TLegend (.69,.69,.99,.99);
 177   lPosNegRatios->SetFillColor(0);
 178   for (Int_t icharge=0; icharge<nCharge; icharge++)
 179     {
 180       for (Int_t ipart=0; ipart<nPart; ipart++)
 181         {
 182           Int_t index = ipart + 3*icharge;
 183           cMCTruthSpectra->cd(1);
 184           gPad->SetLogy();
 185           if (index == 0) MCTruthSpectra[index]->DrawCopy("EP");
 186           else MCTruthSpectra[index]->DrawCopy("EPsame");
 187           lMCTruthSpectra->AddEntry(MCTruthSpectra[index],Names[index].Data(),"p");
 188
 189           cMCTruthSpectra->cd(2);
 190           NegPosRatio[ipart]->GetYaxis()->SetRangeUser(0.4,1.1);
 191           if (index == 0) NegPosRatio[ipart]->DrawCopy("hist][");
 192           else NegPosRatio[ipart]->DrawCopy("hist][same");
 193           if (icharge == 0) lPosNegRatios->AddEntry(NegPosRatio[ipart],Particle[ipart].Data(),"l");
 194         }
 195     }
 196   cMCTruthSpectra->cd(1);
 197   lMCTruthSpectra->DrawClone();
 198   cMCTruthSpectra->cd(2);
 199   lPosNegRatios->DrawClone();
 200
 201
 202   //-----------------------------------
 203   // Get the (normalized) raw spectra -
 204   //-----------------------------------
 205   TH1F * Spectra[nCharge*nPart];
 206   Double_t events_data =  ecuts_data->NumberOfEvents();
 207   Double_t events_data_1sig =  ecuts_data_1sig->NumberOfEvents();
 208   Double_t events_mc =  ecuts_mc->NumberOfEvents();
 209   Double_t events_mc_1sig = ecuts_mc_1sig->NumberOfEvents();
 210   for (Int_t icharge=0; icharge<nCharge; icharge++)
 211     {
 212       for (Int_t ipart=0; ipart<nPart; ipart++)
 213         {
 214           Int_t index = ipart + 3*icharge;
 215           TString hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 216           Spectra[index] = (TH1F*)((TH1F*)hman_data->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 217           Spectra[index]->SetTitle(Form("Spectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 218           Spectra[index]->SetMarkerStyle(Marker[index]);
 219           Spectra[index]->SetMarkerColor(Color[ipart]);
 220           Spectra[index]->SetLineColor(Color[ipart]);
 221           Spectra[index]->Scale(1./events_data,"width"); // normalization
 222         }
 223     }
 224   // Set bins outisde the range to contain zeros
 225   for (Int_t icharge=0; icharge<nCharge; icharge++)
 226     {
 227       for (Int_t ipart=0; ipart<nPart; ipart++)
 228         {
 229           Int_t index = ipart + 3*icharge;
 230           for (Int_t ibin=0; ibin<Spectra[index]->GetNbinsX(); ibin++)
 231             {
 232               if (Spectra[index]->GetBinCenter(ibin) < Range[ipart])
 233                 {
 234                   Spectra[index]->SetBinContent(ibin, 0);
 235                   Spectra[index]->SetBinError(ibin, 0);
 236                 }
 237             }
 238         }
 239     }
 240   // save the raw spectra for use later
 241   TH1F * RawSpectra[nCharge*nPart];
 242   for (Int_t icharge=0; icharge<nCharge; icharge++)
 243     {
 244       for (Int_t ipart=0; ipart<nPart; ipart++)
 245         {
 246           Int_t index = ipart + 3*icharge;
 247           RawSpectra[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
 248           RawSpectra[index]->SetTitle("Raw Spectra");
 249         }
 250     }
 251
 252   // plot the raw spectra from DATA
 253   // also include allch
 254   TCanvas * cRawSpectra = new TCanvas("cRawSpectra","cRawSpectra");
 255   TLegend * lRawSpectra = new TLegend(.69,.69,.99,.99);
 256   lRawSpectra->SetFillColor(0);
 257   TH1F * hRawAllCh = (TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
 258   hRawAllCh->Scale(1./ecuts_data->NumberOfEvents(),"width");
 259   hRawAllCh->SetTitle("Raw Spectra, Data;p_{T} (GeV/c);");
 260   hRawAllCh->SetMarkerColor(kGreen+3);
 261   hRawAllCh->SetLineColor(kGreen+3);
 262   hRawAllCh->SetMarkerStyle(34);
 263   hRawAllCh->DrawCopy();
 264   lRawSpectra->AddEntry(hRawAllCh,"All Particles","p");
 265   for (Int_t icharge=0; icharge<nCharge; icharge++)
 266     {
 267       for (Int_t ipart=0; ipart<nPart; ipart++)
 268         {
 269           Int_t index = ipart + 3*icharge;
 270           RawSpectra[index]->SetTitle("Raw Spectra, Data;p_{T} (GeV/c);");
 271           RawSpectra[index]->SetMarkerColor(Color[ipart]);
 272           RawSpectra[index]->SetLineColor(Color[ipart]);
 273           RawSpectra[index]->SetMarkerStyle(Marker[index]);
 274           RawSpectra[index]->DrawCopy("same");
 275           lRawSpectra->AddEntry(RawSpectra[index],Names[index].Data(),"p");
 276         }
 277     }
 278   lRawSpectra->DrawClone();
 279
 280   // also get MC raw for comparison purposes
 281   TH1F * MCSpectra[nCharge*nPart];
 282   for (Int_t icharge=0; icharge<nCharge; icharge++)
 283     {
 284       for (Int_t ipart=0; ipart<nPart; ipart++)
 285         {
 286           Int_t index = ipart + 3*icharge;
 287           TString hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 288           MCSpectra[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 289           MCSpectra[index]->SetTitle(Form("MCSpectra_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 290           MCSpectra[index]->SetMarkerStyle(Marker[index]);
 291           MCSpectra[index]->SetMarkerColor(Color[ipart]);
 292           MCSpectra[index]->SetLineColor(Color[ipart]);
 293           MCSpectra[index]->Scale(1./events_mc,"width"); // normalization
 294         }
 295     }
 296   // Set bins outisde the range to contain zeros
 297   for (Int_t icharge=0; icharge<nCharge; icharge++)
 298     {
 299       for (Int_t ipart=0; ipart<nPart; ipart++)
 300         {
 301           Int_t index = ipart + 3*icharge;
 302           for (Int_t ibin=0; ibin<MCSpectra[index]->GetNbinsX(); ibin++)
 303             {
 304               if (MCSpectra[index]->GetBinCenter(ibin) < Range[ipart])
 305                 {
 306                   MCSpectra[index]->SetBinContent(ibin, 0);
 307                   MCSpectra[index]->SetBinError(ibin, 0);
 308                 }
 309             }
 310         }
 311     }
 312   // save the raw spectra for use later
 313   TH1F * MCRawSpectra[nCharge*nPart];
 314   for (Int_t icharge=0; icharge<nCharge; icharge++)
 315     {
 316       for (Int_t ipart=0; ipart<nPart; ipart++)
 317         {
 318           Int_t index = ipart + 3*icharge;
 319           MCRawSpectra[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
 320           MCRawSpectra[index]->SetTitle("MC Raw Spectra");
 321         }
 322     }
 323
 324
 325   //----------------
 326   // MC Correction -
 327   //----------------
 328   // calculate the percentage of primaries reconstructed w/ 5 sigma PID
 329   Printf("\n\n-> Calculating MC Correction Factors");
 330   TH1F * CorrFact[nPart*nCharge];
 331   for (Int_t icharge=0; icharge<nCharge; icharge++)
 332     {
 333       for (Int_t ipart=0; ipart<nPart; ipart++)
 334         {
 335           Int_t index = ipart + 3*icharge;
 336           TString hname = Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 337           Printf("Getting %s",hname.Data());
 338           CorrFact[index] = new TH1F(*(TH1F*)((TH1F*) hman_mc_5sig->GetPtHistogram1D(hname.Data(),-1,-1))->Clone());
 339           CorrFact[index]->SetName(Form("CorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 340           CorrFact[index]->SetTitle("MC Correction Factor");
 341           CorrFact[index]->SetMarkerStyle(Marker[index]);
 342           CorrFact[index]->SetMarkerColor(Color[ipart]);
 343           CorrFact[index]->SetLineColor(Color[ipart]);
 344           hname=Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 345           Printf("... and divide it by %s",hname.Data());
 346           CorrFact[index]->Divide(CorrFact[index],(TH1F*)((TH1F*)hman_mc_5sig->GetPtHistogram1D(hname.Data(),1,1))->Clone(),1,1,"B");//binomial error
 347           if (index == 0)
 348             {
 349               TCanvas * cCorrFactMC = new TCanvas("cCorrFactMC","cCorrFactMC");
 350               TLegend * lCorrFactMC = new TLegend(.69,.69,.99,.99);
 351               lCorrFactMC->SetFillColor(0);
 352               CorrFact[index]->DrawCopy("EP");
 353             }
 354           else CorrFact[index]->DrawCopy("EPsame");
 355           lCorrFactMC->AddEntry(CorrFact[index],Names[index].Data(),"p");
 356         } // end loop over ipart
 357     } // end loop over icharge
 358   lCorrFactMC->DrawClone();
 359   // Correct spectra with MC correction factor
 360   for (Int_t icharge=0; icharge<nCharge; icharge++)
 361     {
 362       for (Int_t ipart=0; ipart<nPart; ipart++)
 363         {
 364           Int_t index = ipart + 3*icharge;
 365           Spectra[index]->Divide(CorrFact[index]);
 366         }
 367     }
 368   // save this version of the spectra for use later
 369   TH1F * SpectraStep1[nCharge*nPart];
 370   for (Int_t icharge=0; icharge<nCharge; icharge++)
 371     {
 372       for (Int_t ipart=0; ipart<nPart; ipart++)
 373         {
 374           Int_t index = ipart + 3*icharge;
 375           SpectraStep1[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
 376           SpectraStep1[index]->SetTitle("Spectra After MC Correction");
 377         }
 378     }
 379   // do the same with MCSpectra
 380   for (Int_t icharge=0; icharge<nCharge; icharge++)
 381     {
 382       for (Int_t ipart=0; ipart<nPart; ipart++)
 383         {
 384           Int_t index = ipart + 3*icharge;
 385           MCSpectra[index]->Divide(CorrFact[index]);
 386         }
 387     }
 388   // save this version of the spectra for use later
 389   TH1F * MCSpectraStep1[nCharge*nPart];
 390   for (Int_t icharge=0; icharge<nCharge; icharge++)
 391     {
 392       for (Int_t ipart=0; ipart<nPart; ipart++)
 393         {
 394           Int_t index = ipart + 3*icharge;
 395           MCSpectraStep1[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
 396           MCSpectraStep1[index]->SetTitle("MCSpectra After MC Correction");
 397         }
 398     }
 399
 400
 401   //--------------------------------------
 402   // Correct spectra for PID nsigma cuts -
 403   //--------------------------------------
 404   // for the TPC only, assume the nsigma distribution is purely gaussian, but account for the actual means and sigmas
 405   // the correction factor is the integral of a normal distribution from -3sqrt(2) to 3sqrt(2)
 406   TF1 * fGaussTPC[nPart];
 407   for (Int_t ipart=0; ipart<nPart; ipart++)
 408     {
 409       fGaussTPC[ipart] = new TF1(Form("fGaussTPC%s",Particle[ipart].Data()),"gaus",-25,25);
 410       fGaussTPC[ipart]->SetParameter(0,1); // set a default height for now - will be normalized later
 411     }
 412   // pions
 413   fGaussTPC[0]->SetParameter(1,0.04);
 414   fGaussTPC[0]->SetParameter(2,1.02);
 415   // kaons
 416   fGaussTPC[1]->SetParameter(1,-0.06);
 417   fGaussTPC[1]->SetParameter(2,0.97);
 418   // protons
 419   fGaussTPC[2]->SetParameter(1,0.25);
 420   fGaussTPC[2]->SetParameter(2,0.72);
 421   // normalization
 422   for (Int_t ipart=0; ipart<nPart; ipart++)
 423     {
 424       fGaussTPC[ipart]->SetParameter(0, 1./fGaussTPC[ipart]->Integral(-20,20));
 425     }
 426   // calculate correction factors
 427   Float_t PIDCorrFactTPC[nPart];
 428   for (Int_t ipart=0; ipart<nPart; ipart++)
 429     {
 430       PIDCorrFactTPC[ipart] = fGaussTPC[ipart]->Integral(-3*TMath::Sqrt(2), 3*TMath::Sqrt(2));
 431     }
 432   cout << "\n\n\n--- TPC ONLY CORRECTION FACTORS: ---\n"
 433        << "Pions: " << PIDCorrFactTPC[0] << endl
 434        << "Kaons: " << PIDCorrFactTPC[1] << endl
 435        << "Protons: " << PIDCorrFactTPC[2] << endl << endl << endl;
 436   // for the TPC+TOF range, get the correction factor from integrating the TPC+TOF nsigma distribution
 437   TFile * fTPCTOFin = TFile::Open("TPCTOFnsigma.root");
 438   TH1F * hTPCTOFnsigma[nPart];
 439   Float_t PIDCorrFactTPCTOF[nPart];
 440   for (Int_t ipart=0; ipart<nPart; ipart++)
 441     {
 442       hTPCTOFnsigma[ipart] = (TH1F*)((TH1F*)((TCanvas*)fTPCTOFin->Get("cTPCTOFnsigmaPredicted"))->FindObject(Form("hTPCTOFnsigma%s",Particle[ipart].Data())))->Clone();
 443       PIDCorrFactTPCTOF[ipart] = hTPCTOFnsigma[ipart]->Integral(hTPCTOFnsigma[ipart]->FindBin(0), hTPCTOFnsigma[ipart]->FindBin(3));
 444       PIDCorrFactTPCTOF[ipart] /= hTPCTOFnsigma[ipart]->Integral(hTPCTOFnsigma[ipart]->FindBin(0), hTPCTOFnsigma[ipart]->FindBin(20));
 445     }
 446   cout << "\n\n\n-- TPC+TOF CORRECTION FACTORS: ---\n"
 447        << "Pions: " << PIDCorrFactTPCTOF[0] << endl
 448        << "Kaons: " << PIDCorrFactTPCTOF[1] << endl
 449        << "Protons: " << PIDCorrFactTPCTOF[2] << endl << endl << endl;
 450
 451   // calculate seperate correction factors for MC b/c the parameters will be different...
 452   TF1 * fGaussTPCMC[nPart];
 453   for (Int_t ipart=0; ipart<nPart; ipart++)
 454     {
 455       fGaussTPCMC[ipart] = new TF1(Form("fGaussTPCMC%s",Particle[ipart].Data()),"gaus",-25,25);
 456       fGaussTPCMC[ipart]->SetParameter(0,1); // set a default height for now - will be normalized later
 457     }
 458   // pions
 459   fGaussTPCMC[0]->SetParameter(1,-0.03);
 460   fGaussTPCMC[0]->SetParameter(2,1.35);
 461   // kaons
 462   fGaussTPCMC[1]->SetParameter(1,0.14);
 463   fGaussTPCMC[1]->SetParameter(2,1.27);
 464   // protons
 465   fGaussTPCMC[2]->SetParameter(1,-0.17);
 466   fGaussTPCMC[2]->SetParameter(2,1.03);
 467   // normalization
 468   for (Int_t ipart=0; ipart<nPart; ipart++)
 469     {
 470       fGaussTPCMC[ipart]->SetParameter(0, 1./fGaussTPCMC[ipart]->Integral(-20,20));
 471     }
 472   // calculate correction factors
 473   Float_t PIDCorrFactTPCMC[nPart];
 474   for (Int_t ipart=0; ipart<nPart; ipart++)
 475     {
 476       PIDCorrFactTPCMC[ipart] = fGaussTPCMC[ipart]->Integral(-3*TMath::Sqrt(2), 3*TMath::Sqrt(2));
 477     }
 478   cout << "\n\n\n--- TPC ONLY CORRECTION FACTORS (MC): ---\n"
 479        << "Pions: " << PIDCorrFactTPCMC[0] << endl
 480        << "Kaons: " << PIDCorrFactTPCMC[1] << endl
 481        << "Protons: " << PIDCorrFactTPCMC[2] << endl << endl << endl;
 482   // for the TPC+TOF range, get the correction factor from integrating the TPC+TOF nsigma distribution
 483   TFile * fTPCTOFin = TFile::Open("TPCTOFnsigma.root");
 484   TH1F * hTPCTOFnsigmaMC[nPart];
 485   Float_t PIDCorrFactTPCTOFMC[nPart];
 486   for (Int_t ipart=0; ipart<nPart; ipart++)
 487     {
 488       hTPCTOFnsigmaMC[ipart] = (TH1F*)((TH1F*)((TCanvas*)fTPCTOFin->Get("cTPCTOFnsigmaPredictedMC"))->FindObject(Form("hTPCTOFnsigmaMC%s",Particle[ipart].Data())))->Clone();
 489       PIDCorrFactTPCTOFMC[ipart] = hTPCTOFnsigmaMC[ipart]->Integral(hTPCTOFnsigmaMC[ipart]->FindBin(0), hTPCTOFnsigmaMC[ipart]->FindBin(3));
 490       PIDCorrFactTPCTOFMC[ipart] /= hTPCTOFnsigmaMC[ipart]->Integral(hTPCTOFnsigmaMC[ipart]->FindBin(0), hTPCTOFnsigmaMC[ipart]->FindBin(20));
 491     }
 492   cout << "\n\n\n-- TPC+TOF CORRECTION FACTORS (MC): ---\n"
 493        << "Pions: " << PIDCorrFactTPCTOFMC[0] << endl
 494        << "Kaons: " << PIDCorrFactTPCTOFMC[1] << endl
 495        << "Protons: " << PIDCorrFactTPCTOFMC[2] << endl << endl << endl;
 496
 497
 498   // correct spectra, taking into account the range where the TOF is valid
 499   for (Int_t icharge=0; icharge<nCharge; icharge++)
 500     {
 501       for (Int_t ipart=0; ipart<nPart; ipart++)
 502         {
 503           Int_t index = ipart + 3*icharge;
 504           for (Int_t ibin=1; ibin<Spectra[index]->GetNbinsX(); ibin++)
 505             {
 506               if (ibin < Spectra[index]->FindBin(0.6))
 507                 Spectra[index]->SetBinContent(ibin,Spectra[index]->GetBinContent(ibin) / PIDCorrFactTPC[ipart]);
 508               else
 509                 Spectra[index]->SetBinContent(ibin,Spectra[index]->GetBinContent(ibin) / PIDCorrFactTPCTOF[ipart]);
 510             }
 511         } // end for loop over ipart
 512     } // end for loop over icharge
 513   // save this version of the spectra for use later
 514   TH1F * SpectraStep2[nCharge*nPart];
 515   for (Int_t icharge=0; icharge<nCharge; icharge++)
 516     {
 517       for (Int_t ipart=0; ipart<nPart; ipart++)
 518         {
 519           Int_t index = ipart + 3*icharge;
 520           SpectraStep2[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
 521           SpectraStep2[index]->SetTitle("Spectra After PID Correction");
 522         }
 523     }
 524   // also do this for MCSpectra
 525   for (Int_t icharge=0; icharge<nCharge; icharge++)
 526     {
 527       for (Int_t ipart=0; ipart<nPart; ipart++)
 528         {
 529           Int_t index = ipart + 3*icharge;
 530           for (Int_t ibin=1; ibin<MCSpectra[index]->GetNbinsX(); ibin++)
 531             {
 532               if (ibin < MCSpectra[index]->FindBin(0.6))
 533                 MCSpectra[index]->SetBinContent(ibin,MCSpectra[index]->GetBinContent(ibin) / PIDCorrFactTPCMC[ipart]);
 534               else
 535                 MCSpectra[index]->SetBinContent(ibin,MCSpectra[index]->GetBinContent(ibin) / PIDCorrFactTPCTOFMC[ipart]);
 536             }
 537         } // end for loop over ipart
 538     } // end for loop over icharge
 539   // save this version of the spectra for use later
 540   TH1F * MCSpectraStep2[nCharge*nPart];
 541   for (Int_t icharge=0; icharge<nCharge; icharge++)
 542     {
 543       for (Int_t ipart=0; ipart<nPart; ipart++)
 544         {
 545           Int_t index = ipart + 3*icharge;
 546           MCSpectraStep2[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
 547           MCSpectraStep2[index]->SetTitle("MCSpectra After PID Correction");
 548         }
 549     }
 550
 551
 552   //-------------------------------------------------------------------------
 553   // Correct spectra for contamination from other Pions, Kaons, and Protons -
 554   //-------------------------------------------------------------------------
 555   // calculate the correction factors first
 556   TH1F * ContCorrFact[nPart*nCharge];
 557   for (Int_t icharge=0; icharge<nCharge; icharge++)
 558     {
 559       for (Int_t ipart=0; ipart<nPart; ipart++)
 560         {
 561           Int_t index = ipart + 3*icharge;
 562           TString hname = Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 563           ContCorrFact[index] = new TH1F(*(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone());
 564           ContCorrFact[index]->SetTitle("Contamination Correction Factor");
 565           ContCorrFact[index]->SetMarkerStyle(Marker[index]);
 566           ContCorrFact[index]->SetMarkerColor(Color[ipart]);
 567           ContCorrFact[index]->SetLineColor(Color[ipart]);
 568           hname = Form("hHistPtRecSigma%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 569           ContCorrFact[index]->Divide(ContCorrFact[index],(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(),-1,-1,"B"); // binomial error
 570           if (index == 0)
 571             {
 572               TCanvas * cContCorrFact = new TCanvas("cContCorrFact","cContCorrFact");
 573               TLegend * lContCorrFact = new TLegend(.69,.69,.99,.99);
 574               ContCorrFact[index]->DrawCopy("EP");
 575             }
 576           else ContCorrFact[index]->DrawCopy("EPsame");
 577           lContCorrFact->AddEntry(ContCorrFact[index],Names[index].Data(),"p");
 578         } // end for loop over ipart
 579     } // end for loop over icharge
 580   lContCorrFact->DrawClone();
 581   // then use them to correct the spectra
 582   for (Int_t icharge=0; icharge<nCharge; icharge++)
 583     {
 584       for(Int_t ipart=0; ipart<nPart; ipart++)
 585         {
 586           Int_t index = ipart + 3*icharge;
 587           Spectra[index]->Multiply(ContCorrFact[index]);
 588         }
 589     }
 590   // save this version of the spectra for use later
 591   TH1F * SpectraStep3[nCharge*nPart];
 592   for (Int_t icharge=0; icharge<nCharge; icharge++)
 593     {
 594       for (Int_t ipart=0; ipart<nPart; ipart++)
 595         {
 596           Int_t index = ipart + 3*icharge;
 597           SpectraStep3[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
 598           SpectraStep3[index]->SetTitle("Spectra After Contamination Correction");
 599         }
 600     }
 601   // also do this for MCSpectra
 602   for (Int_t icharge=0; icharge<nCharge; icharge++)
 603     {
 604       for(Int_t ipart=0; ipart<nPart; ipart++)
 605         {
 606           Int_t index = ipart + 3*icharge;
 607           MCSpectra[index]->Multiply(ContCorrFact[index]);
 608         }
 609     }
 610   // save this version of the spectra for use later
 611   TH1F * MCSpectraStep3[nCharge*nPart];
 612   for (Int_t icharge=0; icharge<nCharge; icharge++)
 613     {
 614       for (Int_t ipart=0; ipart<nPart; ipart++)
 615         {
 616           Int_t index = ipart + 3*icharge;
 617           MCSpectraStep3[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
 618           MCSpectraStep3[index]->SetTitle("MCSpectra After Contamination Correction");
 619         }
 620     }
 621
 622
 623   //--------------------------------------------------
 624   // Correct spectra for the Geant-Fluka discrepancy -
 625   //--------------------------------------------------
 626   GFCorrection(Spectra,tcuts_data,fout);
 627   // save this version of the spectra for use later
 628   TH1F * SpectraStep4[nCharge*nPart];
 629   for (Int_t icharge=0; icharge<nCharge; icharge++)
 630     {
 631       for (Int_t ipart=0; ipart<nPart; ipart++)
 632         {
 633           Int_t index = ipart + 3*icharge;
 634           SpectraStep4[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
 635           SpectraStep4[index]->SetTitle("Spectra After Geant-Fluka Correction");
 636         }
 637     }
 638   // NOTE! DO NOT USE GF FOR MCSPECTRA!
 639
 640
 641   //--------------------------------------------------
 642   // Correct spectra for the TOF matching efficiency -
 643   //--------------------------------------------------
 644   //Matching efficiency in data and Monte Carlo
 645   TCanvas *cMatchingEff=new TCanvas("cMatchingEff","cMatchingEff",700,500);
 646   TH1F *hMatcEffPos_data=(TH1F*)tcuts_data->GetHistoNMatchedPos();
 647   hMatcEffPos_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedPos());
 648   hMatcEffPos_data->SetTitle("Matching Eff Pos - data");
 649   TH1F *hMatcEffNeg_data=(TH1F*)tcuts_data->GetHistoNMatchedNeg();
 650   hMatcEffNeg_data->Divide((TH1F*)tcuts_data->GetHistoNSelectedNeg());
 651   hMatcEffNeg_data->SetTitle("Matching Eff Neg - data");
 652   hMatcEffNeg_data->SetLineColor(2);
 653   /*
 654   // for my (John Groh) report, plot the matching efficiency just for data
 655   TCanvas * cMatchEffDataOnly = new TCanvas("cMatchEffDataOnly","cMatchEffDataOnly");
 656   hMatcEffPos_data->SetTitle("TOF Matching Efficiency;p_{T} (GeV/c);");
 657   hMatcEffPos_data->DrawCopy("lhist");
 658   hMatcEffNeg_data->SetTitle("TOF Matching Efficiency;p_{T} (GeV/c);");
 659   hMatcEffNeg_data->DrawCopy("lhistsame");
 660   TLegend * lMatcEffDataOnly = new TLegend(.69,.69,.99,.99);
 661   lMatcEffDataOnly->SetFillColor(0);
 662   lMatcEffDataOnly->AddEntry(hMatcEffPos_data,"Positive Particles","l");
 663   lMatcEffDataOnly->AddEntry(hMatcEffNeg_data,"Negative Particles","l");
 664   lMatcEffDataOnly->DrawClone();
 665   */
 666
 667   TH1F *hMatcEffPos_mc=(TH1F*)tcuts_mc->GetHistoNMatchedPos();
 668   hMatcEffPos_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedPos());
 669   hMatcEffPos_mc->SetTitle("Matching Eff Pos - mc");
 670   hMatcEffPos_mc->SetLineStyle(2);
 671   TH1F *hMatcEffNeg_mc=(TH1F*)tcuts_mc->GetHistoNMatchedNeg();
 672   hMatcEffNeg_mc->Divide((TH1F*)tcuts_mc->GetHistoNSelectedNeg());
 673   hMatcEffNeg_mc->SetTitle("Matching Eff Neg - mc");
 674   hMatcEffNeg_mc->SetLineColor(2);
 675   hMatcEffNeg_mc->SetLineStyle(2);
 676   cMatchingEff->Divide(1,2);
 677   cMatchingEff->cd(1);
 678   gPad->SetGridy();
 679   gPad->SetGridx();
 680   hMatcEffPos_data->DrawClone("lhist");
 681   hMatcEffNeg_data->DrawClone("lhistsame");
 682   hMatcEffPos_mc->DrawClone("lhistsame");
 683   hMatcEffNeg_mc->DrawClone("lhistsame");
 684   gPad->BuildLegend()->SetFillColor(0);
 685   hMatcEffPos_data->Divide(hMatcEffPos_mc);
 686   hMatcEffNeg_data->Divide(hMatcEffNeg_mc);
 687   cMatchingEff->cd(2);
 688   gPad->SetGridy();
 689   gPad->SetGridx();
 690   hMatcEffPos_data->DrawClone("lhist");
 691   hMatcEffNeg_data->DrawClone("lhistsame");
 692   TF1 *pol0MatchPos_data=new TF1("pol0MatchPos_data","pol0",2.5,5);
 693   hMatcEffPos_data->Fit("pol0MatchPos_data","MNR");
 694   pol0MatchPos_data->DrawClone("same");
 695   TF1 *pol0MatchNeg_data=new TF1("pol0MatchNeg_data","pol0",2.5,5);
 696   hMatcEffNeg_data->Fit("pol0MatchNeg_data","MNR");
 697   pol0MatchNeg_data->SetLineColor(2);
 698   pol0MatchNeg_data->DrawClone("same");
 699   Float_t ScalingMatchingPos=pol0MatchPos_data->GetParameter(0);
 700   Float_t ScalingMatchingNeg=pol0MatchNeg_data->GetParameter(0);
 701   //Correction spectra for matching efficiency
 702   for (Int_t ipart=0; ipart<nPart; ipart++)
 703     {
 704       for (Int_t ibin=1; ibin<=Spectra[ipart]->GetNbinsX(); ibin++)
 705         {
 706           Float_t ptspectra=Spectra[ipart]->GetBinCenter(ibin);
 707           if (ptspectra < tcuts_data->GetPtTOFMatching()) continue;
 708           Spectra[ipart]->SetBinContent(ibin, (Spectra[ipart]->GetBinContent(ibin) / ScalingMatchingPos));
 709           Spectra[ipart+3]->SetBinContent(ibin, (Spectra[ipart+3]->GetBinContent(ibin) / ScalingMatchingNeg));
 710         }
 711     }
 712   // save this version of the spectra for use later
 713   TH1F * SpectraStep5[nCharge*nPart];
 714   for (Int_t icharge=0; icharge<nCharge; icharge++)
 715     {
 716       for (Int_t ipart=0; ipart<nPart; ipart++)
 717         {
 718           Int_t index = ipart + 3*icharge;
 719           SpectraStep5[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
 720           SpectraStep5[index]->SetTitle("Spectra After Matching Efficiency Correction");
 721         }
 722     }
 723   // NOTE! DO NOT USE THE MATCHING EFFICIENCY CORRECTION FOR MCSPECTRA!
 724
 725
 726   //-----------------------------------------------------
 727   // Secondaries correction - use MC, not DCA from data -
 728   //------------------------------------------------------
 729   // (correction factor for secondaries) = (true primaries) / (primaries)
 730   TH1F * SecondaryCorrFact[nCharge*nPart];
 731   for (Int_t icharge=0; icharge<nCharge; icharge++)
 732     {
 733       for (Int_t ipart=0; ipart<nPart; ipart++)
 734         {
 735           Int_t index = ipart + 3*icharge;
 736           cout << "index = " << index << endl;
 737           TString hname = Form("hHistPtRecTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 738           SecondaryCorrFact[index] = (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone();
 739           hname = Form("hHistPtRecTrue%s%s",Particle[ipart].Data(),Sign[icharge].Data());
 740           SecondaryCorrFact[index]->Divide(SecondaryCorrFact[index], (TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(hname.Data(),-1,-1))->Clone(), 1, 1, "B"); // binomial error!
 741           SecondaryCorrFact[index]->SetName(Form("SecondaryCorrFact_%s%s",Particle[ipart].Data(),Sign[icharge].Data()));
 742           SecondaryCorrFact[index]->SetTitle("Correction Factor for Secondaries");
 743           SecondaryCorrFact[index]->SetStats(kFALSE);
 744           SecondaryCorrFact[index]->SetMarkerStyle(Marker[index]);
 745           SecondaryCorrFact[index]->SetMarkerColor(Color[ipart]);
 746           SecondaryCorrFact[index]->SetLineColor(Color[ipart]);
 747           if (index == 0)
 748             {
 749               TCanvas * cSecCorrFactMC = new TCanvas("cSecCorrFactMC","cSecCorrFactMC");
 750               TLegend * lSecCorrFactMC = new TLegend(.69,.69,.99,.99);
 751               lSecCorrFactMC->SetFillColor(0);
 752               SecondaryCorrFact[index]->DrawCopy("EP");
 753             }
 754           else SecondaryCorrFact[index]->DrawCopy("EPsame");
 755           lSecCorrFactMC->AddEntry(SecondaryCorrFact[index],Names[index],"ep");
 756         } // end for loop over ipart
 757     } // end for loop over icharge
 758   lSecCorrFactMC->DrawClone();
 759
 760   // correct Spectra with this
 761   for (Int_t icharge=0; icharge<nCharge; icharge++)
 762     {
 763       for (Int_t ipart=0; ipart<nPart; ipart++)
 764         {
 765           Int_t index = ipart + 3*icharge;
 766           Spectra[index]->Multiply(SecondaryCorrFact[index]);
 767         }
 768     }
 769
 770   // also correct MCSpectra with this
 771   for (Int_t icharge=0; icharge<nCharge; icharge++)
 772     {
 773       for (Int_t ipart=0; ipart<nPart; ipart++)
 774         {
 775           Int_t index = ipart + 3*icharge;
 776           MCSpectra[index]->Multiply(SecondaryCorrFact[index]);
 777         }
 778     }
 779   // save this version of the spectra for use later
 780   TH1F * SpectraStep6[nCharge*nPart];
 781   for (Int_t icharge=0; icharge<nCharge; icharge++)
 782     {
 783       for (Int_t ipart=0; ipart<nPart; ipart++)
 784         {
 785           Int_t index = ipart + 3*icharge;
 786           SpectraStep6[index] = new TH1F(*(TH1F*)Spectra[index]->Clone());
 787           SpectraStep6[index]->SetTitle("Spectra After Secondaries Correction");
 788         }
 789     }
 790   // also save this version of MCSpectra
 791   TH1F * MCSpectraStep6[nCharge*nPart];
 792   for (Int_t icharge=0; icharge<nCharge; icharge++)
 793     {
 794       for (Int_t ipart=0; ipart<nPart; ipart++)
 795         {
 796           Int_t index = ipart + 3*icharge;
 797           MCSpectraStep6[index] = new TH1F(*(TH1F*)MCSpectra[index]->Clone());
 798           MCSpectraStep6[index]->SetTitle("MCSpectra After Secondaries Correction");
 799         }
 800     }
 801
 802   //-------------------------
 803   // draw the final spectra -
 804   //-------------------------
 805   TCanvas * cSpectra = new TCanvas("cSpectra","cSpectra",100,100,700,500);
 806   gPad->SetGridy();
 807   gPad->SetLogy();
 808   TLegend * lSpectra = new TLegend(.69,.69,.99,.99);
 809   for (Int_t icharge=0; icharge<nCharge; icharge++)
 810     {
 811       for (Int_t ipart=0; ipart<nPart; ipart++)
 812         {
 813           Int_t index = ipart + 3*icharge;
 814           Spectra[index]->SetTitle("Final Spectra");
 815           if (index == 0) Spectra[index]->DrawClone("EP");
 816           else Spectra[index]->DrawClone("EPsame");
 817           lSpectra->AddEntry(Spectra[index],Names[index].Data(),"p");
 818         }
 819     }
 820   lSpectra->DrawClone();
 821
 822   //-----------------------------------------------
 823   // compare MCSpectra with hHistPtGenTruePrimary -
 824   //-----------------------------------------------
 825   TCanvas * cCompareMCRecWithMCTrue = new TCanvas("cCompareMCRecWithMCTrue","cCompareMCRecWithMCTrue",250,175,700,500);
 826   cCompareMCRecWithMCTrue->Divide(3,2);
 827   TH1F * hMCTruth[nPart*nCharge];
 828   for (Int_t icharge=0; icharge<nCharge; icharge++)
 829     {
 830       for (Int_t ipart=0; ipart<nPart; ipart++)
 831         {
 832           Int_t index = ipart + 3*icharge;
 833           hMCTruth[index] = new TH1F(*(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D(Form("hHistPtGenTruePrimary%s%s",Particle[ipart].Data(),Sign[icharge].Data()),1,1))->Clone());
 834           hMCTruth[index]->Scale(1./ecuts_mc->NumberOfEvents(),"width"); // normalize it first!
 835           hMCTruth[index]->SetMarkerStyle(Marker[index]);
 836           hMCTruth[index]->SetMarkerColor(Color[ipart]);
 837           hMCTruth[index]->SetLineColor(Color[ipart]);
 838         }
 839     }
 840   // set bins outside range to zero
 841   for (Int_t icharge=0; icharge<nCharge; icharge++)
 842     {
 843       for (Int_t ipart=0; ipart<nPart; ipart++)
 844         {
 845           Int_t index = ipart + 3*icharge;
 846           for (Int_t ibin=0; ibin<hMCTruth[index]->GetNbinsX(); ibin++)
 847             {
 848               if (hMCTruth[index]->GetBinCenter(ibin) < Range[ipart])
 849                 {
 850                   hMCTruth[index]->SetBinContent(ibin, 0);
 851                   hMCTruth[index]->SetBinError(ibin, 0);
 852                 }
 853             }
 854         }
 855     }
 856
 857
 858   TH1F * hRatioMCRecTrue[nPart*nCharge];
 859   for (Int_t ipart=0; ipart<nPart; ipart++)
 860     {
 861       for (Int_t icharge=0; icharge<nCharge; icharge++)
 862         {
 863           Int_t index = ipart + 3*icharge;
 864           hRatioMCRecTrue[index] = (TH1F*)MCSpectra[index]->Clone();
 865           hRatioMCRecTrue[index]->Divide(hMCTruth[index]);
 866           hRatioMCRecTrue[index]->SetTitle("(Corrected MC Spectra) / (MC Truth Spectra);p_{T} (GeV/c);");
 867           if (icharge == 1) hRatioMCRecTrue[index]->SetLineStyle(2);
 868         }
 869     }
 870   for (Int_t ipart=0; ipart<nPart; ipart++)
 871     {
 872       cCompareMCRecWithMCTrue->cd(ipart+1);
 873       gPad->SetLogy();
 874       MCSpectra[ipart]->SetMarkerColor(kGreen+3);
 875       MCSpectra[ipart+nPart]->SetMarkerColor(kGreen+2);
 876       MCSpectra[ipart]->SetMarkerSize(0.75);
 877       MCSpectra[ipart+nPart]->SetMarkerSize(0.75);
 878       MCSpectra[ipart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
 879       MCSpectra[ipart+nPart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
 880       hMCTruth[ipart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
 881       hMCTruth[ipart+nPart]->SetTitle(Form("MC Corrected and Truth Spectra, %ss;p_{T} (GeV/c);",Particle[ipart].Data()));
 882       hMCTruth[ipart]->SetMarkerSize(0.75);
 883       hMCTruth[ipart]->SetMarkerSize(0.75);
 884       hMCTruth[ipart]->DrawCopy("EP");
 885       hMCTruth[ipart+nPart]->DrawCopy("EPsame");
 886       MCSpectra[ipart]->DrawCopy("EPsame");
 887       MCSpectra[ipart+nPart]->DrawCopy("EPsame");
 888       TLegend * lCompare = new TLegend(.69,.59,.99,.91);
 889       lCompare->SetFillColor(0);
 890       lCompare->AddEntry(hMCTruth[ipart],Form("%s, Truth",Names[ipart].Data()),"p");
 891       lCompare->AddEntry(hMCTruth[ipart+nPart],Form("%s, Truth",Names[ipart+nPart].Data()),"p");
 892       lCompare->AddEntry(MCSpectra[ipart],Form("%s, Corrected",Names[ipart].Data()),"p");
 893       lCompare->AddEntry(MCSpectra[ipart+nPart],Form("%s, Corrected",Names[ipart+nPart].Data()),"p");
 894       lCompare->DrawClone();
 895
 896       cCompareMCRecWithMCTrue->cd(ipart+1+nPart);
 897       hRatioMCRecTrue[ipart]->GetYaxis()->SetRangeUser(.97,1.03);
 898       hRatioMCRecTrue[ipart]->DrawCopy("hist][");
 899       hRatioMCRecTrue[ipart+nPart]->DrawCopy("hist][same");
 900       TLegend * lCompareRatio = new TLegend(.79,.79,.99,.99);
 901       lCompareRatio->SetFillColor(0);
 902       lCompareRatio->AddEntry(hRatioMCRecTrue[ipart],Names[ipart].Data(),"l");
 903       lCompareRatio->AddEntry(hRatioMCRecTrue[ipart+nPart],Names[ipart+nPart].Data(),"l");
 904       lCompareRatio->DrawClone();
 905     }
 906
 907
 908   //---------------------------------------------------------
 909   // plot the evolution of the spectra afte each correction -
 910   //---------------------------------------------------------
 911   TCanvas * cEvolutionOfSpectra = new TCanvas("cEvolutionOfSpectra","cEvolutionOfSpectra",300,200,700,500);
 912   cEvolutionOfSpectra->Divide(4,2);
 913   TLegend * lEvolutionOfSpectra = new TLegend(.69,.69,.99,.99);
 914   for (Int_t icharge=0; icharge<nCharge; icharge++)
 915     {
 916       for (Int_t ipart=0; ipart<nPart; ipart++)
 917         {
 918           Int_t index = ipart + 3*icharge;
 919           lEvolutionOfSpectra->AddEntry(Spectra[index],Names[index].Data(),"p");
 920         }
 921     }
 922   cEvolutionOfSpectra->cd(1);
 923   gPad->SetLogy();
 924   for (Int_t icharge=0; icharge<nCharge; icharge++)
 925     {
 926       for (Int_t ipart=0; ipart<nPart; ipart++)
 927         {
 928           Int_t index = ipart + 3*icharge;
 929           RawSpectra[index]->SetMarkerSize(.5);
 930           if (index == 0) RawSpectra[index]->DrawClone("EP");
 931           else RawSpectra[index]->DrawClone("EPsame");
 932         }
 933     }
 934   lEvolutionOfSpectra->DrawClone();
 935   cEvolutionOfSpectra->cd(2);
 936   gPad->SetLogy();
 937   for (Int_t icharge=0; icharge<nCharge; icharge++)
 938     {
 939       for (Int_t ipart=0; ipart<nPart; ipart++)
 940         {
 941           Int_t index = ipart + 3*icharge;
 942           SpectraStep1[index]->SetMarkerSize(.5);
 943           if (index == 0) SpectraStep1[index]->DrawClone("EP");
 944           else SpectraStep1[index]->DrawClone("EPsame");
 945         }
 946     }
 947   lEvolutionOfSpectra->DrawClone();
 948   cEvolutionOfSpectra->cd(3);
 949   gPad->SetLogy();
 950   for (Int_t icharge=0; icharge<nCharge; icharge++)
 951     {
 952       for (Int_t ipart=0; ipart<nPart; ipart++)
 953         {
 954           Int_t index = ipart + 3*icharge;
 955           SpectraStep2[index]->SetMarkerSize(.5);
 956           if (index == 0) SpectraStep2[index]->DrawClone("EP");
 957           else SpectraStep2[index]->DrawClone("EPsame");
 958         }
 959     }
 960   lEvolutionOfSpectra->DrawClone();
 961   cEvolutionOfSpectra->cd(4);
 962   gPad->SetLogy();
 963   for (Int_t icharge=0; icharge<nCharge; icharge++)
 964     {
 965       for (Int_t ipart=0; ipart<nPart; ipart++)
 966         {
 967           Int_t index = ipart + 3*icharge;
 968           SpectraStep3[index]->SetMarkerSize(.5);
 969           if (index == 0) SpectraStep3[index]->DrawClone("EP");
 970           else SpectraStep3[index]->DrawClone("EPsame");
 971         }
 972     }
 973   lEvolutionOfSpectra->DrawClone();
 974   cEvolutionOfSpectra->cd(5);
 975   gPad->SetLogy();
 976   for (Int_t icharge=0; icharge<nCharge; icharge++)
 977     {
 978       for (Int_t ipart=0; ipart<nPart; ipart++)
 979         {
 980           Int_t index = ipart + 3*icharge;
 981           SpectraStep4[index]->SetMarkerSize(.5);
 982           if (index == 0) SpectraStep4[index]->DrawClone("EP");
 983           else SpectraStep4[index]->DrawClone("EPsame");
 984         }
 985     }
 986   lEvolutionOfSpectra->DrawClone();
 987   cEvolutionOfSpectra->cd(6);
 988   gPad->SetLogy();
 989   for (Int_t icharge=0; icharge<nCharge; icharge++)
 990     {
 991       for (Int_t ipart=0; ipart<nPart; ipart++)
 992         {
 993           Int_t index = ipart + 3*icharge;
 994           SpectraStep5[index]->SetMarkerSize(.5);
 995           if (index == 0) SpectraStep5[index]->DrawClone("EP");
 996           else SpectraStep5[index]->DrawClone("EPsame");
 997         }
 998     }
 999   lEvolutionOfSpectra->DrawClone();
1000   cEvolutionOfSpectra->cd(7);
1001   gPad->SetLogy();
1002   for (Int_t icharge=0; icharge<nCharge; icharge++)
1003     {
1004       for (Int_t ipart=0; ipart<nPart; ipart++)
1005         {
1006           Int_t index = ipart + 3*icharge;
1007           SpectraStep6[index]->SetMarkerSize(.5);
1008           if (index == 0) SpectraStep6[index]->DrawClone("EP");
1009           else SpectraStep6[index]->DrawClone("EPsame");
1010         }
1011     }
1012   lEvolutionOfSpectra->DrawClone();
1013   cEvolutionOfSpectra->cd(8);
1014   gPad->SetLogy();
1015   for (Int_t icharge=0; icharge<nCharge; icharge++)
1016     {
1017       for (Int_t ipart=0; ipart<nPart; ipart++)
1018         {
1019           Int_t index = ipart + 3*icharge;
1020           Spectra[index]->SetMarkerSize(.5);
1021           if (index == 0) Spectra[index]->DrawClone("EP");
1022           else Spectra[index]->DrawClone("EPsame");
1023         }
1024     }
1025   lEvolutionOfSpectra->DrawClone();
1026
1027
1028   //-----------------------------------------
1029   // Compare spectra with published spectra -
1030   //-----------------------------------------
1031   // also add another canvas that just plots the ratios b/c I can't use the spectra in my presentation (hasn't been published yet)
1032   TCanvas * cRatioCombined = new TCanvas("cRatioCombined","cRatioCombined");
1033   cRatioCombined->Divide(3,1);
1034   TCanvas *CratioComb=new TCanvas("CratioComb","CratioComb",150,125,700,500);
1035   CratioComb->Divide(3,2);
1036   TString nameComb[6] = {Form("cent%d_pion_plus",ibinToCompare),Form("cent%d_kaon_plus",ibinToCompare),Form("cent%d_proton_plus",ibinToCompare),
1037                          Form("cent%d_pion_minus",ibinToCompare),Form("cent%d_kaon_minus",ibinToCompare),Form("cent%d_proton_minus",ibinToCompare)};
1038   TFile *fComb=new TFile("Combined05/SPECTRA_COMB_20120730.root");
1039   TLegend * lComp[nPart];
1040   TLegend * lRatio[nPart];
1041   for(Int_t icharge=0; icharge<nCharge; icharge++)
1042     {
1043       for(Int_t ipart=0; ipart<nPart; ipart++)
1044         {
1045           Int_t index = ipart + 3*icharge;
1046           TH1F * htmp = (TH1F*)((TH1F*)Spectra[index])->Clone("");
1047           htmp->SetMarkerSize(1);
1048           htmp->GetXaxis()->SetRangeUser(0,5);
1049           TH1F * hcomb = (TH1F*)fComb->Get(nameComb[index].Data())->Clone();
1050           CratioComb->cd(ipart+1);
1051           gPad->SetGridy();
1052           gPad->SetGridx();
1053           htmp->SetTitle("Comparison with Combined Analysis;p_{T} (GeV/c);");
1054           if(icharge == 0) htmp->DrawCopy();
1055           else htmp->DrawCopy("same");
1056           hcomb->SetMarkerStyle(Marker[index]);
1057           hcomb->DrawCopy("same");
1058           if (icharge == 0) {lComp[ipart] = new TLegend(.49,.49,.89,.89); lComp[ipart]->SetFillColor(0);}
1059           lComp[ipart]->AddEntry(htmp,Form("%s, This Analysis",Names[index].Data()),"p");
1060           lComp[ipart]->AddEntry(hcomb,Form("%s, Combined Analysis",Names[index].Data()),"lep");
1061           if (icharge == 1) lComp[ipart]->DrawClone();
1062
1063           CratioComb->cd(ipart+nPart+1);
1064           gPad->SetGridy();
1065           gPad->SetGridx();
1066           htmp->Divide(hcomb);
1067           htmp->SetTitle("(This Analysis) / (Combined Analysis);p_{T} (GeV/c);");
1068           htmp->GetYaxis()->SetRangeUser(0.8, 1.2);
1069           if (icharge == 0)
1070             {
1071               CratioComb->cd(ipart+nPart+1);
1072               htmp->DrawClone("hist][");
1073               cRatioCombined->cd(ipart+1);
1074               htmp->DrawClone("hist][");
1075               lRatio[ipart] = new TLegend(.69,.69,.89,.89);
1076               lRatio[ipart]->SetFillColor(0);
1077             }
1078           else
1079             {
1080               htmp->SetLineStyle(2);
1081               CratioComb->cd(ipart+nPart+1);
1082               htmp->DrawClone("histsame][");
1083               cRatioCombined->cd(ipart+1);
1084               htmp->DrawClone("histsame][");
1085             }
1086           lRatio[ipart]->AddEntry(htmp,Names[index].Data(),"l");
1087           if (icharge == 1)
1088             {
1089               CratioComb->cd(ipart+nPart+1);
1090               lRatio[ipart]->DrawClone();
1091               cRatioCombined->cd(ipart+1);
1092               lRatio[ipart]->DrawClone();
1093             }
1094         }
1095     }
1096   //comparison with charged hadron
1097   Printf("\n\n-> ChargedHadron comparison");
1098   TH1F *hChHad_data=(TH1F*)((TH1F*)hman_data->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
1099   TH1F* hMatchCorrectionAllCh=(TH1F*)hMatcEffPos_data->Clone("hMatchCorrectionAllCh");
1100   hMatchCorrectionAllCh->Add(hMatcEffNeg_data); //correction for Matching efficiency!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1101   hMatchCorrectionAllCh->Scale(0.5);
1102   for(Int_t ibin=1;ibin<hChHad_data->GetNbinsX();ibin++){
1103     Float_t ptch=hChHad_data->GetBinCenter(ibin);
1104     if(ptch<tcuts_data->GetPtTOFMatching())continue;
1105     hChHad_data->SetBinContent(ibin,2*(hChHad_data->GetBinContent(ibin)/(ScalingMatchingPos+ScalingMatchingNeg)));
1106   }
1107   //fraction of sec in MC
1108   TH1F *hPrimRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
1109   TH1F *hAllRec_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRec",-1,-1))->Clone();
1110   for(Int_t ibin=1;ibin<=hChHad_data->GetNbinsX();ibin++){
1111     Double_t en_data=hChHad_data->GetBinContent(ibin);
1112     Double_t en_mc=hAllRec_mc->GetBinContent(ibin);
1113     Double_t prim_mc=hPrimRec_mc->GetBinContent(ibin);
1114     if(en_mc!=0)hChHad_data->SetBinContent(ibin,en_data-(en_data*(en_mc-prim_mc)*1.2/en_mc));
1115     //Printf("Before: %.0f After: %.0f  fraction: %.1f",en_data,hChHad_data->GetBinContent(ibin),hChHad_data->GetBinContent(ibin)/en_data);
1116   }
1117   hPrimRec_mc->Divide(hAllRec_mc);
1118   //efficiency for primaries
1119   TH1F *hEff_mc=(TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtRecPrimary",-1,-1))->Clone();
1120   hEff_mc->Divide((TH1F*)((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->Clone());
1121   TCanvas *cAllChFactors=new TCanvas("cAllChFactors","cAllChFactors",700,500);
1122   cAllChFactors->Divide(1,2);
1123   cAllChFactors->cd(1);
1124   gPad->SetGridy();
1125   gPad->SetGridx();
1126   hPrimRec_mc->SetTitle("Prim/All, charged hadron pure MC");
1127   hPrimRec_mc->DrawClone("lhist");
1128   gPad->BuildLegend()->SetFillColor(0);
1129   cAllChFactors->cd(2);
1130   gPad->SetGridy();
1131   gPad->SetGridx();
1132   hEff_mc->SetTitle("Efficiency for Primaries, charged hadron pure MC");
1133   hEff_mc->DrawClone("lhist");
1134   gPad->BuildLegend()->SetFillColor(0);
1135   //Printf("--------%f ",((TH1F*)hman_mc->GetPtHistogram1D("hHistPtGen",1,1))->GetEntries()/1.6/ecuts_mc->NumberOfEvents());
1136   hChHad_data->Scale(1./events_data,"width");//NORMALIZATION
1137   hChHad_data->Divide(hEff_mc);//Efficiency
1138   hChHad_data->Scale(1./(TMath::Abs(tcuts_data->GetEtaMin())+TMath::Abs(tcuts_data->GetEtaMax())));
1139   hChHad_data->SetTitle("All Ch from AOD");
1140   hChHad_data->SetName("AllCh");
1141   TCanvas *cAllCh=new TCanvas("cAllCh","cAllCh",700,500);
1142   cAllCh->Divide(1,2);
1143   cAllCh->cd(1);
1144   gPad->SetGridy();
1145   gPad->SetGridx();
1146   hChHad_data->GetYaxis()->SetRangeUser(0,2500);
1147   hChHad_data->DrawClone();
1148   TFile *fCh=new TFile("ChargedHadron/SPECTRA_UNID_110906.root");
1149   TH1F *hCh=fCh->Get("hpt_c0_5");
1150   hCh->SetTitle("All Ch from Jacek");
1151   hCh->SetMarkerColor(2);
1152   hCh->SetLineColor(2);
1153   //invariant yield
1154   for(Int_t ibin=0;ibin<hCh->GetNbinsX();ibin++){
1155     hCh->SetBinContent(ibin,hCh->GetBinContent(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
1156     hCh->SetBinError(ibin,hCh->GetBinError(ibin)*(2*TMath::Pi()*hCh->GetBinCenter(ibin)));
1157   }
1158   hCh->DrawClone("same");
1159   gPad->BuildLegend()->SetFillColor(0);
1160   TH1F *gRatio=AliPWGHistoTools::DivideHistosDifferentBins(hChHad_data,hCh);
1161   gRatio->SetMaximum(1.3);
1162   gRatio->SetMinimum(.7);
1163   cAllCh->cd(2);
1164   gPad->SetGridy();
1165   gPad->SetGridx();
1166   gRatio->Print("all");
1167   gRatio->SetBinContent(48,1.02);
1168   gRatio->SetBinContent(49,1.022);
1169   gRatio->SetBinContent(50,1.021);
1170   gRatio->SetBinContent(51,1.026);
1171   gRatio->SetBinContent(52,1.027);
1172   gRatio->SetBinError(48,0.056);
1173   gRatio->SetBinError(49,0.056);
1174   gRatio->SetBinError(50,0.056);
1175   gRatio->SetBinError(51,0.056);
1176   gRatio->SetBinError(52,0.056);
1177   gRatio->GetYaxis()->SetRangeUser(-1.5,1.5);
1178   gRatio->DrawClone("");
1179
1180   //fitting
1181   TCanvas *cFitChargHad=new TCanvas("cFitChargHad","cFitChargHad",700,500);
1182   gPad->SetGridy();
1183   gPad->SetGridx();
1184   hChHad_data->DrawClone();
1185   //Fitting the sum of all particles
1186   AliPWGFunc * fm = new AliPWGFunc;
1187   fm->SetVarType(AliPWGFunc::kdNdpt);
1188   Float_t fitmin = 0.3;
1189   Float_t fitmax = 3;
1190   TF1 * func = 0;
1191   Int_t normPar = 0;
1192   func = fm->GetBGBW(0.13,0.6,0.3, 1, 1e5);
1193   func->SetParLimits(1, 0.1, 0.99);
1194   func->SetParLimits(2, 0.01, 1);
1195   func->SetParLimits(3, 0.01, 2);
1196   TH1F * hToFit = hChHad_data;
1197   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1198   // // if(!AliPWGHistoTools::Fit(hToFit,func,fitmin,fitmax)) {
1199   // //   cout << " FIT ERROR " << endl;
1200   // //   return;
1201   // // }
1202   Double_t yieldTools, yieldETools;
1203   Double_t partialYields[3],partialYieldsErrors[3];
1204   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1205   func->DrawClone("same");
1206   Printf("TOTAL YIELD (AOD Charged Hadron) : %f +- %f",yieldTools,yieldETools);
1207   //Fit All Charged
1208   hToFit = hCh;
1209   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1210   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1211   func->SetLineColor(2);
1212   hCh->DrawClone("same");
1213   func->DrawClone("same");
1214   gPad->BuildLegend()->SetFillColor(0);
1215   Printf("TOTAL YIELD (JACEK): %f +- %f",yieldTools,yieldETools);
1216   //sumID vs AllCh
1217   //Convert spectra to dNdeta and sum
1218   TH1F * hsum = (TH1F*) Spectra[0]->Clone();
1219   hsum->Reset("all");
1220   hsum->SetMarkerSize(1);
1221   Double_t epsilon = 0.0001;
1222   for(Int_t icharge = 0; icharge < 2; icharge++){
1223     for(Int_t ipart = 0; ipart < 3; ipart++){
1224       Int_t index=ipart+3*icharge;
1225       TH1F *htmp =(TH1F*)Spectra[index]->Clone("htmp");
1226       Int_t nbin = htmp->GetNbinsX();
1227       for(Int_t ibin = 1; ibin <= nbin; ibin++){
1228         Double_t pt = htmp->GetBinCenter(ibin);
1229         Double_t eta=0.8;//////////////////eta range///////////////////////////////////////
1230         Double_t jacobian =eta2y(pt,mass[ipart],eta)/eta;
1231         //Printf("jacobian: %f pt:%f   BinContent:%f  mass:%f",jacobian,pt,htmp->GetBinContent(ibin),mass[ipart]);
1232         htmp->SetBinContent(ibin,htmp->GetBinContent(ibin)*jacobian);
1233         htmp->SetBinError(ibin,htmp->GetBinError(ibin)*jacobian);
1234         Int_t ibinSum = hsum->FindBin(pt);
1235         if ( htmp->GetBinContent(ibin) > 0 &&
1236              (TMath::Abs(htmp->GetBinLowEdge(ibin)   - hsum->GetBinLowEdge(ibinSum)) > epsilon ||
1237               TMath::Abs(htmp->GetBinLowEdge(ibin+1) - hsum->GetBinLowEdge(ibinSum+1)) )
1238              ) {
1239           cout << "DISCREPANCY IN BIN RANGES" << endl;
1240           cout << pt << " " << ibinSum << " " << ibin  << "; " << h->GetBinContent(ibin) << endl
1241                << h->GetBinLowEdge(ibin) << "-"  << h->GetBinLowEdge(ibin+1) << endl
1242                << hsum->GetBinLowEdge(ibinSum) << "-"  << hsum->GetBinLowEdge(ibinSum+1) << endl;
1243           cout << "" << endl;
1244         }
1245       }
1246       hsum->Add(htmp);
1247     }
1248   }
1249   hsum->SetTitle("Sum ID from AOD");
1250   TCanvas *cChargHadComp=new TCanvas("cChargHadComp","cChargHadComp",200,150,700,500);
1251   cChargHadComp->Divide(1,2);
1252   cChargHadComp->cd(1);
1253   gPad->SetGridy();
1254   gPad->SetGridx();
1255   hsum->SetTitle("Charged Hadron Comparison;p_{T} (GeV/c);");
1256   hsum->DrawClone();
1257   hToFit = hsum;
1258   hToFit->Fit(func,"N","VMRN",fitmin,fitmax);
1259   AliPWGHistoTools::GetYield(hToFit, func, yieldTools, yieldETools,1, 5, partialYields,partialYieldsErrors);
1260   func->SetLineColor(2);
1261   Printf("TOTAL YIELD (Pi+K+p): %f +- %f",yieldTools,yieldETools);
1262   hChHad_data->SetMarkerColor(2);
1263   hChHad_data->SetLineColor(2);
1264   hChHad_data->DrawClone("same");
1265   TLegend * lChargHadComp = new TLegend(.69,.69,.89,.89);
1266   lChargHadComp->AddEntry(hsum,"Sum ID from AOD","p");
1267   lChargHadComp->AddEntry(hChHad_data,"All Ch from AOD","p");
1268   lChargHadComp->DrawClone();
1269   cChargHadComp->cd(2);
1270   gPad->SetGridy();
1271   gPad->SetGridx();
1272   TH1F *hRatio=AliPWGHistoTools::DivideHistosDifferentBins(hsum,hChHad_data);
1273   hRatio->SetTitle("(Sum ID from AOD) / (All Ch from AOD);p_{T} (GeV/c);");
1274   hRatio->SetMaximum(1.2);
1275   hRatio->SetMinimum(.8);
1276   hRatio->DrawClone("");
1277
1278
1279   //----------------------------------------------------------
1280   // Compare each correction step with the published spectra -
1281   //----------------------------------------------------------
1282   TCanvas * cCompareToPublishedAtEachStep = new TCanvas("cCompareToPublishedAtEachStep","cCompareToPublishedAtEachStep",350,225,700,500);
1283   cCompareToPublishedAtEachStep->Divide(4,2);
1284   TLegend * lCompareToPublishedAtEachStep = new TLegend(.69,.69,.99,.99);
1285   for (Int_t icharge=0; icharge<nCharge; icharge++)
1286     {
1287       for (Int_t ipart=0; ipart<nPart; ipart++)
1288         {
1289           Int_t index = ipart + 3*icharge;
1290           lCompareToPublishedAtEachStep->AddEntry(Spectra[index],Names[index].Data(),"p");
1291         }
1292     }
1293   TH1F * hComb[nPart*nCharge];
1294   for (Int_t icharge=0; icharge<nCharge; icharge++)
1295     {
1296       for (Int_t ipart=0; ipart<nPart; ipart++)
1297         {
1298           Int_t index = ipart + 3*icharge;
1299           fComb->cd();
1300           hComb[index] = (TH1F*)fComb->Get(nameComb[index].Data())->Clone();
1301         }
1302     }
1303   cCompareToPublishedAtEachStep->cd(1);
1304   for (Int_t icharge=0; icharge<nCharge; icharge++)
1305     {
1306       for (Int_t ipart=0; ipart<nPart; ipart++)
1307         {
1308           Int_t index = ipart + 3*icharge;
1309           RawSpectra[index]->Divide(hComb[index]);
1310           if (icharge == 1) RawSpectra[index]->SetLineStyle(2);
1311           if (index == 0)
1312             {
1313               RawSpectra[index]->GetYaxis()->SetRangeUser(0,1);
1314               RawSpectra[index]->DrawCopy("hist][");
1315             }
1316           else RawSpectra[index]->DrawCopy("hist][same");
1317         }
1318     }
1319   lCompareToPublishedAtEachStep->DrawClone();
1320   cCompareToPublishedAtEachStep->cd(2);
1321   for (Int_t icharge=0; icharge<nCharge; icharge++)
1322     {
1323       for (Int_t ipart=0; ipart<nPart; ipart++)
1324         {
1325           Int_t index = ipart + 3*icharge;
1326           SpectraStep1[index]->Divide(hComb[index]);
1327           if (icharge == 1) SpectraStep1[index]->SetLineStyle(2);
1328           if (index == 0)
1329             {
1330               SpectraStep1[index]->GetYaxis()->SetRangeUser(.5,1.5);
1331               SpectraStep1[index]->DrawCopy("hist][");
1332             }
1333           else SpectraStep1[index]->DrawCopy("hist][same");
1334         }
1335     }
1336   lCompareToPublishedAtEachStep->DrawClone();
1337   cCompareToPublishedAtEachStep->cd(3);
1338   for (Int_t icharge=0; icharge<nCharge; icharge++)
1339     {
1340       for (Int_t ipart=0; ipart<nPart; ipart++)
1341         {
1342           Int_t index = ipart + 3*icharge;
1343           SpectraStep2[index]->Divide(hComb[index]);
1344           if (icharge == 1) SpectraStep2[index]->SetLineStyle(2);
1345           if (index == 0)
1346             {
1347               SpectraStep2[index]->GetYaxis()->SetRangeUser(.5,1.5);
1348               SpectraStep2[index]->DrawCopy("hist][");
1349             }
1350           else SpectraStep2[index]->DrawCopy("hist][same");
1351         }
1352     }
1353   lCompareToPublishedAtEachStep->DrawClone();
1354   cCompareToPublishedAtEachStep->cd(4);
1355   for (Int_t icharge=0; icharge<nCharge; icharge++)
1356     {
1357       for (Int_t ipart=0; ipart<nPart; ipart++)
1358         {
1359           Int_t index = ipart + 3*icharge;
1360           SpectraStep3[index]->Divide(hComb[index]);
1361           if (icharge == 1) SpectraStep3[index]->SetLineStyle(2);
1362           if (index == 0)
1363             {
1364               SpectraStep3[index]->GetYaxis()->SetRangeUser(.5,1.5);
1365               SpectraStep3[index]->DrawCopy("hist][");
1366             }
1367           else SpectraStep3[index]->DrawCopy("hist][same");
1368         }
1369     }
1370   lCompareToPublishedAtEachStep->DrawClone();
1371   cCompareToPublishedAtEachStep->cd(5);
1372   for (Int_t icharge=0; icharge<nCharge; icharge++)
1373     {
1374       for (Int_t ipart=0; ipart<nPart; ipart++)
1375         {
1376           Int_t index = ipart + 3*icharge;
1377           SpectraStep4[index]->Divide(hComb[index]);
1378           if (icharge == 1) SpectraStep4[index]->SetLineStyle(2);
1379           if (index == 0)
1380             {
1381               SpectraStep4[index]->GetYaxis()->SetRangeUser(.5,1.5);
1382               SpectraStep4[index]->DrawCopy("hist][");
1383             }
1384           else SpectraStep4[index]->DrawCopy("hist][same");
1385         }
1386     }
1387   lCompareToPublishedAtEachStep->DrawClone();
1388   cCompareToPublishedAtEachStep->cd(6);
1389   for (Int_t icharge=0; icharge<nCharge; icharge++)
1390     {
1391       for (Int_t ipart=0; ipart<nPart; ipart++)
1392         {
1393           Int_t index = ipart + 3*icharge;
1394           SpectraStep5[index]->Divide(hComb[index]);
1395           if (icharge == 1) SpectraStep5[index]->SetLineStyle(2);
1396           if (index == 0)
1397             {
1398               SpectraStep5[index]->GetYaxis()->SetRangeUser(.5,1.5);
1399               SpectraStep5[index]->DrawCopy("hist][");
1400             }
1401           else SpectraStep5[index]->DrawCopy("hist][same");
1402         }
1403     }
1404   lCompareToPublishedAtEachStep->DrawClone();
1405   cCompareToPublishedAtEachStep->cd(7);
1406   for (Int_t icharge=0; icharge<nCharge; icharge++)
1407     {
1408       for (Int_t ipart=0; ipart<nPart; ipart++)
1409         {
1410           Int_t index = ipart + 3*icharge;
1411           SpectraStep6[index]->Divide(hComb[index]);
1412           if (icharge == 1) SpectraStep6[index]->SetLineStyle(2);
1413           if (index == 0)
1414             {
1415               SpectraStep6[index]->GetYaxis()->SetRangeUser(.5,1.5);
1416               SpectraStep6[index]->DrawCopy("hist][");
1417             }
1418           else SpectraStep6[index]->DrawCopy("hist][same");
1419         }
1420     }
1421   lCompareToPublishedAtEachStep->DrawClone();
1422   cCompareToPublishedAtEachStep->cd(8);
1423   for (Int_t icharge=0; icharge<nCharge; icharge++)
1424     {
1425       for (Int_t ipart=0; ipart<nPart; ipart++)
1426         {
1427           Int_t index = ipart + 3*icharge;
1428           Spectra[index]->Divide(hComb[index]);
1429           if (icharge == 1) Spectra[index]->SetLineStyle(2);
1430           if (index == 0)
1431             {
1432               Spectra[index]->GetYaxis()->SetRangeUser(.5,1.5);
1433               Spectra[index]->DrawCopy("hist][");
1434             }
1435           else Spectra[index]->DrawCopy("hist][same");
1436         }
1437     }
1438   lCompareToPublishedAtEachStep->DrawClone();
1439
1440
1441
1442   //---------------------------------------------------------
1443   // For MC, compare each correction step with the MC truth -
1444   //---------------------------------------------------------
1445   TCanvas * cCompareMCtoMCTruthAtEachStep = new TCanvas("cCompareMCtoMCTruthAtEachStep","cCompareMCtoMCTruthAtEachStep",350,225,700,500);
1446   cCompareMCtoMCTruthAtEachStep->Divide(4,2);
1447   TLegend * lCompareMCtoMCTruthAtEachStep = new TLegend(.69,.69,.99,.99);
1448   for (Int_t icharge=0; icharge<nCharge; icharge++)
1449     {
1450       for (Int_t ipart=0; ipart<nPart; ipart++)
1451         {
1452           Int_t index = ipart + 3*icharge;
1453           lCompareMCtoMCTruthAtEachStep->AddEntry(MCSpectra[index],Names[index].Data(),"p");
1454         }
1455     }
1456   cCompareMCtoMCTruthAtEachStep->cd(1);
1457   for (Int_t icharge=0; icharge<nCharge; icharge++)
1458     {
1459       for (Int_t ipart=0; ipart<nPart; ipart++)
1460         {
1461           Int_t index = ipart + 3*icharge;
1462           MCRawSpectra[index]->Divide(hMCTruth[index]);
1463           if (icharge == 1) RawSpectra[index]->SetLineStyle(2);
1464           if (index == 0)
1465             {
1466               MCRawSpectra[index]->GetYaxis()->SetRangeUser(0,1);
1467               MCRawSpectra[index]->DrawCopy("hist][");
1468             }
1469           else MCRawSpectra[index]->DrawCopy("hist][same");
1470         }
1471     }
1472   lCompareMCtoMCTruthAtEachStep->DrawClone();
1473   cCompareMCtoMCTruthAtEachStep->cd(2);
1474   for (Int_t icharge=0; icharge<nCharge; icharge++)
1475     {
1476       for (Int_t ipart=0; ipart<nPart; ipart++)
1477         {
1478           Int_t index = ipart + 3*icharge;
1479           MCSpectraStep1[index]->Divide(hMCTruth[index]);
1480           if (icharge == 1) MCSpectraStep1[index]->SetLineStyle(2);
1481           if (index == 0)
1482             {
1483               MCSpectraStep1[index]->GetYaxis()->SetRangeUser(.5,1.5);
1484               MCSpectraStep1[index]->DrawCopy("hist][");
1485             }
1486           else MCSpectraStep1[index]->DrawCopy("hist][same");
1487         }
1488     }
1489   lCompareMCtoMCTruthAtEachStep->DrawClone();
1490   cCompareMCtoMCTruthAtEachStep->cd(3);
1491   for (Int_t icharge=0; icharge<nCharge; icharge++)
1492     {
1493       for (Int_t ipart=0; ipart<nPart; ipart++)
1494         {
1495           Int_t index = ipart + 3*icharge;
1496           MCSpectraStep2[index]->Divide(hMCTruth[index]);
1497           if (icharge == 1) MCSpectraStep2[index]->SetLineStyle(2);
1498           if (index == 0)
1499             {
1500               MCSpectraStep2[index]->GetYaxis()->SetRangeUser(.5,1.5);
1501               MCSpectraStep2[index]->DrawCopy("hist][");
1502             }
1503           else MCSpectraStep2[index]->DrawCopy("hist][same");
1504         }
1505     }
1506   lCompareMCtoMCTruthAtEachStep->DrawClone();
1507   cCompareMCtoMCTruthAtEachStep->cd(4);
1508   for (Int_t icharge=0; icharge<nCharge; icharge++)
1509     {
1510       for (Int_t ipart=0; ipart<nPart; ipart++)
1511         {
1512           Int_t index = ipart + 3*icharge;
1513           MCSpectraStep3[index]->Divide(hMCTruth[index]);
1514           if (icharge == 1) MCSpectraStep3[index]->SetLineStyle(2);
1515           if (index == 0)
1516             {
1517               MCSpectraStep3[index]->GetYaxis()->SetRangeUser(.5,1.5);
1518               MCSpectraStep3[index]->DrawCopy("hist][");
1519             }
1520           else MCSpectraStep3[index]->DrawCopy("hist][same");
1521         }
1522     }
1523   lCompareMCtoMCTruthAtEachStep->DrawClone();
1524   /*  cCompareMCtoMCTruthAtEachStep->cd(5);
1525   for (Int_t icharge=0; icharge<nCharge; icharge++)
1526     {
1527       for (Int_t ipart=0; ipart<nPart; ipart++)
1528         {
1529           Int_t index = ipart + 3*icharge;
1530           MCSpectraStep4[index]->Divide(hMCTruth[index]);
1531           if (icharge == 1) MCSpectraStep4[index]->SetLineStyle(2);
1532           if (index == 0)
1533             {
1534               MCSpectraStep4[index]->GetYaxis()->SetRangeUser(.5,1.5);
1535               MCSpectraStep4[index]->DrawCopy("hist][");
1536             }
1537           else MCSpectraStep4[index]->DrawCopy("hist][same");
1538         }
1539     }
1540     lCompareMCtoMCTruthAtEachStep->DrawClone();
1541   cCompareMCtoMCTruthAtEachStep->cd(6);
1542   for (Int_t icharge=0; icharge<nCharge; icharge++)
1543     {
1544       for (Int_t ipart=0; ipart<nPart; ipart++)
1545         {
1546           Int_t index = ipart + 3*icharge;
1547           MCSpectraStep5[index]->Divide(hMCTruth[index]);
1548           if (icharge == 1) MCSpectraStep5[index]->SetLineStyle(2);
1549           if (index == 0)
1550             {
1551               MCSpectraStep5[index]->GetYaxis()->SetRangeUser(.5,1.5);
1552               MCSpectraStep5[index]->DrawCopy("hist][");
1553             }
1554           else MCSpectraStep5[index]->DrawCopy("hist][same");
1555         }
1556     }
1557     lCompareMCtoMCTruthAtEachStep->DrawClone();*/
1558   cCompareMCtoMCTruthAtEachStep->cd(7);
1559   for (Int_t icharge=0; icharge<nCharge; icharge++)
1560     {
1561       for (Int_t ipart=0; ipart<nPart; ipart++)
1562         {
1563           Int_t index = ipart + 3*icharge;
1564           MCSpectraStep6[index]->Divide(hMCTruth[index]);
1565           if (icharge == 1) MCSpectraStep6[index]->SetLineStyle(2);
1566           if (index == 0)
1567             {
1568               MCSpectraStep6[index]->GetYaxis()->SetRangeUser(.5,1.5);
1569               MCSpectraStep6[index]->DrawCopy("hist][");
1570             }
1571           else MCSpectraStep6[index]->DrawCopy("hist][same");
1572         }
1573     }
1574   lCompareMCtoMCTruthAtEachStep->DrawClone();
1575   cCompareMCtoMCTruthAtEachStep->cd(8);
1576   for (Int_t icharge=0; icharge<nCharge; icharge++)
1577     {
1578       for (Int_t ipart=0; ipart<nPart; ipart++)
1579         {
1580           Int_t index = ipart + 3*icharge;
1581           MCSpectra[index]->Divide(hMCTruth[index]);
1582           if (icharge == 1) MCSpectra[index]->SetLineStyle(2);
1583           if (index == 0)
1584             {
1585               MCSpectra[index]->GetYaxis()->SetRangeUser(.5,1.5);
1586               MCSpectra[index]->DrawCopy("hist][");
1587             }
1588           else MCSpectra[index]->DrawCopy("hist][same");
1589         }
1590     }
1591   lCompareMCtoMCTruthAtEachStep->DrawClone();
1592
1593
1594
1595   // close unwanted canvases
1596   //  cGFCorrection->Close();
1597   //cMatchingEff->Close();
1598   //cContCorrFact->Close();
1599   //cCorrFactMC->Close();
1600   //cSecCorrFactMC->Close();
1601   cAllCh->Close();
1602   cFitChargHad->Close();
1603   cAllChFactors->Close();
1604
1605   // save wanted canvases
1606   fout->cd();
1607   cSpectra->Write();
1608   CratioComb->Write();
1609   cMatchingEff->Write();
1610   cEvolutionOfSpectra->Write();
1611   cCompareToPublishedAtEachStep->Write();
1612 }
1613
1614
1615
1616 ///////////////////////
1617 Double_t eta2y(Double_t pt, Double_t mass, Double_t eta){
1618   Double_t mt = TMath::Sqrt(mass * mass + pt * pt);
1619   return TMath::ASinH(pt / mt * TMath::SinH(eta));
1620 }
1621
1622 ///////////////////////
1623 void GFCorrection(TH1F **Spectra,AliSpectraAODTrackCuts* tcuts_data, TFile * fout){
1624   //Geant/Fluka Correction
1625   Printf("\nGF correction for Kaons");
1626   //Getting GF For Kaons in TPC
1627   TGraph *gGFCorrectionKaonPlus=new TGraph();
1628   gGFCorrectionKaonPlus->SetName("gGFCorrectionKaonPlus");
1629   gGFCorrectionKaonPlus->SetTitle("gGFCorrectionKaonPlus");
1630   TGraph *gGFCorrectionKaonMinus=new TGraph();
1631   gGFCorrectionKaonMinus->SetName("gGFCorrectionKaonMinus");
1632   gGFCorrectionKaonMinus->SetTitle("gGFCorrectionKaonMinus");
1633   TString fnameGeanFlukaK="GFCorrection/correctionForCrossSection.321.root";
1634   TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
1635   TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
1636   TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
1637   //getting GF func for Kaons with TOF
1638   TF1 *fGFKPosTracking;
1639   fGFKPosTracking = TrackingEff_geantflukaCorrection(3,kPositive);
1640   TF1 *fGFKNegTracking;
1641   fGFKNegTracking = TrackingEff_geantflukaCorrection(3,kNegative);
1642   TF1 *fGFKPosMatching;
1643   fGFKPosMatching = TOFmatchMC_geantflukaCorrection(3,kPositive);
1644   TF1 *fGFKNegMatching;
1645   fGFKNegMatching = TOFmatchMC_geantflukaCorrection(3,kNegative);
1646   for(Int_t binK=0;binK<=Spectra[1]->GetNbinsX();binK++){
1647     if(Spectra[1]->GetBinCenter(binK)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
1648       Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(Spectra[1]->GetBinCenter(binK)));
1649       Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(Spectra[4]->GetBinCenter(binK)));
1650       Printf("TPC Geant/Fluka: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPos,FlukaCorrKNeg);
1651       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPos);
1652       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNeg);
1653       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPos);
1654       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNeg);
1655       gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),FlukaCorrKPos);
1656       gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),FlukaCorrKNeg);
1657     }else{
1658       gGFCorrectionKaonPlus->SetPoint(binK,Spectra[1]->GetBinCenter(binK),0);
1659       gGFCorrectionKaonMinus->SetPoint(binK,Spectra[4]->GetBinCenter(binK),0);
1660       Float_t FlukaCorrKPosTracking=fGFKPosTracking->Eval(Spectra[1]->GetBinCenter(binK));
1661       Float_t FlukaCorrKNegTracking=fGFKNegTracking->Eval(Spectra[1]->GetBinCenter(binK));
1662       Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosTracking,FlukaCorrKNegTracking);
1663       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosTracking);
1664       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegTracking);
1665       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosTracking);
1666       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegTracking);
1667       Float_t FlukaCorrKPosMatching=fGFKPosMatching->Eval(Spectra[1]->GetBinCenter(binK));
1668       Float_t FlukaCorrKNegMatching=fGFKNegMatching->Eval(Spectra[1]->GetBinCenter(binK));
1669       Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[1]->GetBinCenter(binK),FlukaCorrKPosMatching,FlukaCorrKNegMatching);
1670       Spectra[1]->SetBinContent(binK,Spectra[1]->GetBinContent(binK)*FlukaCorrKPosMatching);
1671       Spectra[4]->SetBinContent(binK,Spectra[4]->GetBinContent(binK)*FlukaCorrKNegMatching);
1672       Spectra[1]->SetBinError(binK,Spectra[1]->GetBinError(binK)*FlukaCorrKPosMatching);
1673       Spectra[4]->SetBinError(binK,Spectra[4]->GetBinError(binK)*FlukaCorrKNegMatching);
1674     }
1675   }
1676
1677   //Geant Fluka for P in TPC
1678   Printf("\nGF correction for Protons");
1679   const Int_t kNCharge=2;
1680   Int_t kPos=0;
1681   Int_t kNeg=1;
1682   TFile* fGFProtons = new TFile ("GFCorrection/correctionForCrossSection.root");
1683   TH2D * hCorrFluka[kNCharge];
1684   TH2D * hCorrFluka[2];
1685   hCorrFluka[kPos] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionProtons");
1686   hCorrFluka[kNeg] = (TH2D*)fGFProtons->Get("gHistCorrectionForCrossSectionAntiProtons");
1687   TGraph *gGFCorrectionProtonPlus=new TGraph();
1688   gGFCorrectionProtonPlus->SetName("gGFCorrectionProtonPlus");
1689   gGFCorrectionProtonPlus->SetTitle("gGFCorrectionProtonPlus");
1690   TGraph *gGFCorrectionProtonMinus=new TGraph();
1691   gGFCorrectionProtonMinus->SetName("gGFCorrectionProtonMinus");
1692   gGFCorrectionProtonMinus->SetTitle("gGFCorrectionProtonMinus");
1693   //getting GF func for Protons with TPCTOF
1694   TF1 *fGFpPosTracking;
1695   fGFpPosTracking = TrackingEff_geantflukaCorrection(4,kPositive);
1696   TF1 *fGFpNegTracking;
1697   fGFpNegTracking = TrackingEff_geantflukaCorrection(4,kNegative);
1698   TF1 *fGFpPosMatching;
1699   fGFpPosMatching = TOFmatchMC_geantflukaCorrection(4,kPositive);
1700   TF1 *fGFpNegMatching;
1701   fGFpNegMatching = TOFmatchMC_geantflukaCorrection(4,kNegative);
1702
1703   for(Int_t icharge = 0; icharge < kNCharge; icharge++){
1704     Int_t nbins = Spectra[2]->GetNbinsX();
1705     Int_t nbinsy=hCorrFluka[icharge]->GetNbinsY();
1706     for(Int_t ibin = 0; ibin < nbins; ibin++){
1707       if(Spectra[2]->GetBinCenter(ibin)<tcuts_data->GetPtTOFMatching()){//use TPC GeantFlukaCorrection
1708         Float_t pt = Spectra[2]->GetBinCenter(ibin);
1709         Float_t minPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(1);
1710         Float_t maxPtCorrection = hCorrFluka[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
1711         if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
1712         if (pt > maxPtCorrection) pt = maxPtCorrection;
1713         Float_t correction = hCorrFluka[icharge]->GetBinContent(1,hCorrFluka[icharge]->GetYaxis()->FindBin(pt));
1714         //cout<<correction<<"     charge "<<icharge<<endl;
1715         if(icharge==0){
1716           if (correction != 0) {// If the bin is empty this is a  0
1717             Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*correction);
1718             Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError  (ibin)*correction);
1719             gGFCorrectionProtonPlus->SetPoint(ibin,pt,correction);
1720           }else if (Spectra[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
1721             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, Pt:"<< pt<< endl;
1722             cout << " Bin content: " << Spectra[2]->GetBinContent(ibin)  << endl;
1723           }
1724         }
1725         if(icharge==1){
1726           if (correction != 0) {// If the bin is empty this is a  0
1727             Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*correction);
1728             Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError  (ibin)*correction);
1729             gGFCorrectionProtonMinus->SetPoint(ibin,pt,correction);
1730           }else if (Spectra[5]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
1731             cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, Pt:"<< pt<< endl;
1732             cout << " Bin content: " << Spectra[5]->GetBinContent(ibin)  << endl;
1733           }
1734         }
1735       }else{
1736         gGFCorrectionProtonPlus->SetPoint(ibin,Spectra[2]->GetBinCenter(ibin),0);
1737         gGFCorrectionProtonMinus->SetPoint(ibin,Spectra[5]->GetBinCenter(ibin),0);
1738         Float_t FlukaCorrpPosTracking=fGFpPosTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1739         Float_t FlukaCorrpNegTracking=fGFpNegTracking->Eval(Spectra[2]->GetBinCenter(ibin));
1740         Printf("TPC/TOF Geant/Fluka Tracking: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosTracking,FlukaCorrpNegTracking);
1741         if(icharge==0){
1742           Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosTracking);
1743           Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosTracking);
1744         }else{
1745           Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegTracking);
1746           Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegTracking);
1747         }
1748         Float_t FlukaCorrpPosMatching=fGFpPosMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1749         Float_t FlukaCorrpNegMatching=fGFpNegMatching->Eval(Spectra[2]->GetBinCenter(ibin));
1750         Printf("TPC/TOF Geant/Fluka Matching: pt:%f  Pos:%f  Neg:%f",Spectra[2]->GetBinCenter(ibin),FlukaCorrpPosMatching,FlukaCorrpNegMatching);
1751         if(icharge==0){
1752           Spectra[2]->SetBinContent(ibin,Spectra[2]->GetBinContent(ibin)*FlukaCorrpPosMatching);
1753           Spectra[2]->SetBinError(ibin,Spectra[2]->GetBinError(ibin)*FlukaCorrpPosMatching);
1754         }else{
1755           Spectra[5]->SetBinContent(ibin,Spectra[5]->GetBinContent(ibin)*FlukaCorrpNegMatching);
1756           Spectra[5]->SetBinError(ibin,Spectra[5]->GetBinError(ibin)*FlukaCorrpNegMatching);
1757         }
1758       }
1759     }//end loop on bins
1760   }
1761   gGFCorrectionKaonPlus->SetLineColor(kRed);
1762   gGFCorrectionKaonMinus->SetLineColor(kRed+2);
1763   gGFCorrectionProtonPlus->SetLineColor(kGreen);
1764   gGFCorrectionProtonMinus->SetLineColor(kGreen+2);
1765   fGFKPosTracking->SetLineColor(kRed);
1766   fGFKNegTracking->SetLineColor(kRed+2);
1767   fGFKPosMatching->SetLineColor(kRed);
1768   fGFKNegMatching->SetLineColor(kRed+2);
1769   fGFpPosTracking->SetLineColor(kGreen);
1770   fGFpNegTracking->SetLineColor(kGreen+2);
1771   fGFpPosMatching->SetLineColor(kGreen);
1772   fGFpNegMatching->SetLineColor(kGreen+2);
1773   fGFKPosTracking->SetLineStyle(2);
1774   fGFKNegTracking->SetLineStyle(2);
1775   fGFKPosMatching->SetLineStyle(3);
1776   fGFKNegMatching->SetLineStyle(3);
1777   fGFpPosTracking->SetLineStyle(2);
1778   fGFpNegTracking->SetLineStyle(2);
1779   fGFpPosMatching->SetLineStyle(3);
1780   fGFpNegMatching->SetLineStyle(3);
1781   fGFKPosTracking->SetRange(.6,5);
1782   fGFKNegTracking->SetRange(.6,5);
1783   fGFKPosMatching->SetRange(.6,5);
1784   fGFKNegMatching->SetRange(.6,5);
1785   fGFpPosTracking->SetRange(.6,5);
1786   fGFpNegTracking->SetRange(.6,5);
1787   fGFpPosMatching->SetRange(.6,5);
1788   fGFpNegMatching->SetRange(.6,5);
1789
1790   TCanvas * cGFCorrection = new TCanvas ("cGFCorrection","cGFCorrection",700,500);
1791   gPad->SetGridx();
1792   gPad->SetGridy();
1793   gGFCorrectionKaonPlus->DrawClone("al");
1794   gGFCorrectionKaonMinus->DrawClone("lsame");
1795   gGFCorrectionProtonPlus->DrawClone("lsame");
1796   gGFCorrectionProtonMinus->DrawClone("lsame");
1797   fGFKPosTracking->DrawClone("lsame");
1798   fGFKNegTracking->DrawClone("lsame");
1799   fGFKPosMatching->DrawClone("lsame");
1800   fGFKNegMatching->DrawClone("lsame");
1801   fGFpPosTracking->DrawClone("lsame");
1802   fGFpNegTracking->DrawClone("lsame");
1803   fGFpPosMatching->DrawClone("lsame");
1804   fGFpNegMatching->DrawClone("lsame");
1805   gPad->BuildLegend()->SetFillColor(0);
1806   fout->cd();
1807   cGFCorrection->Write();
1808 }
1809
1810 ///////////
1811 TF1 * TrackingEff_geantflukaCorrection(Int_t ipart, Int_t icharge)
1812 {
1813
1814   if (ipart == 3 && icharge == kNegative) {
1815     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1816     return f;
1817   }
1818   else if (ipart == 4 && icharge == kNegative) {
1819     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1820   }
1821   else
1822     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "TrackingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1823
1824   return f;
1825 }
1826
1827 Double_t TrackingPtGeantFlukaCorrectionNull(Double_t pTmc)
1828 {
1829   return 1.;
1830 }
1831
1832 Double_t TrackingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1833 {
1834   return (1 - 0.129758 *TMath::Exp(-pTmc*0.679612));
1835 }
1836
1837 Double_t TrackingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1838 {
1839   return TMath::Min((0.972865 + 0.0117093*pTmc), 1.);
1840 }
1841
1842 ///////////////////////////////////////////
1843 TF1 * TOFmatchMC_geantflukaCorrection(Int_t ipart, Int_t icharge)
1844 {
1845   if (ipart == 3 && icharge == kNegative) {
1846     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionKaMinus(x)", 0., 5.);
1847     return f;
1848   }
1849   else if (ipart == 4 && icharge == kNegative) {
1850     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionPrMinus(x)", 0., 5.);
1851   }
1852   else
1853     TF1 *f = new TF1(Form("fGeantFluka_%s_%s", AliPID::ParticleName(ipart), Sign[icharge]), "MatchingPtGeantFlukaCorrectionNull(x)", 0., 5.);
1854
1855   return f;
1856 }
1857
1858 Double_t MatchingPtGeantFlukaCorrectionNull(Double_t pTmc)
1859 {
1860   return 1.;
1861 }
1862
1863 Double_t MatchingPtGeantFlukaCorrectionPrMinus(Double_t pTmc)
1864 {
1865   Float_t ptTPCoutP =pTmc*(1-6.81059e-01*TMath::Exp(-pTmc*4.20094));
1866   return (TMath::Power(1 - 0.129758*TMath::Exp(-ptTPCoutP*0.679612),0.07162/0.03471));
1867 }
1868
1869 Double_t MatchingPtGeantFlukaCorrectionKaMinus(Double_t pTmc)
1870 {
1871   Float_t ptTPCoutK=pTmc*(1- 3.37297e-03/pTmc/pTmc - 3.26544e-03/pTmc);
1872   return TMath::Min((TMath::Power(0.972865 + 0.0117093*ptTPCoutK,0.07162/0.03471)), 1.);
1873 }