--- /dev/null
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TF1.h>
+#include <TPad.h>
+#include <TGraphErrors.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TGrid.h>
+#include <TGridResult.h>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TLatex.h>
+#include <TImage.h>
+#include <TPaveText.h>
+#include <TLine.h>
+#endif
+
+void RatioPlot(TH1F **num,TH1F **den,TString t1,TString t2,TString opt,Double_t Low[],Double_t Up[]);
+void SetDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1);
+void ResetOutOfRange(TH1F *histo, Int_t ipart, Double_t lowRange[3], Double_t upRange[3]);
+
+void MakeCorrectedSpectraITSsaNsigma( TString period="LHC10d",
+ TString MCname="LHC10f6a",
+ Float_t DCAcut=7,
+ Int_t lowmult=-1,
+ Int_t upmult=-1
+ )
+{
+
+ // parameters
+
+ Bool_t MakeMCPlots=1;
+ Bool_t MakeSystErr=0;
+ Bool_t NormalizeData=1;
+ Bool_t NormalizeToINEL=0;
+ Bool_t SaveOnlyAsymm=0;
+ Bool_t CorrectOnlyForEfficiency=0;
+ Bool_t MarekNorm=0;
+ Bool_t DCACorrection=1;
+ Bool_t GFCorrection=1;
+
+ //MC Histos
+ TH1F *fHistPrimMCposBefEvSel[3];
+ TH1F *fHistPrimMCnegBefEvSel[3];
+ TH1F *fHistPrimMCpos[3];
+ TH1F *fHistPrimMCneg[3];
+ TH1F *hHistPosNSigma[3];
+ TH1F *hHistNegNSigma[3];
+ TH1F *hHistPosNSigmaMean[3];
+ TH1F *hHistNegNSigmaMean[3];
+ TH1F *hHistPosNSigmaPrim[3];
+ TH1F *hHistNegNSigmaPrim[3];
+ TH1F *hHistPosNSigmaPrimMean[3];
+ TH1F *hHistNegNSigmaPrimMean[3];
+ TH1F *hHistPosNSigmaPrimMC[3];
+ TH1F *hHistNegNSigmaPrimMC[3];
+ TH1F *hHistPosNSigmaPrimMCMean[3];
+ TH1F *hHistNegNSigmaPrimMCMean[3];
+ TH1F *hHistPosNSigmaMC[3];
+ TH1F *hHistNegNSigmaMC[3];
+ TH1F *hHistPosNSigmaMCMean[3];
+ TH1F *hHistNegNSigmaMCMean[3];
+ TH1F *hCorrFactPos[3];
+ TH1F *hCorrFactNeg[3];
+ TH1F *hCorrFactPosMean[3];
+ TH1F *hCorrFactNegMean[3];
+ TH1F *hEffPos[3];
+ TH1F *hEffNeg[3];
+ TH1F *hEffPosMean[3];
+ TH1F *hEffNegMean[3];
+
+ //DATA Histos
+ TH1F *hHistPosNSigmaDATA[3];
+ TH1F *hHistNegNSigmaDATA[3];
+ TH1F *hHistPosNSigmaMeanDATA[3];
+ TH1F *hHistNegNSigmaMeanDATA[3];
+ TH1F *hITSsaPos[3];
+ TH1F *hITSsaNeg[3];
+ TH1F *hITSsaPosMean[3];
+ TH1F *hITSsaNegMean[3];
+ TH1F *hITSsaRawPos[3];
+ TH1F *hITSsaRawNeg[3];
+ TH1F *hITSsaRawPosMean[3];
+ TH1F *hITSsaRawNegMean[3];
+
+ // Float_t nmin=1.5;
+ Double_t Low[3]={0.10,0.2,0.3};
+ Double_t Up[3]={0.4,0.5,0.7};
+ Double_t lowRange[3]={0.10,0.2,0.3};
+ Double_t upRange[3]={0.6,0.6,0.6};
+
+
+
+ const Int_t nbinspt=22;
+ Double_t xbins[nbinspt+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
+
+ TString fname_MC=Form("./%s/AnalysisResults.root",MCname.Data());
+ cout<<fname_MC.Data()<<endl;
+
+ TString foutMC=Form("outSpectraMC_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
+ TFile *outMC=new TFile(foutMC.Data(),"recreate");
+ outMC->Close();
+ delete outMC;
+
+ TString foutDATA;
+ if(CorrectOnlyForEfficiency){
+ foutDATA=Form("outSpectraData_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
+ }else{
+ foutDATA=Form("outCorrSpectraData_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
+ }
+
+ TString fnameDATA=Form("%s/AnalysisResults.root",period.Data());
+ cout<<fnameDATA.Data()<<endl;
+ TFile *outDATA=new TFile(foutDATA.Data(),"recreate");
+ outDATA->Close();
+ delete outDATA;
+
+
+ TString lnameMC=Form("clistITSsaMult%ito%i",lowmult,upmult);
+ Printf("\n------------------ READING MC %s -------------------\n",lnameMC.Data());
+ TFile *finMC = new TFile(fname_MC.Data());
+ TDirectory *dirFileMC=(TDirectory*)finMC->Get("PWG2SpectraITSsa");
+ TList * cOutputMC = (TList*)dirFileMC->Get(lnameMC.Data());
+ //cOutputMC->Print();
+
+ TH1F *fHistNEventsMC = (TH1F*)cOutputMC->FindObject("fHistNEvents");
+ Float_t normMC=(Float_t)1/fHistNEventsMC->GetBinContent(fHistNEventsMC->FindBin(1.));
+ printf("- NormalizationMC (NEvents): %.0f \n",fHistNEventsMC->GetBinContent(fHistNEventsMC->FindBin(1.)));
+
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ fHistPrimMCposBefEvSel[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCposBefEvSel%i",ipart));//primaries generated (before event selection)
+ fHistPrimMCnegBefEvSel[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCnegBefEvSel%i",ipart));
+ fHistPrimMCpos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCpos%i",ipart)); //primaries generated (after event selection)
+ fHistPrimMCneg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCneg%i",ipart));
+ hHistPosNSigma[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigma%i",ipart));//NSigma recontructed, no MCtruth
+ hHistNegNSigma[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigma%i",ipart));
+ hHistPosNSigmaPrimMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMC%i",ipart));//NSigma recontructed, MCtruth
+ hHistNegNSigmaPrimMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMC%i",ipart));
+ hHistPosNSigmaMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMC%i",ipart));//NSigma recontructed, MCtruth
+ hHistNegNSigmaMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMC%i",ipart));
+ hHistPosNSigmaPrim[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrim%i",ipart));//NSigma recontructed, truth
+ hHistNegNSigmaPrim[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrim%i",ipart));
+ hHistPosNSigmaMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMean%i",ipart));//NSigmaMean recontructed, no MCtruth
+ hHistNegNSigmaMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMean%i",ipart));
+ hHistPosNSigmaPrimMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMCMean%i",ipart));//NSigmaMean recontructed, no MCtruth
+ hHistNegNSigmaPrimMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMCMean%i",ipart));
+ hHistPosNSigmaMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMCMean%i",ipart));//NSigmaMean recontructed, no MCtruth
+ hHistNegNSigmaMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMCMean%i",ipart));
+ hHistPosNSigmaPrimMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMean%i",ipart));//NSigmaMean recontructed, no MCtruth
+ hHistNegNSigmaPrimMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMean%i",ipart));
+ }
+
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ hCorrFactPos[ipart]=new TH1F(Form("hCorrFactPos%i",ipart),Form("hCorrFactPos%i",ipart),nbinspt,xbins);
+ hCorrFactNeg[ipart]=new TH1F(Form("hCorrFactNeg%i",ipart),Form("hCorrFactNeg%i",ipart),nbinspt,xbins);
+ hCorrFactPosMean[ipart]=new TH1F(Form("hCorrFactPosMean%i",ipart),Form("hCorrFactPosMean%i",ipart),nbinspt,xbins);
+ hCorrFactNegMean[ipart]=new TH1F(Form("hCorrFactNegMean%i",ipart),Form("hCorrFactNegMean%i",ipart),nbinspt,xbins);
+ hEffPos[ipart]=new TH1F(Form("hEffPos%i",ipart),Form("hEffPos%i",ipart),nbinspt,xbins);
+ hEffNeg[ipart]=new TH1F(Form("hEffNeg%i",ipart),Form("hEffNeg%i",ipart),nbinspt,xbins);
+ hEffPosMean[ipart]=new TH1F(Form("hEffPosMean%i",ipart),Form("hEffPosMean%i",ipart),nbinspt,xbins);
+ hEffNegMean[ipart]=new TH1F(Form("hEffNegMean%i",ipart),Form("hEffNegMean%i",ipart),nbinspt,xbins);
+ }
+
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ hHistPosNSigma[ipart]->Sumw2();
+ hHistNegNSigma[ipart]->Sumw2();
+ hHistPosNSigmaMean[ipart]->Sumw2();
+ hHistNegNSigmaMean[ipart]->Sumw2();
+ hHistPosNSigmaPrimMC[ipart]->Sumw2();
+ hHistNegNSigmaPrimMC[ipart]->Sumw2();
+ hHistPosNSigmaMC[ipart]->Sumw2();
+ hHistNegNSigmaMC[ipart]->Sumw2();
+ hHistPosNSigmaPrimMCMean[ipart]->Sumw2();
+ hHistNegNSigmaPrimMCMean[ipart]->Sumw2();
+ hHistPosNSigmaPrim[ipart]->Sumw2();
+ hHistNegNSigmaPrim[ipart]->Sumw2();
+ hHistPosNSigmaPrimMean[ipart]->Sumw2();
+ hHistNegNSigmaPrimMean[ipart]->Sumw2();
+ hHistPosNSigmaMCMean[ipart]->Sumw2();
+ hHistNegNSigmaMCMean[ipart]->Sumw2();
+
+ hCorrFactPos[ipart]->Divide(hHistPosNSigma[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+ hCorrFactNeg[ipart]->Divide(hHistNegNSigma[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+ hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+ hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+
+ // hCorrFactPos[ipart]->Divide(hHistPosNSigmaPrim[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+ // hCorrFactNeg[ipart]->Divide(hHistNegNSigmaPrim[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+ // hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaPrimMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+ // hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaPrimMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+
+ //hCorrFactPos[ipart]->Divide(hHistPosNSigmaPrimMC[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+ //hCorrFactNeg[ipart]->Divide(hHistNegNSigmaPrimMC[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+ //hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaPrimMCMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
+ //hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaPrimMCMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
+
+ hCorrFactPos[ipart]->SetTitle(Form("Correction Factor NSigma Pos%i",ipart));
+ hCorrFactNeg[ipart]->SetTitle(Form("Correction Factor NSigma Neg%i",ipart));
+ hCorrFactPosMean[ipart]->SetTitle(Form("Corr Fact NSigma Pos%i",ipart));
+ hCorrFactNegMean[ipart]->SetTitle(Form("Corr Fact NSigma Neg%i",ipart));
+ hCorrFactPos[ipart]->SetName(Form("hCorrFactPos%i",ipart));
+ hCorrFactNeg[ipart]->SetName(Form("hCorrFactNeg%i",ipart));
+ hCorrFactPosMean[ipart]->SetName(Form("hCorrFactPosMean%i",ipart));
+ hCorrFactNegMean[ipart]->SetName(Form("hCorrFactNegMean%i",ipart));
+ SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hCorrFactPos[ipart]);
+ SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hCorrFactNeg[ipart]);
+ SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hCorrFactPosMean[ipart]);
+ SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hCorrFactNegMean[ipart]);
+ ResetOutOfRange(hCorrFactPos[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hCorrFactNeg[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hCorrFactPosMean[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hCorrFactNegMean[ipart],ipart,lowRange,upRange);
+
+ hHistPosNSigmaPrimMC[ipart]->Sumw2();
+ hHistNegNSigmaPrimMC[ipart]->Sumw2();
+ hHistPosNSigmaPrimMCMean[ipart]->Sumw2();
+ hHistNegNSigmaPrimMCMean[ipart]->Sumw2();
+ hEffPos[ipart]->Divide(hHistPosNSigmaPrimMC[ipart],fHistPrimMCpos[ipart],1,1,"B");
+ hEffNeg[ipart]->Divide(hHistNegNSigmaPrimMC[ipart],fHistPrimMCneg[ipart],1,1,"B");
+ hEffPosMean[ipart]->Divide(hHistPosNSigmaPrimMCMean[ipart],fHistPrimMCpos[ipart],1,1,"B");
+ hEffNegMean[ipart]->Divide(hHistNegNSigmaPrimMCMean[ipart],fHistPrimMCneg[ipart],1,1,"B");
+ hEffPos[ipart]->SetTitle(Form("Efficiency NSigma Pos%i",ipart));
+ hEffNeg[ipart]->SetTitle(Form("Efficiency NSigma Neg%i",ipart));
+ hEffPosMean[ipart]->SetTitle(Form("Efficiency NSigma Asymm Pos%i",ipart));
+ hEffNegMean[ipart]->SetTitle(Form("Efficiency NSigma Asymm Neg%i",ipart));
+ hEffPos[ipart]->SetName(Form("hEffPos%i",ipart));
+ hEffNeg[ipart]->SetName(Form("hEffNeg%i",ipart));
+ hEffPosMean[ipart]->SetName(Form("hEffPosMean%i",ipart));
+ hEffNegMean[ipart]->SetName(Form("hEffNegMean%i",ipart));
+
+ SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hEffPos[ipart]);
+ SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hEffNeg[ipart]);
+ SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hEffPosMean[ipart]);
+ SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hEffNegMean[ipart]);
+ ResetOutOfRange(hEffPos[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hEffNeg[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hEffPosMean[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hEffNegMean[ipart],ipart,lowRange,upRange);
+
+
+ hCorrFactPos[ipart]->SetMinimum(0.);
+ hCorrFactNeg[ipart]->SetMinimum(0.);
+ hCorrFactPosMean[ipart]->SetMinimum(0.);
+ hCorrFactNegMean[ipart]->SetMinimum(0.);
+ hEffPos[ipart]->SetMinimum(0.);
+ hEffNeg[ipart]->SetMinimum(0.);
+ hEffPosMean[ipart]->SetMinimum(0.);
+ hEffNegMean[ipart]->SetMinimum(0.);
+ hCorrFactPos[ipart]->SetMaximum(1.);
+ hCorrFactNeg[ipart]->SetMaximum(1.);
+ hCorrFactPosMean[ipart]->SetMaximum(1.);
+ hCorrFactNegMean[ipart]->SetMaximum(1.);
+ hEffPos[ipart]->SetMaximum(1.);
+ hEffNeg[ipart]->SetMaximum(1.);
+ hEffPosMean[ipart]->SetMaximum(1.);
+ hEffNegMean[ipart]->SetMaximum(1.);
+ }
+
+ if(MakeMCPlots){
+ TCanvas *cCorr=new TCanvas("cCorr","cCorr");
+ cCorr->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hCorrFactPos[ipart]->Draw("P");
+ else hCorrFactPos[ipart]->Draw("PSAME");
+ hCorrFactNeg[ipart]->Draw("PSAME");
+ }
+ cCorr->BuildLegend();
+
+ TCanvas *cCorrMean=new TCanvas("cCorrMean","cCorrMean");
+ cCorrMean->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0){
+ TH1F * hempty=(TH1F*)hCorrFactPosMean[ipart]->Clone();
+ hempty->SetXTitle("P_{T} (GeV/c)");
+ hempty->SetYTitle("Acc x #epsilon x Contamination");
+ hempty->Reset("all");
+ hempty->Draw();
+
+ }
+ hCorrFactPosMean[ipart]->Draw("PSAME");
+ hCorrFactNegMean[ipart]->Draw("PSAME");
+ }
+ cCorrMean->BuildLegend();
+
+ TCanvas *cCorrcfPos=new TCanvas("cCorrcfPos","cCorrcfPos");
+ cCorrcfPos->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hCorrFactPos[ipart]->Draw("P");
+ else hCorrFactPos[ipart]->Draw("PSAME");
+ hCorrFactPosMean[ipart]->Draw("PSAME");
+ }
+ cCorrcfPos->BuildLegend();
+
+ TCanvas *cCorrcfNeg=new TCanvas("cCorrcfNeg","cCorrcfNeg");
+ cCorrcfNeg->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hCorrFactNeg[ipart]->Draw("P");
+ else hCorrFactNeg[ipart]->Draw("PSAME");
+ hCorrFactNegMean[ipart]->Draw("PSAME");
+ }
+ cCorrcfNeg->BuildLegend();
+
+ TCanvas *cEff=new TCanvas("cEff","cEff");
+ cEff->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hEffPos[ipart]->Draw("P");
+ else hEffPos[ipart]->Draw("PSAME");
+ hEffNeg[ipart]->Draw("PSAME");
+ }
+ cEff->BuildLegend();
+
+ TCanvas *cEffMean=new TCanvas("cEffMean","cEffMean");
+ cEffMean->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hEffPosMean[ipart]->Draw("P");
+ else hEffPosMean[ipart]->Draw("PSAME");
+ hEffNegMean[ipart]->Draw("PSAME");
+ }
+ cEffMean->BuildLegend();
+
+ RatioPlot(hEffNegMean,hEffPosMean,"Neg","Pos","Efficiency",Low,Up);
+
+
+ TCanvas *cEffcfPos=new TCanvas("cEffcfPos","cEffcfPos");
+ cEffcfPos->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hEffPos[ipart]->Draw("P");
+ else hEffPos[ipart]->Draw("PSAME");
+ hEffPosMean[ipart]->Draw("PSAME");
+ }
+ cEffcfPos->BuildLegend();
+
+ TCanvas *cEffcfNeg=new TCanvas("cEffcfNeg","cEffcfNeg");
+ cEffcfNeg->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hEffNeg[ipart]->Draw("P");
+ else hEffNeg[ipart]->Draw("PSAME");
+ hEffNegMean[ipart]->Draw("PSAME");
+ }
+ cEffcfNeg->BuildLegend();
+ }
+
+
+
+
+ for(Int_t ipart=0;ipart<=2;ipart++) {
+ hHistPosNSigma[ipart]->Scale(normMC,"width");
+ hHistNegNSigma[ipart]->Scale(normMC,"width");
+ hHistPosNSigmaMean[ipart]->Scale(normMC,"width");
+ hHistNegNSigmaMean[ipart]->Scale(normMC,"width");
+ }
+
+
+ //////////////////////////////Prim/All Sec/All SecStr/All
+ TH1F *hPrimPos[3];
+ TH1F *hPrimNeg[3];
+ TH1F *hSecPos[3];
+ TH1F *hSecNeg[3];
+ TH1F *hSecStrPos[3];
+ TH1F *hSecStrNeg[3];
+ TH1F *hAllPos[3];
+ TH1F *hAllNeg[3];
+ TH1F *hAllSecPos[3];
+ TH1F *hAllSecNeg[3];
+
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ hPrimPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCposReco%i",ipart));
+ hPrimNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCnegReco%i",ipart));
+ hSecPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecMatMCposReco%i",ipart));
+ hSecNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecMatMCnegReco%i",ipart));
+ hSecStrPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecStrMCposReco%i",ipart));
+ hSecStrNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecStrMCnegReco%i",ipart));
+
+ hAllSecPos[ipart]=(TH1F*)hSecPos[ipart]->Clone(Form("hAllSecPos%i",ipart));
+ hAllSecPos[ipart]->Add(hSecStrPos[ipart]);
+ hAllSecNeg[ipart]=(TH1F*)hSecNeg[ipart]->Clone(Form("hAllSecNeg%i",ipart));
+ hAllSecNeg[ipart]->Add(hSecStrNeg[ipart]);
+ hAllPos[ipart]=(TH1F*)hPrimPos[ipart]->Clone(Form("hAllPos%i",ipart));
+ hAllPos[ipart]->Add(hSecPos[ipart]);
+ hAllPos[ipart]->Add(hSecStrPos[ipart]);
+ hAllNeg[ipart]=(TH1F*)hPrimNeg[ipart]->Clone(Form("hAllNeg%i",ipart));
+ hAllNeg[ipart]->Add(hSecNeg[ipart]);
+ hAllNeg[ipart]->Add(hSecStrNeg[ipart]);
+
+ hPrimPos[ipart]->Divide(hAllPos[ipart]);
+ hSecPos[ipart]->Divide(hAllPos[ipart]);
+ hSecStrPos[ipart]->Divide(hAllPos[ipart]);
+ hPrimNeg[ipart]->Divide(hAllNeg[ipart]);
+ hSecNeg[ipart]->Divide(hAllNeg[ipart]);
+ hSecStrNeg[ipart]->Divide(hAllNeg[ipart]);
+
+ hAllSecPos[ipart]->Divide(hAllPos[ipart]);
+ hAllSecNeg[ipart]->Divide(hAllNeg[ipart]);
+ for(Int_t binnum=1;binnum<=hAllSecNeg[ipart]->GetNbinsX();binnum++){
+ hAllSecPos[ipart]->SetBinContent(binnum,1-hAllSecPos[ipart]->GetBinContent(binnum));
+ hAllSecNeg[ipart]->SetBinContent(binnum,1-hAllSecNeg[ipart]->GetBinContent(binnum));
+ }
+ }
+
+ /////////////////////////////////////////////
+ TFile *outMC2=new TFile(foutMC.Data(),"update");
+ TList *lMC=new TList();
+ lMC->SetOwner();
+ lMC->SetName(Form("MC_Mult%ito%i",lowmult,upmult));
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ lMC->Add(hCorrFactPos[ipart]);
+ lMC->Add(hCorrFactNeg[ipart]);
+ lMC->Add(hCorrFactPosMean[ipart]);
+ lMC->Add(hCorrFactNegMean[ipart]);
+ lMC->Add(hEffPos[ipart]);
+ lMC->Add(hEffNeg[ipart]);
+ lMC->Add(hEffPosMean[ipart]);
+ lMC->Add(hEffNegMean[ipart]);
+ lMC->Add(hHistPosNSigma[ipart]);
+ lMC->Add(hHistNegNSigma[ipart]);
+ lMC->Add(hHistPosNSigmaMean[ipart]);
+ lMC->Add(hHistNegNSigmaMean[ipart]);
+ lMC->Add(hPrimPos[ipart]);
+ lMC->Add(hPrimNeg[ipart]);
+ lMC->Add(hSecPos[ipart]);
+ lMC->Add(hSecNeg[ipart]);
+ lMC->Add(hSecStrPos[ipart]);
+ lMC->Add(hSecStrNeg[ipart]);
+ lMC->Add(hAllSecPos[ipart]);
+ lMC->Add(hAllSecNeg[ipart]);
+ }
+ lMC->Write(Form("MC_Mult%ito%i",lowmult,upmult),1);
+ outMC2->Close();
+ delete outMC2;
+
+
+ TString lnameDATA=Form("clistITSsaMult%ito%i",lowmult,upmult);
+ Printf("\n\n------------------ READING DATA %s -------------------\n",lnameDATA.Data());
+ TFile *finDATA = new TFile(fnameDATA.Data());
+ TDirectory *dirFileDATA=(TDirectory*)finDATA->Get("PWG2SpectraITSsa");
+ TList *cOutputDATA = (TList*)dirFileDATA->Get(lnameDATA.Data());
+
+
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ hHistPosNSigmaDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistPosNSigma%i",ipart));
+ hHistNegNSigmaDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistNegNSigma%i",ipart));
+ hHistPosNSigmaMeanDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistPosNSigmaMean%i",ipart));
+ hHistNegNSigmaMeanDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistNegNSigmaMean%i",ipart));
+
+ hITSsaPos[ipart]=(TH1F*)hHistPosNSigmaDATA[ipart]->Clone(Form("hITSsaPos%i",ipart));
+ hITSsaNeg[ipart]=(TH1F*)hHistNegNSigmaDATA[ipart]->Clone(Form("hITSsaNeg%i",ipart));
+ hITSsaPosMean[ipart]=(TH1F*)hHistPosNSigmaMeanDATA[ipart]->Clone(Form("hITSsaPosAsymm%i",ipart));
+ hITSsaNegMean[ipart]=(TH1F*)hHistNegNSigmaMeanDATA[ipart]->Clone(Form("hITSsaNegAsymm%i",ipart));
+ hITSsaRawPos[ipart]=(TH1F*)hHistPosNSigmaDATA[ipart]->Clone(Form("hITSsaRawPos%i",ipart));
+ hITSsaRawNeg[ipart]=(TH1F*)hHistNegNSigmaDATA[ipart]->Clone(Form("hITSsaRawNeg%i",ipart));
+ hITSsaRawPosMean[ipart]=(TH1F*)hHistPosNSigmaMeanDATA[ipart]->Clone(Form("hITSsaRawPosAsymm%i",ipart));
+ hITSsaRawNegMean[ipart]=(TH1F*)hHistNegNSigmaMeanDATA[ipart]->Clone(Form("hITSsaRawNegAsymm%i",ipart));
+
+ hITSsaPos[ipart]->Sumw2();
+ hITSsaNeg[ipart]->Sumw2();
+ hITSsaPosMean[ipart]->Sumw2();
+ hITSsaNegMean[ipart]->Sumw2();
+ hITSsaRawPos[ipart]->Sumw2();
+ hITSsaRawNeg[ipart]->Sumw2();
+ hITSsaRawPosMean[ipart]->Sumw2();
+ hITSsaRawNegMean[ipart]->Sumw2();
+
+ hITSsaPos[ipart]->Divide(hCorrFactPos[ipart]);
+ hITSsaNeg[ipart]->Divide(hCorrFactNeg[ipart]);
+ hITSsaPosMean[ipart]->Divide(hCorrFactPosMean[ipart]);
+ hITSsaNegMean[ipart]->Divide(hCorrFactNegMean[ipart]);
+
+ SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hITSsaPos[ipart]);
+ SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hITSsaNeg[ipart]);
+ SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hITSsaPosMean[ipart]);
+ SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hITSsaNegMean[ipart]);
+ SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hITSsaRawPos[ipart]);
+ SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hITSsaRawNeg[ipart]);
+ SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hITSsaRawPosMean[ipart]);
+ SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hITSsaRawNegMean[ipart]);
+ }
+
+ if(!CorrectOnlyForEfficiency){
+
+
+ if(DCACorrection){
+ printf("- Correction from DCA\n");
+
+ // Correction factor based on fit to DCA distr for P and Pibar
+ TString fnameDCAPosP=Form("./DCACorr%s_%s_%.0fDCA_PosP_TFraction.root",period.Data(),MCname.Data(),DCAcut);
+ TFile *fDCAPosP= new TFile(fnameDCAPosP.Data());
+ TH1F * hDCAPosP= (TH1F*)fDCAPosP->Get("hPrimAllDATAMC");
+ //TH1F * hDCAPosP= (TH1F*)fDCAPosP->Get("hPrimAllDATA");
+ TString fnameDCANegP=Form("./DCACorr%s_%s_%.0fDCA_NegP_TFraction.root",period.Data(),MCname.Data(),DCAcut);
+ TFile *fDCANegP= new TFile(fnameDCANegP.Data());
+ TH1F * hDCANegP= (TH1F*)fDCANegP->Get("hPrimAllDATAMC");
+ //TH1F * hDCANegP= (TH1F*)fDCANegP->Get("hPrimAllDATA");
+
+ for(Int_t pbin=0;pbin<=hDCAPosP->GetNbinsX();pbin++)hDCAPosP->SetBinError(pbin,0);
+ for(Int_t pbin=0;pbin<=hDCANegP->GetNbinsX();pbin++)hDCANegP->SetBinError(pbin,0);
+
+ hITSsaPos[2]->Multiply(hDCAPosP);
+ hITSsaNeg[2]->Multiply(hDCANegP);
+ hITSsaPosMean[2]->Multiply(hDCAPosP);
+ hITSsaNegMean[2]->Multiply(hDCANegP);
+
+ hDCAPosP->Draw("");
+ hDCANegP->Draw("same");
+
+ }
+ else{
+ Printf("NO DCA correction for protons");
+ }
+
+ if(GFCorrection){
+ printf("- Geant/FLUKA correction\n");
+ //GeantFluka correction for Pi and K
+ TString fnameGeanFlukaPi="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSection.211.root";
+ TFile *fGeanFlukaPi= new TFile(fnameGeanFlukaPi.Data());
+ TH1F *hGeantFlukaPiPos=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionParticles");
+ TH1F *hGeantFlukaPiNeg=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionAntiParticles");
+ for(Int_t binPi=0;binPi<=hITSsaPos[0]->GetNbinsX();binPi++){
+ Float_t FlukaCorrPiPos=hGeantFlukaPiPos->GetBinContent(hGeantFlukaPiPos->FindBin(hITSsaPos[0]->GetBinCenter(binPi)));
+ Float_t FlukaCorrPiNeg=hGeantFlukaPiNeg->GetBinContent(hGeantFlukaPiNeg->FindBin(hITSsaNeg[0]->GetBinCenter(binPi)));
+ hITSsaPos[0]->SetBinContent(binPi,hITSsaPos[0]->GetBinContent(binPi)*FlukaCorrPiPos);
+ hITSsaPos[0]->SetBinError(binPi,hITSsaPos[0]->GetBinError(binPi)*FlukaCorrPiPos);
+ hITSsaNeg[0]->SetBinContent(binPi,hITSsaNeg[0]->GetBinContent(binPi)*FlukaCorrPiNeg);
+ hITSsaNeg[0]->SetBinError(binPi,hITSsaNeg[0]->GetBinError(binPi)*FlukaCorrPiNeg);
+ hITSsaPosMean[0]->SetBinContent(binPi,hITSsaPosMean[0]->GetBinContent(binPi)*FlukaCorrPiPos);
+ hITSsaPosMean[0]->SetBinError(binPi,hITSsaPosMean[0]->GetBinError(binPi)*FlukaCorrPiPos);
+ hITSsaNegMean[0]->SetBinContent(binPi,hITSsaNegMean[0]->GetBinContent(binPi)*FlukaCorrPiNeg);
+ hITSsaNegMean[0]->SetBinError(binPi,hITSsaNegMean[0]->GetBinError(binPi)*FlukaCorrPiNeg);
+ }
+
+ TString fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSection.321.root";
+ TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
+ TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
+ TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
+ for(Int_t binK=0;binK<=hITSsaPos[1]->GetNbinsX();binK++){
+ Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(hITSsaPos[1]->GetBinCenter(binK)));
+ Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(hITSsaNeg[1]->GetBinCenter(binK)));
+ hITSsaPos[1]->SetBinContent(binK,hITSsaPos[1]->GetBinContent(binK)*FlukaCorrKPos);
+ hITSsaNeg[1]->SetBinContent(binK,hITSsaNeg[1]->GetBinContent(binK)*FlukaCorrKNeg);
+ hITSsaPos[1]->SetBinError(binK,hITSsaPos[1]->GetBinError(binK)*FlukaCorrKPos);
+ hITSsaNeg[1]->SetBinError(binK,hITSsaNeg[1]->GetBinError(binK)*FlukaCorrKNeg);
+ hITSsaPosMean[1]->SetBinContent(binK,hITSsaPosMean[1]->GetBinContent(binK)*FlukaCorrKPos);
+ hITSsaNegMean[1]->SetBinContent(binK,hITSsaNegMean[1]->GetBinContent(binK)*FlukaCorrKNeg);
+ hITSsaPosMean[1]->SetBinError(binK,hITSsaPosMean[1]->GetBinError(binK)*FlukaCorrKPos);
+ hITSsaNegMean[1]->SetBinError(binK,hITSsaNegMean[1]->GetBinError(binK)*FlukaCorrKNeg);
+ }
+
+ //Geant Fluka for P
+ // ITS specific file for protons/antiprotons
+ Int_t kNCharge=2;
+ Int_t kPos=0;
+ Int_t kNeg=1;
+ TFile* fITS = new TFile ("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSectionITS_20100719.root");
+ // TH2D * hCorrFlukaITS[kNCharge];
+ TH2D * hCorrFlukaITS[2];
+ hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
+ hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
+ for(Int_t icharge = 0; icharge < kNCharge; icharge++){
+ Int_t nbins = hITSsaPos[2]->GetNbinsX();
+ Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
+ //for(Int_t binP=0;binP<=hITSsaPos[2]->GetNbinsX();binP++){
+ for(Int_t ibin = 0; ibin < nbins; ibin++){
+ Float_t pt = hITSsaPos[2]->GetBinCenter(ibin);
+ Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
+ Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
+ if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
+ if (pt > maxPtCorrection) pt = maxPtCorrection;
+ Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
+ cout<<correction<<" charge "<<icharge<<endl;
+ if(icharge==0){
+ if (correction != 0) {// If the bin is empty this is a 0
+ hITSsaPos[2]->SetBinContent(ibin,hITSsaPos[2]->GetBinContent(ibin)*correction);
+ hITSsaPos[2]->SetBinError(ibin,hITSsaPos[2]->GetBinError (ibin)*correction);
+ hITSsaPosMean[2]->SetBinContent(ibin,hITSsaPosMean[2]->GetBinContent(ibin)*correction);
+ hITSsaPosMean[2]->SetBinError(ibin,hITSsaPosMean[2]->GetBinError (ibin)*correction);
+
+ }else if (hITSsaPos[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+ cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << endl;
+ cout << " Bin content: " << hITSsaPos[2]->GetBinContent(ibin) << endl;
+ }
+ }
+ if(icharge==1){
+ if (correction != 0) {// If the bin is empty this is a 0
+ hITSsaNeg[2]->SetBinContent(ibin,hITSsaNeg[2]->GetBinContent(ibin)*correction);
+ hITSsaNeg[2]->SetBinError(ibin,hITSsaNeg[2]->GetBinError (ibin)*correction);
+ hITSsaNegMean[2]->SetBinContent(ibin,hITSsaNegMean[2]->GetBinContent(ibin)*correction);
+ hITSsaNegMean[2]->SetBinError(ibin,hITSsaNegMean[2]->GetBinError (ibin)*correction);
+ }else if (hITSsaNeg[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
+ cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, ITS, " << endl;
+ cout << " Bin content: " << hITSsaNeg[2]->GetBinContent(ibin) << endl;
+ }
+ }
+
+ }
+ }
+ }else{
+ Printf("SPECTRA ARE NOT CORRECTED FOR GEANT/FLUKA");
+
+ }
+
+ }else{
+ Printf("SPECTRA ARE NOT CORRECTED FOR DCA FIT + GEANT/FLUKA");
+ }
+
+ Float_t norm;
+ TH1F *fHistNEvents =0x0;
+ if(MarekNorm){
+ Printf("Normalization to the number of events with good vertex");
+ TFile *finMAREK = new TFile(fnameDATA.Data());
+ TDirectory *dirFileMAREK=(TDirectory*)finMAREK->Get("PWG2SpectraITSTPC");
+ if(!dirFileMAREK)Printf("!dirFileMAREK");
+ TString lnameMAREK=Form("outputlow%dup%dHI0",lowmult,upmult);
+ TList *cOutputMAREK = (TList*)dirFileMAREK->Get(lnameMAREK.Data());
+ if(!cOutputMAREK)Printf("!cOutputMAREK");
+ TH1F *StatsHistMAREK=(TH1F*)cOutputMAREK->FindObject("StatsHist");
+ norm=(Float_t)1/StatsHistMAREK->GetBinContent(3);
+ }
+ else{
+ Printf("Normalization to the number of events after phys sel");
+ fHistNEvents = (TH1F*)cOutputDATA->FindObject("fHistNEvents");
+ norm=(Float_t)1/fHistNEvents->GetBinContent(fHistNEvents->FindBin(1.));
+ }
+ if(NormalizeData)printf("- Normalization (NEvents): %.0f \n",1/norm);
+ else printf("- Normalization (NEvents): %.0f !!!!!!!SPECTRA ARE NOT NORMALIZED!!!!!!\n",fHistNEvents->GetBinContent(fHistNEvents->FindBin(1.)));
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(NormalizeData){
+ hITSsaPos[ipart]->Scale(norm,"width");
+ hITSsaNeg[ipart]->Scale(norm,"width");
+ hITSsaPosMean[ipart]->Scale(norm,"width");
+ hITSsaNegMean[ipart]->Scale(norm,"width");
+ hITSsaRawPos[ipart]->Scale(norm,"width");
+ hITSsaRawNeg[ipart]->Scale(norm,"width");
+ hITSsaRawPosMean[ipart]->Scale(norm,"width");
+ hITSsaRawNegMean[ipart]->Scale(norm,"width");
+ }
+ else{
+ hITSsaPos[ipart]->Scale(1,"width");
+ hITSsaNeg[ipart]->Scale(1,"width");
+ hITSsaPosMean[ipart]->Scale(1,"width");
+ hITSsaNegMean[ipart]->Scale(1,"width");
+ hITSsaRawPos[ipart]->Scale(1,"width");
+ hITSsaRawNeg[ipart]->Scale(1,"width");
+ hITSsaRawPosMean[ipart]->Scale(1,"width");
+ hITSsaRawNegMean[ipart]->Scale(1,"width");
+ }
+ if(NormalizeToINEL){
+ Float_t INEL=0.86;
+ printf("Normalization to INEL : %.2f \n",INEL);
+ hITSsaPos[ipart]->Scale(INEL,"");
+ hITSsaNeg[ipart]->Scale(INEL,"");
+ hITSsaPosMean[ipart]->Scale(INEL,"");
+ hITSsaNegMean[ipart]->Scale(INEL,"");
+ hITSsaRawPos[ipart]->Scale(INEL,"");
+ hITSsaRawNeg[ipart]->Scale(INEL,"");
+ hITSsaRawPosMean[ipart]->Scale(INEL,"");
+ hITSsaRawNegMean[ipart]->Scale(INEL,"");
+ }
+ else printf("Normalization to INEL ignored \n");
+
+ hITSsaPos[ipart]->SetMinimum(0.00001);
+ hITSsaNeg[ipart]->SetMinimum(0.00001);
+ hITSsaPosMean[ipart]->SetMinimum(0.00001);
+ hITSsaNegMean[ipart]->SetMinimum(0.00001);
+ hITSsaRawPos[ipart]->SetMinimum(0.00001);
+ hITSsaRawNeg[ipart]->SetMinimum(0.00001);
+ hITSsaRawPosMean[ipart]->SetMinimum(0.00001);
+ hITSsaRawNegMean[ipart]->SetMinimum(0.00001);
+ hITSsaPos[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+ hITSsaNeg[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+ hITSsaPosMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+ hITSsaNegMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+ hITSsaRawPos[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+ hITSsaRawNeg[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+ hITSsaRawPosMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+ hITSsaRawNegMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
+ hITSsaPos[ipart]->SetTitle(Form("ITSsa NSigma Symm Pos%i",ipart));
+ hITSsaNeg[ipart]->SetTitle(Form("ITSsa NSigma Symm Neg%i",ipart));
+ hITSsaPosMean[ipart]->SetTitle(Form("ITSsa NSigma Asymm Pos%i",ipart));
+ hITSsaNegMean[ipart]->SetTitle(Form("ITSsa NSigma Asymm Neg%i",ipart));
+ hITSsaRawPos[ipart]->SetTitle(Form("ITSsaRaw NSigma Symm Pos%i",ipart));
+ hITSsaRawNeg[ipart]->SetTitle(Form("ITSsaRaw NSigma Symm Neg%i",ipart));
+ hITSsaRawPosMean[ipart]->SetTitle(Form("ITSsaRaw NSigma Asymm Pos%i",ipart));
+ hITSsaRawNegMean[ipart]->SetTitle(Form("ITSsaRaw NSigma Asymm Neg%i",ipart));
+
+
+ ResetOutOfRange(hITSsaPos[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hITSsaNeg[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hITSsaPosMean[ipart],ipart,lowRange,upRange);
+ ResetOutOfRange(hITSsaNegMean[ipart],ipart,lowRange,upRange);
+
+ }
+ TH1F *hSystematicsITSsaPos[3];
+ TH1F *hSystematicsITSsaNeg[3];
+ TH1F *hSystematicsITSsaPosMean[3];
+ TH1F *hSystematicsITSsaNegMean[3];
+
+ if(MakeSystErr){
+ printf("- making Systematic Error\n");
+ //preliminary systematic error
+ TFile *fsys = new TFile("RootFilesSystError/ITSsa-systematics-20101014.root");
+ //fsys->Print("all");
+ //cout<<fsys->GetSize()<<endl;
+ TH1F* hSystPos[3];
+ TH1F* hSystNeg[3];
+ hSystPos[0]=(TH1F*)fsys->Get("hSystTotPosPion");
+ hSystNeg[0]=(TH1F*)fsys->Get("hSystTotNegPion");
+ hSystPos[1]=(TH1F*)fsys->Get("hSystTotPosKaon");
+ hSystNeg[1]=(TH1F*)fsys->Get("hSystTotNegKaon");
+ hSystPos[2]=(TH1F*)fsys->Get("hSystTotPosProton");
+ hSystNeg[2]=(TH1F*)fsys->Get("hSystTotNegProton");
+
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ hSystematicsITSsaPos[ipart]=(TH1F*)hITSsaPos[ipart]->Clone(Form("hSystematicsITSsaPos%i",ipart));
+ hSystematicsITSsaNeg[ipart]=(TH1F*)hITSsaNeg[ipart]->Clone(Form("hSystematicsITSsaNeg%i",ipart));
+ hSystematicsITSsaPosMean[ipart]=(TH1F*)hITSsaPosMean[ipart]->Clone(Form("hSystematicsITSsaPosAsymm%i",ipart));
+ hSystematicsITSsaNegMean[ipart]=(TH1F*)hITSsaNegMean[ipart]->Clone(Form("hSystematicsITSsaNegAsymm%i",ipart));
+
+ for(Int_t ibin=0;ibin<=hSystematicsITSsaPos[ipart]->GetNbinsX();ibin++){
+ Float_t syserr=hSystPos[ipart]->GetBinContent(hSystPos[ipart]->FindBin(hSystematicsITSsaPos[ipart]->GetBinCenter(ibin)));
+ hSystematicsITSsaPos[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaPos[ipart]->GetBinContent(ibin));
+ hSystematicsITSsaPosMean[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaPosMean[ipart]->GetBinContent(ibin));
+ }
+ for(Int_t ibin=0;ibin<=hSystematicsITSsaNeg[ipart]->GetNbinsX();ibin++){
+ Float_t syserr=hSystNeg[ipart]->GetBinContent(hSystNeg[ipart]->FindBin(hSystematicsITSsaNeg[ipart]->GetBinCenter(ibin)));
+ hSystematicsITSsaNeg[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaNeg[ipart]->GetBinContent(ibin));
+ hSystematicsITSsaNegMean[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaNegMean[ipart]->GetBinContent(ibin));
+ }
+ }
+ }
+
+ TCanvas *cSpectraSymm=new TCanvas("cSpectraSymm","cSpectraSymm");
+ cSpectraSymm->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hITSsaPos[ipart]->DrawCopy("P");
+ else hITSsaPos[ipart]->DrawCopy("PSAME");
+ hITSsaNeg[ipart]->DrawCopy("PSAME");
+ }
+ cSpectraSymm->BuildLegend();
+
+ TCanvas *cSpectraAsymm=new TCanvas("cSpectraAsymm","cSpectraAsymm");
+ cSpectraAsymm->SetGridy();
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ if(ipart==0)hITSsaPosMean[ipart]->DrawCopy("P");
+ else hITSsaPosMean[ipart]->DrawCopy("PSAME");
+ hITSsaNegMean[ipart]->DrawCopy("PSAME");
+ }
+ cSpectraAsymm->BuildLegend();
+
+
+
+ /////////////////////////// SCALING SPECTRA FOR TRACKING EFFICIENCY !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ Float_t ScalingFactor=1.015;
+ Printf("SCALING FACTOR: %f",ScalingFactor);
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ hITSsaPosMean[ipart]->Scale(ScalingFactor);
+ Printf("Scaling %s by %f",hITSsaPosMean[ipart]->GetName(),ScalingFactor);
+ hITSsaNegMean[ipart]->Scale(ScalingFactor);
+ Printf("Scaling %s by %f",hITSsaNegMean[ipart]->GetName(),ScalingFactor);
+ hITSsaPos[ipart]->Scale(ScalingFactor);
+ Printf("Scaling %s by %f",hITSsaPos[ipart]->GetName(),ScalingFactor);
+ hITSsaNeg[ipart]->Scale(ScalingFactor);
+ Printf("Scaling %s by %f",hITSsaNeg[ipart]->GetName(),ScalingFactor);
+ }
+
+ TFile *outDATA2=new TFile(foutDATA.Data(),"update");
+ TList *lDATA=new TList();
+ lDATA->SetOwner();
+ lDATA->SetName(Form("DATA_Mult%ito%i",lowmult,upmult));
+ if(SaveOnlyAsymm){
+ for(Int_t ipart=0;ipart<=2;ipart++) lDATA->Add(hITSsaPosMean[ipart]);
+ for(Int_t ipart=0;ipart<=2;ipart++) lDATA->Add(hITSsaNegMean[ipart]);
+ }else{
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ lDATA->Add(fHistPrimMCposBefEvSel[ipart]);
+ lDATA->Add(fHistPrimMCnegBefEvSel[ipart]);
+ lDATA->Add(hHistPosNSigmaPrimMean[ipart]);
+ lDATA->Add(hHistNegNSigmaPrimMean[ipart]);
+ lDATA->Add(hHistPosNSigmaPrimMCMean[ipart]);
+ lDATA->Add(hHistNegNSigmaPrimMCMean[ipart]);
+ lDATA->Add(hITSsaPos[ipart]);
+ lDATA->Add(hITSsaNeg[ipart]);
+ lDATA->Add(hITSsaPosMean[ipart]);
+ lDATA->Add(hITSsaNegMean[ipart]);
+ lDATA->Add(hITSsaRawPos[ipart]);
+ lDATA->Add(hITSsaRawNeg[ipart]);
+ lDATA->Add(hITSsaRawPosMean[ipart]);
+ lDATA->Add(hITSsaRawNegMean[ipart]);
+ }
+ if(MakeSystErr){
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ lDATA->Add(hSystematicsITSsaPos[ipart]);
+ lDATA->Add(hSystematicsITSsaNeg[ipart]);
+ lDATA->Add( hSystematicsITSsaPosMean[ipart]);
+ lDATA->Add( hSystematicsITSsaNegMean[ipart]);
+ }
+ }
+ }
+ lDATA->Write(Form("DATA_Mult%ito%i",lowmult,upmult),1);
+ outDATA2->Close();
+ delete outDATA2;
+
+
+ RatioPlot(hITSsaNeg,hITSsaPos,"Neg","Pos","Symm NSigma",Low,Up);
+ RatioPlot(hITSsaNegMean,hITSsaPosMean,"Neg","Pos","Assymm NSigma",Low,Up);
+ RatioPlot(hITSsaPos,hITSsaPosMean,"Symm","Asymm","Positive",Low,Up);
+ RatioPlot(hITSsaNeg,hITSsaNegMean,"Symm","Asymm","Negative",Low,Up);
+
+
+} //end main
+
+
+void ResetOutOfRange(TH1F *histo, Int_t ipart, Double_t lowRange[3], Double_t upRange[3]){
+ // set to -1 the bin contents out of the selcted pt range
+ for(Int_t ibin=0;ibin<=histo->GetNbinsX();ibin++){
+ if(histo->GetBinCenter(ibin)<lowRange[ipart] || histo->GetBinCenter(ibin)>upRange[ipart]) histo->SetBinContent(ibin,-1);
+ }
+}
+
+void RatioPlot(TH1F **num,TH1F **den,TString t1,TString t2,TString opt,Double_t Low[],Double_t Up[])
+{
+ gROOT->SetStyle("Plain");
+ gStyle->SetOptTitle(0);
+ TH1F *Num[3];
+ TString title=Form("%s/%s - %s",t1.Data(),t2.Data(),opt.Data());
+ TCanvas *c=new TCanvas(title.Data(),title.Data());
+ c->Divide(3,2);
+ for(Int_t ipart=0;ipart<=2;ipart++){
+ num[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
+ num[ipart]->SetMarkerColor(2);
+ num[ipart]->SetLineColor(2);
+ num[ipart]->SetStats(kFALSE);
+ den[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
+ den[ipart]->SetMarkerColor(1);
+ den[ipart]->SetLineColor(1);
+ den[ipart]->SetStats(kFALSE);
+ c->cd(ipart+1);
+ gPad->SetGridy();
+ //gPad->SetLogy();
+ num[ipart]->Draw("");
+ den[ipart]->Draw("same");
+ TPaveText *tpave=new TPaveText(0.1,0.2,0.5,0.29,"brNDC");
+ tpave->SetBorderSize(0);
+ tpave->SetFillStyle(0);
+ tpave->SetFillColor(0);
+ tpave->SetTextColor(2);
+ tpave->SetTextSize(.04);
+ TText *txt1=tpave->AddText(t1.Data());
+ txt1->SetTextFont(62);
+ txt1->SetTextColor(2);
+ TText *txt2=tpave->AddText(t2.Data());
+ txt2->SetTextFont(62);
+ txt2->SetTextColor(1);
+ tpave->Draw();
+ Num[ipart]=(TH1F*)num[ipart]->Clone(title.Data());
+ Num[ipart]->SetTitle(title.Data());
+ Num[ipart]->Divide(den[ipart]);
+ Num[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
+ Num[ipart]->SetMarkerColor(4);
+ Num[ipart]->SetLineColor(4);
+ Num[ipart]->SetMaximum(1.2);
+ Num[ipart]->SetMinimum(.8);
+ c->cd(ipart+4);
+ gPad->SetGridy();
+ Num[ipart]->SetStats(kFALSE);
+ Num[ipart]->Draw();
+ TPaveText *tpave_1=new TPaveText(0.3,0.7,0.7,0.99,"brNDC");
+ tpave_1->SetBorderSize(0);
+ tpave_1->SetFillStyle(0);
+ tpave_1->SetFillColor(0);
+ tpave_1->SetTextColor(2);
+ tpave_1->SetTextSize(.04);
+ TText *txt1_1=tpave_1->AddText(title.Data());
+ txt1_1->SetTextFont(62);
+ tpave_1->Draw();
+ // TLine line = TLine(Low[ipart],1,Up[ipart],1);
+ // line.SetLineColor(2);
+ // line.SetLineStyle(2);
+ // line.SetLineWidth(2);
+ // line.DrawClone("same");
+ }
+}
+
+
+
+
+ void SetDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1)
+{
+ h1->SetMarkerStyle(markerstyle);
+ h1->SetMarkerColor(markercolor);
+ h1->SetMarkerSize(markersize);
+ h1->SetLineColor(linecolor);
+ h1->SetLineWidth(linewidth);
+}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff