PWGLF/SPECTRA/PiKaPr/ITSsa/MakeCorrectedSpectraITSsaNsigma.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TH1F.h>
   3 #include <TH2F.h>
   4 #include <TF1.h>
   5 #include <TPad.h>
   6 #include <TGraphErrors.h>
   7 #include <TROOT.h>
   8 #include <TFile.h>
   9 #include <TTree.h>
  10 #include <TGrid.h>
  11 #include <TGridResult.h>
  12 #include <TMath.h>
  13 #include <TCanvas.h>
  14 #include <TStyle.h>
  15 #include <TLatex.h>
  16 #include <TImage.h>
  17 #include <TPaveText.h>
  18 #include <TLine.h>
  19 #endif
  20
  21 void RatioPlot(TH1F **num,TH1F **den,TString t1,TString t2,TString opt,Double_t Low[],Double_t Up[]);
  22 void SetDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1);
  23 void ResetOutOfRange(TH1F *histo, Int_t ipart, Double_t lowRange[3], Double_t upRange[3]);
  24
  25 void MakeCorrectedSpectraITSsaNsigma( TString period="LHC10d",
  26                                       TString MCname="LHC10f6a",
  27                                       Float_t DCAcut=7,
  28                                       Int_t lowmult=-1,
  29                                       Int_t upmult=-1
  30                                       )
  31 {
  32
  33   // parameters
  34
  35   Bool_t MakeMCPlots=1;
  36   Bool_t MakeSystErr=0;
  37   Bool_t NormalizeData=1;
  38   Bool_t NormalizeToINEL=0;
  39   Bool_t SaveOnlyAsymm=0;
  40   Bool_t CorrectOnlyForEfficiency=0;
  41   Bool_t MarekNorm=0;
  42   Bool_t DCACorrection=1;
  43   Bool_t GFCorrection=1;
  44
  45   //MC Histos
  46   TH1F *fHistPrimMCposBefEvSel[3];
  47   TH1F *fHistPrimMCnegBefEvSel[3];
  48   TH1F *fHistPrimMCpos[3];
  49   TH1F *fHistPrimMCneg[3];
  50   TH1F *hHistPosNSigma[3];
  51   TH1F *hHistNegNSigma[3];
  52   TH1F *hHistPosNSigmaMean[3];
  53   TH1F *hHistNegNSigmaMean[3];
  54   TH1F *hHistPosNSigmaPrim[3];
  55   TH1F *hHistNegNSigmaPrim[3];
  56   TH1F *hHistPosNSigmaPrimMean[3];
  57   TH1F *hHistNegNSigmaPrimMean[3];
  58   TH1F *hHistPosNSigmaPrimMC[3];
  59   TH1F *hHistNegNSigmaPrimMC[3];
  60   TH1F *hHistPosNSigmaPrimMCMean[3];
  61   TH1F *hHistNegNSigmaPrimMCMean[3];
  62   TH1F *hHistPosNSigmaMC[3];
  63   TH1F *hHistNegNSigmaMC[3];
  64   TH1F *hHistPosNSigmaMCMean[3];
  65   TH1F *hHistNegNSigmaMCMean[3];
  66   TH1F *hCorrFactPos[3];
  67   TH1F *hCorrFactNeg[3];
  68   TH1F *hCorrFactPosMean[3];
  69   TH1F *hCorrFactNegMean[3];
  70   TH1F *hEffPos[3];
  71   TH1F *hEffNeg[3];
  72   TH1F *hEffPosMean[3];
  73   TH1F *hEffNegMean[3];
  74
  75   //DATA Histos
  76   TH1F *hHistPosNSigmaDATA[3];
  77   TH1F *hHistNegNSigmaDATA[3];
  78   TH1F *hHistPosNSigmaMeanDATA[3];
  79   TH1F *hHistNegNSigmaMeanDATA[3];
  80   TH1F *hITSsaPos[3];
  81   TH1F *hITSsaNeg[3];
  82   TH1F *hITSsaPosMean[3];
  83   TH1F *hITSsaNegMean[3];
  84   TH1F *hITSsaRawPos[3];
  85   TH1F *hITSsaRawNeg[3];
  86   TH1F *hITSsaRawPosMean[3];
  87   TH1F *hITSsaRawNegMean[3];
  88
  89   //  Float_t nmin=1.5;
  90   Double_t Low[3]={0.10,0.2,0.3};
  91   Double_t Up[3]={0.4,0.5,0.7};
  92   Double_t lowRange[3]={0.10,0.2,0.3};
  93   Double_t upRange[3]={0.6,0.6,0.6};
  94
  95
  96
  97   const Int_t nbinspt=22;
  98   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};
  99
 100   TString fname_MC=Form("./%s/AnalysisResults.root",MCname.Data());
 101   cout<<fname_MC.Data()<<endl;
 102
 103   TString foutMC=Form("outSpectraMC_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
 104   TFile *outMC=new TFile(foutMC.Data(),"recreate");
 105   outMC->Close();
 106   delete outMC;
 107
 108   TString foutDATA;
 109   if(CorrectOnlyForEfficiency){
 110     foutDATA=Form("outSpectraData_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
 111   }else{
 112     foutDATA=Form("outCorrSpectraData_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
 113   }
 114
 115   TString fnameDATA=Form("%s/AnalysisResults.root",period.Data());
 116   cout<<fnameDATA.Data()<<endl;
 117   TFile *outDATA=new TFile(foutDATA.Data(),"recreate");
 118   outDATA->Close();
 119   delete outDATA;
 120
 121
 122   TString lnameMC=Form("clistITSsaMult%ito%i",lowmult,upmult);
 123   Printf("\n------------------ READING MC  %s -------------------\n",lnameMC.Data());
 124   TFile *finMC = new TFile(fname_MC.Data());
 125   TDirectory *dirFileMC=(TDirectory*)finMC->Get("PWG2SpectraITSsa");
 126   TList * cOutputMC = (TList*)dirFileMC->Get(lnameMC.Data());
 127   //cOutputMC->Print();
 128
 129   TH1F *fHistNEventsMC = (TH1F*)cOutputMC->FindObject("fHistNEvents");
 130   Float_t normMC=(Float_t)1/fHistNEventsMC->GetBinContent(fHistNEventsMC->FindBin(1.));
 131   printf("- NormalizationMC (NEvents): %.0f \n",fHistNEventsMC->GetBinContent(fHistNEventsMC->FindBin(1.)));
 132
 133   for(Int_t ipart=0;ipart<=2;ipart++){
 134     fHistPrimMCposBefEvSel[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCposBefEvSel%i",ipart));//primaries generated (before event selection)
 135     fHistPrimMCnegBefEvSel[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCnegBefEvSel%i",ipart));
 136     fHistPrimMCpos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCpos%i",ipart)); //primaries generated (after event selection)
 137     fHistPrimMCneg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCneg%i",ipart));
 138     hHistPosNSigma[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigma%i",ipart));//NSigma recontructed, no MCtruth
 139     hHistNegNSigma[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigma%i",ipart));
 140     hHistPosNSigmaPrimMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMC%i",ipart));//NSigma recontructed, MCtruth
 141     hHistNegNSigmaPrimMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMC%i",ipart));
 142     hHistPosNSigmaMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMC%i",ipart));//NSigma recontructed, MCtruth
 143     hHistNegNSigmaMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMC%i",ipart));
 144     hHistPosNSigmaPrim[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrim%i",ipart));//NSigma recontructed, truth
 145     hHistNegNSigmaPrim[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrim%i",ipart));
 146     hHistPosNSigmaMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMean%i",ipart));//NSigmaMean recontructed, no MCtruth
 147     hHistNegNSigmaMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMean%i",ipart));
 148     hHistPosNSigmaPrimMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMCMean%i",ipart));//NSigmaMean recontructed, no MCtruth
 149     hHistNegNSigmaPrimMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMCMean%i",ipart));
 150     hHistPosNSigmaMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMCMean%i",ipart));//NSigmaMean recontructed, no MCtruth
 151     hHistNegNSigmaMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMCMean%i",ipart));
 152     hHistPosNSigmaPrimMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMean%i",ipart));//NSigmaMean recontructed, no MCtruth
 153     hHistNegNSigmaPrimMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMean%i",ipart));
 154   }
 155
 156   for(Int_t ipart=0;ipart<=2;ipart++){
 157     hCorrFactPos[ipart]=new TH1F(Form("hCorrFactPos%i",ipart),Form("hCorrFactPos%i",ipart),nbinspt,xbins);
 158     hCorrFactNeg[ipart]=new TH1F(Form("hCorrFactNeg%i",ipart),Form("hCorrFactNeg%i",ipart),nbinspt,xbins);
 159     hCorrFactPosMean[ipart]=new TH1F(Form("hCorrFactPosMean%i",ipart),Form("hCorrFactPosMean%i",ipart),nbinspt,xbins);
 160     hCorrFactNegMean[ipart]=new TH1F(Form("hCorrFactNegMean%i",ipart),Form("hCorrFactNegMean%i",ipart),nbinspt,xbins);
 161     hEffPos[ipart]=new TH1F(Form("hEffPos%i",ipart),Form("hEffPos%i",ipart),nbinspt,xbins);
 162     hEffNeg[ipart]=new TH1F(Form("hEffNeg%i",ipart),Form("hEffNeg%i",ipart),nbinspt,xbins);
 163     hEffPosMean[ipart]=new TH1F(Form("hEffPosMean%i",ipart),Form("hEffPosMean%i",ipart),nbinspt,xbins);
 164     hEffNegMean[ipart]=new TH1F(Form("hEffNegMean%i",ipart),Form("hEffNegMean%i",ipart),nbinspt,xbins);
 165   }
 166
 167   for(Int_t ipart=0;ipart<=2;ipart++){
 168     hHistPosNSigma[ipart]->Sumw2();
 169     hHistNegNSigma[ipart]->Sumw2();
 170     hHistPosNSigmaMean[ipart]->Sumw2();
 171     hHistNegNSigmaMean[ipart]->Sumw2();
 172     hHistPosNSigmaPrimMC[ipart]->Sumw2();
 173     hHistNegNSigmaPrimMC[ipart]->Sumw2();
 174     hHistPosNSigmaMC[ipart]->Sumw2();
 175     hHistNegNSigmaMC[ipart]->Sumw2();
 176     hHistPosNSigmaPrimMCMean[ipart]->Sumw2();
 177     hHistNegNSigmaPrimMCMean[ipart]->Sumw2();
 178     hHistPosNSigmaPrim[ipart]->Sumw2();
 179     hHistNegNSigmaPrim[ipart]->Sumw2();
 180     hHistPosNSigmaPrimMean[ipart]->Sumw2();
 181     hHistNegNSigmaPrimMean[ipart]->Sumw2();
 182     hHistPosNSigmaMCMean[ipart]->Sumw2();
 183     hHistNegNSigmaMCMean[ipart]->Sumw2();
 184
 185     hCorrFactPos[ipart]->Divide(hHistPosNSigma[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
 186     hCorrFactNeg[ipart]->Divide(hHistNegNSigma[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
 187     hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
 188     hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
 189
 190     // hCorrFactPos[ipart]->Divide(hHistPosNSigmaPrim[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
 191     // hCorrFactNeg[ipart]->Divide(hHistNegNSigmaPrim[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
 192     // hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaPrimMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
 193     // hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaPrimMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
 194
 195     //hCorrFactPos[ipart]->Divide(hHistPosNSigmaPrimMC[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
 196     //hCorrFactNeg[ipart]->Divide(hHistNegNSigmaPrimMC[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
 197     //hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaPrimMCMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
 198     //hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaPrimMCMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
 199
 200     hCorrFactPos[ipart]->SetTitle(Form("Correction Factor NSigma Pos%i",ipart));
 201     hCorrFactNeg[ipart]->SetTitle(Form("Correction Factor NSigma Neg%i",ipart));
 202     hCorrFactPosMean[ipart]->SetTitle(Form("Corr Fact NSigma Pos%i",ipart));
 203     hCorrFactNegMean[ipart]->SetTitle(Form("Corr Fact NSigma Neg%i",ipart));
 204     hCorrFactPos[ipart]->SetName(Form("hCorrFactPos%i",ipart));
 205     hCorrFactNeg[ipart]->SetName(Form("hCorrFactNeg%i",ipart));
 206     hCorrFactPosMean[ipart]->SetName(Form("hCorrFactPosMean%i",ipart));
 207     hCorrFactNegMean[ipart]->SetName(Form("hCorrFactNegMean%i",ipart));
 208     SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hCorrFactPos[ipart]);
 209     SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hCorrFactNeg[ipart]);
 210     SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hCorrFactPosMean[ipart]);
 211     SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hCorrFactNegMean[ipart]);
 212     ResetOutOfRange(hCorrFactPos[ipart],ipart,lowRange,upRange);
 213     ResetOutOfRange(hCorrFactNeg[ipart],ipart,lowRange,upRange);
 214     ResetOutOfRange(hCorrFactPosMean[ipart],ipart,lowRange,upRange);
 215     ResetOutOfRange(hCorrFactNegMean[ipart],ipart,lowRange,upRange);
 216
 217     hHistPosNSigmaPrimMC[ipart]->Sumw2();
 218     hHistNegNSigmaPrimMC[ipart]->Sumw2();
 219     hHistPosNSigmaPrimMCMean[ipart]->Sumw2();
 220     hHistNegNSigmaPrimMCMean[ipart]->Sumw2();
 221     hEffPos[ipart]->Divide(hHistPosNSigmaPrimMC[ipart],fHistPrimMCpos[ipart],1,1,"B");
 222     hEffNeg[ipart]->Divide(hHistNegNSigmaPrimMC[ipart],fHistPrimMCneg[ipart],1,1,"B");
 223     hEffPosMean[ipart]->Divide(hHistPosNSigmaPrimMCMean[ipart],fHistPrimMCpos[ipart],1,1,"B");
 224     hEffNegMean[ipart]->Divide(hHistNegNSigmaPrimMCMean[ipart],fHistPrimMCneg[ipart],1,1,"B");
 225     hEffPos[ipart]->SetTitle(Form("Efficiency NSigma Pos%i",ipart));
 226     hEffNeg[ipart]->SetTitle(Form("Efficiency NSigma Neg%i",ipart));
 227     hEffPosMean[ipart]->SetTitle(Form("Efficiency NSigma Asymm Pos%i",ipart));
 228     hEffNegMean[ipart]->SetTitle(Form("Efficiency NSigma Asymm Neg%i",ipart));
 229     hEffPos[ipart]->SetName(Form("hEffPos%i",ipart));
 230     hEffNeg[ipart]->SetName(Form("hEffNeg%i",ipart));
 231     hEffPosMean[ipart]->SetName(Form("hEffPosMean%i",ipart));
 232     hEffNegMean[ipart]->SetName(Form("hEffNegMean%i",ipart));
 233
 234     SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hEffPos[ipart]);
 235     SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hEffNeg[ipart]);
 236     SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hEffPosMean[ipart]);
 237     SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hEffNegMean[ipart]);
 238     ResetOutOfRange(hEffPos[ipart],ipart,lowRange,upRange);
 239     ResetOutOfRange(hEffNeg[ipart],ipart,lowRange,upRange);
 240     ResetOutOfRange(hEffPosMean[ipart],ipart,lowRange,upRange);
 241     ResetOutOfRange(hEffNegMean[ipart],ipart,lowRange,upRange);
 242
 243
 244     hCorrFactPos[ipart]->SetMinimum(0.);
 245     hCorrFactNeg[ipart]->SetMinimum(0.);
 246     hCorrFactPosMean[ipart]->SetMinimum(0.);
 247     hCorrFactNegMean[ipart]->SetMinimum(0.);
 248     hEffPos[ipart]->SetMinimum(0.);
 249     hEffNeg[ipart]->SetMinimum(0.);
 250     hEffPosMean[ipart]->SetMinimum(0.);
 251     hEffNegMean[ipart]->SetMinimum(0.);
 252     hCorrFactPos[ipart]->SetMaximum(1.);
 253     hCorrFactNeg[ipart]->SetMaximum(1.);
 254     hCorrFactPosMean[ipart]->SetMaximum(1.);
 255     hCorrFactNegMean[ipart]->SetMaximum(1.);
 256     hEffPos[ipart]->SetMaximum(1.);
 257     hEffNeg[ipart]->SetMaximum(1.);
 258     hEffPosMean[ipart]->SetMaximum(1.);
 259     hEffNegMean[ipart]->SetMaximum(1.);
 260   }
 261
 262   if(MakeMCPlots){
 263     TCanvas *cCorr=new TCanvas("cCorr","cCorr");
 264     cCorr->SetGridy();
 265     for(Int_t ipart=0;ipart<=2;ipart++){
 266       if(ipart==0)hCorrFactPos[ipart]->Draw("P");
 267       else hCorrFactPos[ipart]->Draw("PSAME");
 268       hCorrFactNeg[ipart]->Draw("PSAME");
 269     }
 270     cCorr->BuildLegend();
 271
 272     TCanvas *cCorrMean=new TCanvas("cCorrMean","cCorrMean");
 273     cCorrMean->SetGridy();
 274     for(Int_t ipart=0;ipart<=2;ipart++){
 275       if(ipart==0){
 276         TH1F * hempty=(TH1F*)hCorrFactPosMean[ipart]->Clone();
 277         hempty->SetXTitle("P_{T} (GeV/c)");
 278         hempty->SetYTitle("Acc x #epsilon x Contamination");
 279         hempty->Reset("all");
 280         hempty->Draw();
 281
 282       }
 283       hCorrFactPosMean[ipart]->Draw("PSAME");
 284       hCorrFactNegMean[ipart]->Draw("PSAME");
 285     }
 286     cCorrMean->BuildLegend();
 287
 288     TCanvas *cCorrcfPos=new TCanvas("cCorrcfPos","cCorrcfPos");
 289     cCorrcfPos->SetGridy();
 290     for(Int_t ipart=0;ipart<=2;ipart++){
 291       if(ipart==0)hCorrFactPos[ipart]->Draw("P");
 292       else hCorrFactPos[ipart]->Draw("PSAME");
 293       hCorrFactPosMean[ipart]->Draw("PSAME");
 294     }
 295     cCorrcfPos->BuildLegend();
 296
 297     TCanvas *cCorrcfNeg=new TCanvas("cCorrcfNeg","cCorrcfNeg");
 298     cCorrcfNeg->SetGridy();
 299     for(Int_t ipart=0;ipart<=2;ipart++){
 300       if(ipart==0)hCorrFactNeg[ipart]->Draw("P");
 301       else hCorrFactNeg[ipart]->Draw("PSAME");
 302       hCorrFactNegMean[ipart]->Draw("PSAME");
 303     }
 304     cCorrcfNeg->BuildLegend();
 305
 306     TCanvas *cEff=new TCanvas("cEff","cEff");
 307     cEff->SetGridy();
 308     for(Int_t ipart=0;ipart<=2;ipart++){
 309       if(ipart==0)hEffPos[ipart]->Draw("P");
 310       else hEffPos[ipart]->Draw("PSAME");
 311       hEffNeg[ipart]->Draw("PSAME");
 312     }
 313     cEff->BuildLegend();
 314
 315     TCanvas *cEffMean=new TCanvas("cEffMean","cEffMean");
 316     cEffMean->SetGridy();
 317     for(Int_t ipart=0;ipart<=2;ipart++){
 318       if(ipart==0)hEffPosMean[ipart]->Draw("P");
 319       else hEffPosMean[ipart]->Draw("PSAME");
 320       hEffNegMean[ipart]->Draw("PSAME");
 321     }
 322     cEffMean->BuildLegend();
 323
 324     RatioPlot(hEffNegMean,hEffPosMean,"Neg","Pos","Efficiency",Low,Up);
 325
 326
 327     TCanvas *cEffcfPos=new TCanvas("cEffcfPos","cEffcfPos");
 328     cEffcfPos->SetGridy();
 329     for(Int_t ipart=0;ipart<=2;ipart++){
 330       if(ipart==0)hEffPos[ipart]->Draw("P");
 331       else hEffPos[ipart]->Draw("PSAME");
 332       hEffPosMean[ipart]->Draw("PSAME");
 333     }
 334     cEffcfPos->BuildLegend();
 335
 336     TCanvas *cEffcfNeg=new TCanvas("cEffcfNeg","cEffcfNeg");
 337     cEffcfNeg->SetGridy();
 338     for(Int_t ipart=0;ipart<=2;ipart++){
 339       if(ipart==0)hEffNeg[ipart]->Draw("P");
 340       else hEffNeg[ipart]->Draw("PSAME");
 341       hEffNegMean[ipart]->Draw("PSAME");
 342     }
 343     cEffcfNeg->BuildLegend();
 344   }
 345
 346
 347
 348
 349   for(Int_t ipart=0;ipart<=2;ipart++) {
 350     hHistPosNSigma[ipart]->Scale(normMC,"width");
 351     hHistNegNSigma[ipart]->Scale(normMC,"width");
 352     hHistPosNSigmaMean[ipart]->Scale(normMC,"width");
 353     hHistNegNSigmaMean[ipart]->Scale(normMC,"width");
 354   }
 355
 356
 357   //////////////////////////////Prim/All Sec/All SecStr/All
 358   TH1F *hPrimPos[3];
 359   TH1F *hPrimNeg[3];
 360   TH1F *hSecPos[3];
 361   TH1F *hSecNeg[3];
 362   TH1F *hSecStrPos[3];
 363   TH1F *hSecStrNeg[3];
 364   TH1F *hAllPos[3];
 365   TH1F *hAllNeg[3];
 366   TH1F *hAllSecPos[3];
 367   TH1F *hAllSecNeg[3];
 368
 369   for(Int_t ipart=0;ipart<=2;ipart++){
 370     hPrimPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCposReco%i",ipart));
 371     hPrimNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCnegReco%i",ipart));
 372     hSecPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecMatMCposReco%i",ipart));
 373     hSecNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecMatMCnegReco%i",ipart));
 374     hSecStrPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecStrMCposReco%i",ipart));
 375     hSecStrNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecStrMCnegReco%i",ipart));
 376
 377     hAllSecPos[ipart]=(TH1F*)hSecPos[ipart]->Clone(Form("hAllSecPos%i",ipart));
 378     hAllSecPos[ipart]->Add(hSecStrPos[ipart]);
 379     hAllSecNeg[ipart]=(TH1F*)hSecNeg[ipart]->Clone(Form("hAllSecNeg%i",ipart));
 380     hAllSecNeg[ipart]->Add(hSecStrNeg[ipart]);
 381     hAllPos[ipart]=(TH1F*)hPrimPos[ipart]->Clone(Form("hAllPos%i",ipart));
 382     hAllPos[ipart]->Add(hSecPos[ipart]);
 383     hAllPos[ipart]->Add(hSecStrPos[ipart]);
 384     hAllNeg[ipart]=(TH1F*)hPrimNeg[ipart]->Clone(Form("hAllNeg%i",ipart));
 385     hAllNeg[ipart]->Add(hSecNeg[ipart]);
 386     hAllNeg[ipart]->Add(hSecStrNeg[ipart]);
 387
 388     hPrimPos[ipart]->Divide(hAllPos[ipart]);
 389     hSecPos[ipart]->Divide(hAllPos[ipart]);
 390     hSecStrPos[ipart]->Divide(hAllPos[ipart]);
 391     hPrimNeg[ipart]->Divide(hAllNeg[ipart]);
 392     hSecNeg[ipart]->Divide(hAllNeg[ipart]);
 393     hSecStrNeg[ipart]->Divide(hAllNeg[ipart]);
 394
 395     hAllSecPos[ipart]->Divide(hAllPos[ipart]);
 396     hAllSecNeg[ipart]->Divide(hAllNeg[ipart]);
 397     for(Int_t binnum=1;binnum<=hAllSecNeg[ipart]->GetNbinsX();binnum++){
 398       hAllSecPos[ipart]->SetBinContent(binnum,1-hAllSecPos[ipart]->GetBinContent(binnum));
 399       hAllSecNeg[ipart]->SetBinContent(binnum,1-hAllSecNeg[ipart]->GetBinContent(binnum));
 400     }
 401   }
 402
 403   /////////////////////////////////////////////
 404   TFile *outMC2=new TFile(foutMC.Data(),"update");
 405   TList *lMC=new TList();
 406   lMC->SetOwner();
 407   lMC->SetName(Form("MC_Mult%ito%i",lowmult,upmult));
 408   for(Int_t ipart=0;ipart<=2;ipart++){
 409     lMC->Add(hCorrFactPos[ipart]);
 410     lMC->Add(hCorrFactNeg[ipart]);
 411     lMC->Add(hCorrFactPosMean[ipart]);
 412     lMC->Add(hCorrFactNegMean[ipart]);
 413     lMC->Add(hEffPos[ipart]);
 414     lMC->Add(hEffNeg[ipart]);
 415     lMC->Add(hEffPosMean[ipart]);
 416     lMC->Add(hEffNegMean[ipart]);
 417     lMC->Add(hHistPosNSigma[ipart]);
 418     lMC->Add(hHistNegNSigma[ipart]);
 419     lMC->Add(hHistPosNSigmaMean[ipart]);
 420     lMC->Add(hHistNegNSigmaMean[ipart]);
 421     lMC->Add(hPrimPos[ipart]);
 422     lMC->Add(hPrimNeg[ipart]);
 423     lMC->Add(hSecPos[ipart]);
 424     lMC->Add(hSecNeg[ipart]);
 425     lMC->Add(hSecStrPos[ipart]);
 426     lMC->Add(hSecStrNeg[ipart]);
 427     lMC->Add(hAllSecPos[ipart]);
 428     lMC->Add(hAllSecNeg[ipart]);
 429   }
 430   lMC->Write(Form("MC_Mult%ito%i",lowmult,upmult),1);
 431   outMC2->Close();
 432   delete outMC2;
 433
 434
 435   TString lnameDATA=Form("clistITSsaMult%ito%i",lowmult,upmult);
 436   Printf("\n\n------------------ READING DATA  %s -------------------\n",lnameDATA.Data());
 437   TFile *finDATA = new TFile(fnameDATA.Data());
 438   TDirectory *dirFileDATA=(TDirectory*)finDATA->Get("PWG2SpectraITSsa");
 439   TList *cOutputDATA = (TList*)dirFileDATA->Get(lnameDATA.Data());
 440
 441
 442   for(Int_t ipart=0;ipart<=2;ipart++){
 443     hHistPosNSigmaDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistPosNSigma%i",ipart));
 444     hHistNegNSigmaDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistNegNSigma%i",ipart));
 445     hHistPosNSigmaMeanDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistPosNSigmaMean%i",ipart));
 446     hHistNegNSigmaMeanDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistNegNSigmaMean%i",ipart));
 447
 448     hITSsaPos[ipart]=(TH1F*)hHistPosNSigmaDATA[ipart]->Clone(Form("hITSsaPos%i",ipart));
 449     hITSsaNeg[ipart]=(TH1F*)hHistNegNSigmaDATA[ipart]->Clone(Form("hITSsaNeg%i",ipart));
 450     hITSsaPosMean[ipart]=(TH1F*)hHistPosNSigmaMeanDATA[ipart]->Clone(Form("hITSsaPosAsymm%i",ipart));
 451     hITSsaNegMean[ipart]=(TH1F*)hHistNegNSigmaMeanDATA[ipart]->Clone(Form("hITSsaNegAsymm%i",ipart));
 452     hITSsaRawPos[ipart]=(TH1F*)hHistPosNSigmaDATA[ipart]->Clone(Form("hITSsaRawPos%i",ipart));
 453     hITSsaRawNeg[ipart]=(TH1F*)hHistNegNSigmaDATA[ipart]->Clone(Form("hITSsaRawNeg%i",ipart));
 454     hITSsaRawPosMean[ipart]=(TH1F*)hHistPosNSigmaMeanDATA[ipart]->Clone(Form("hITSsaRawPosAsymm%i",ipart));
 455     hITSsaRawNegMean[ipart]=(TH1F*)hHistNegNSigmaMeanDATA[ipart]->Clone(Form("hITSsaRawNegAsymm%i",ipart));
 456
 457     hITSsaPos[ipart]->Sumw2();
 458     hITSsaNeg[ipart]->Sumw2();
 459     hITSsaPosMean[ipart]->Sumw2();
 460     hITSsaNegMean[ipart]->Sumw2();
 461     hITSsaRawPos[ipart]->Sumw2();
 462     hITSsaRawNeg[ipart]->Sumw2();
 463     hITSsaRawPosMean[ipart]->Sumw2();
 464     hITSsaRawNegMean[ipart]->Sumw2();
 465
 466     hITSsaPos[ipart]->Divide(hCorrFactPos[ipart]);
 467     hITSsaNeg[ipart]->Divide(hCorrFactNeg[ipart]);
 468     hITSsaPosMean[ipart]->Divide(hCorrFactPosMean[ipart]);
 469     hITSsaNegMean[ipart]->Divide(hCorrFactNegMean[ipart]);
 470
 471     SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hITSsaPos[ipart]);
 472     SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hITSsaNeg[ipart]);
 473     SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hITSsaPosMean[ipart]);
 474     SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hITSsaNegMean[ipart]);
 475     SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hITSsaRawPos[ipart]);
 476     SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hITSsaRawNeg[ipart]);
 477     SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hITSsaRawPosMean[ipart]);
 478     SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hITSsaRawNegMean[ipart]);
 479   }
 480
 481   if(!CorrectOnlyForEfficiency){
 482
 483
 484     if(DCACorrection){
 485       printf("- Correction from DCA\n");
 486
 487       // Correction factor based on fit to DCA distr for P and Pibar
 488       TString fnameDCAPosP=Form("./DCACorr%s_%s_%.0fDCA_PosP_TFraction.root",period.Data(),MCname.Data(),DCAcut);
 489       TFile *fDCAPosP= new TFile(fnameDCAPosP.Data());
 490       TH1F * hDCAPosP= (TH1F*)fDCAPosP->Get("hPrimAllDATAMC");
 491       //TH1F * hDCAPosP= (TH1F*)fDCAPosP->Get("hPrimAllDATA");
 492       TString fnameDCANegP=Form("./DCACorr%s_%s_%.0fDCA_NegP_TFraction.root",period.Data(),MCname.Data(),DCAcut);
 493       TFile *fDCANegP= new TFile(fnameDCANegP.Data());
 494       TH1F * hDCANegP= (TH1F*)fDCANegP->Get("hPrimAllDATAMC");
 495       //TH1F * hDCANegP= (TH1F*)fDCANegP->Get("hPrimAllDATA");
 496
 497       for(Int_t pbin=0;pbin<=hDCAPosP->GetNbinsX();pbin++)hDCAPosP->SetBinError(pbin,0);
 498       for(Int_t pbin=0;pbin<=hDCANegP->GetNbinsX();pbin++)hDCANegP->SetBinError(pbin,0);
 499
 500       hITSsaPos[2]->Multiply(hDCAPosP);
 501       hITSsaNeg[2]->Multiply(hDCANegP);
 502       hITSsaPosMean[2]->Multiply(hDCAPosP);
 503       hITSsaNegMean[2]->Multiply(hDCANegP);
 504
 505       hDCAPosP->Draw("");
 506       hDCANegP->Draw("same");
 507
 508     }
 509     else{
 510       Printf("NO DCA correction for protons");
 511     }
 512
 513     if(GFCorrection){
 514       printf("- Geant/FLUKA correction\n");
 515       //GeantFluka correction for Pi and K
 516       TString fnameGeanFlukaPi="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSection.211.root";
 517       TFile *fGeanFlukaPi= new TFile(fnameGeanFlukaPi.Data());
 518       TH1F *hGeantFlukaPiPos=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionParticles");
 519       TH1F *hGeantFlukaPiNeg=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionAntiParticles");
 520       for(Int_t binPi=0;binPi<=hITSsaPos[0]->GetNbinsX();binPi++){
 521         Float_t FlukaCorrPiPos=hGeantFlukaPiPos->GetBinContent(hGeantFlukaPiPos->FindBin(hITSsaPos[0]->GetBinCenter(binPi)));
 522         Float_t FlukaCorrPiNeg=hGeantFlukaPiNeg->GetBinContent(hGeantFlukaPiNeg->FindBin(hITSsaNeg[0]->GetBinCenter(binPi)));
 523         hITSsaPos[0]->SetBinContent(binPi,hITSsaPos[0]->GetBinContent(binPi)*FlukaCorrPiPos);
 524         hITSsaPos[0]->SetBinError(binPi,hITSsaPos[0]->GetBinError(binPi)*FlukaCorrPiPos);
 525         hITSsaNeg[0]->SetBinContent(binPi,hITSsaNeg[0]->GetBinContent(binPi)*FlukaCorrPiNeg);
 526         hITSsaNeg[0]->SetBinError(binPi,hITSsaNeg[0]->GetBinError(binPi)*FlukaCorrPiNeg);
 527         hITSsaPosMean[0]->SetBinContent(binPi,hITSsaPosMean[0]->GetBinContent(binPi)*FlukaCorrPiPos);
 528         hITSsaPosMean[0]->SetBinError(binPi,hITSsaPosMean[0]->GetBinError(binPi)*FlukaCorrPiPos);
 529         hITSsaNegMean[0]->SetBinContent(binPi,hITSsaNegMean[0]->GetBinContent(binPi)*FlukaCorrPiNeg);
 530         hITSsaNegMean[0]->SetBinError(binPi,hITSsaNegMean[0]->GetBinError(binPi)*FlukaCorrPiNeg);
 531       }
 532
 533       TString fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSection.321.root";
 534       TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
 535       TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
 536       TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
 537       for(Int_t binK=0;binK<=hITSsaPos[1]->GetNbinsX();binK++){
 538         Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(hITSsaPos[1]->GetBinCenter(binK)));
 539         Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(hITSsaNeg[1]->GetBinCenter(binK)));
 540         hITSsaPos[1]->SetBinContent(binK,hITSsaPos[1]->GetBinContent(binK)*FlukaCorrKPos);
 541         hITSsaNeg[1]->SetBinContent(binK,hITSsaNeg[1]->GetBinContent(binK)*FlukaCorrKNeg);
 542         hITSsaPos[1]->SetBinError(binK,hITSsaPos[1]->GetBinError(binK)*FlukaCorrKPos);
 543         hITSsaNeg[1]->SetBinError(binK,hITSsaNeg[1]->GetBinError(binK)*FlukaCorrKNeg);
 544         hITSsaPosMean[1]->SetBinContent(binK,hITSsaPosMean[1]->GetBinContent(binK)*FlukaCorrKPos);
 545         hITSsaNegMean[1]->SetBinContent(binK,hITSsaNegMean[1]->GetBinContent(binK)*FlukaCorrKNeg);
 546         hITSsaPosMean[1]->SetBinError(binK,hITSsaPosMean[1]->GetBinError(binK)*FlukaCorrKPos);
 547         hITSsaNegMean[1]->SetBinError(binK,hITSsaNegMean[1]->GetBinError(binK)*FlukaCorrKNeg);
 548       }
 549
 550       //Geant Fluka for P
 551       // ITS specific file for protons/antiprotons
 552       Int_t kNCharge=2;
 553       Int_t kPos=0;
 554       Int_t kNeg=1;
 555       TFile* fITS = new TFile ("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSectionITS_20100719.root");
 556       // TH2D * hCorrFlukaITS[kNCharge];
 557       TH2D * hCorrFlukaITS[2];
 558       hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
 559       hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
 560       for(Int_t icharge = 0; icharge < kNCharge; icharge++){
 561         Int_t nbins = hITSsaPos[2]->GetNbinsX();
 562         Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
 563         //for(Int_t binP=0;binP<=hITSsaPos[2]->GetNbinsX();binP++){
 564         for(Int_t ibin = 0; ibin < nbins; ibin++){
 565           Float_t pt = hITSsaPos[2]->GetBinCenter(ibin);
 566           Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
 567           Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
 568           if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
 569           if (pt > maxPtCorrection) pt = maxPtCorrection;
 570           Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
 571           cout<<correction<<"     charge "<<icharge<<endl;
 572           if(icharge==0){
 573             if (correction != 0) {// If the bin is empty this is a  0
 574               hITSsaPos[2]->SetBinContent(ibin,hITSsaPos[2]->GetBinContent(ibin)*correction);
 575               hITSsaPos[2]->SetBinError(ibin,hITSsaPos[2]->GetBinError  (ibin)*correction);
 576               hITSsaPosMean[2]->SetBinContent(ibin,hITSsaPosMean[2]->GetBinContent(ibin)*correction);
 577               hITSsaPosMean[2]->SetBinError(ibin,hITSsaPosMean[2]->GetBinError  (ibin)*correction);
 578
 579             }else if (hITSsaPos[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 580               cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << endl;
 581               cout << " Bin content: " << hITSsaPos[2]->GetBinContent(ibin)  << endl;
 582             }
 583           }
 584           if(icharge==1){
 585             if (correction != 0) {// If the bin is empty this is a  0
 586               hITSsaNeg[2]->SetBinContent(ibin,hITSsaNeg[2]->GetBinContent(ibin)*correction);
 587               hITSsaNeg[2]->SetBinError(ibin,hITSsaNeg[2]->GetBinError  (ibin)*correction);
 588               hITSsaNegMean[2]->SetBinContent(ibin,hITSsaNegMean[2]->GetBinContent(ibin)*correction);
 589               hITSsaNegMean[2]->SetBinError(ibin,hITSsaNegMean[2]->GetBinError  (ibin)*correction);
 590             }else if (hITSsaNeg[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
 591               cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, ITS, " << endl;
 592               cout << " Bin content: " << hITSsaNeg[2]->GetBinContent(ibin)  << endl;
 593             }
 594           }
 595
 596         }
 597       }
 598     }else{
 599       Printf("SPECTRA ARE NOT CORRECTED FOR GEANT/FLUKA");
 600
 601     }
 602
 603   }else{
 604     Printf("SPECTRA ARE NOT CORRECTED FOR DCA FIT + GEANT/FLUKA");
 605   }
 606
 607   Float_t norm;
 608   TH1F *fHistNEvents =0x0;
 609   if(MarekNorm){
 610     Printf("Normalization to the number of events with good vertex");
 611     TFile *finMAREK = new TFile(fnameDATA.Data());
 612     TDirectory *dirFileMAREK=(TDirectory*)finMAREK->Get("PWG2SpectraITSTPC");
 613     if(!dirFileMAREK)Printf("!dirFileMAREK");
 614     TString lnameMAREK=Form("outputlow%dup%dHI0",lowmult,upmult);
 615     TList *cOutputMAREK = (TList*)dirFileMAREK->Get(lnameMAREK.Data());
 616     if(!cOutputMAREK)Printf("!cOutputMAREK");
 617     TH1F *StatsHistMAREK=(TH1F*)cOutputMAREK->FindObject("StatsHist");
 618     norm=(Float_t)1/StatsHistMAREK->GetBinContent(3);
 619   }
 620   else{
 621     Printf("Normalization to the number of events after phys sel");
 622     fHistNEvents = (TH1F*)cOutputDATA->FindObject("fHistNEvents");
 623     norm=(Float_t)1/fHistNEvents->GetBinContent(fHistNEvents->FindBin(1.));
 624   }
 625   if(NormalizeData)printf("- Normalization (NEvents): %.0f \n",1/norm);
 626   else printf("- Normalization (NEvents): %.0f !!!!!!!SPECTRA ARE NOT NORMALIZED!!!!!!\n",fHistNEvents->GetBinContent(fHistNEvents->FindBin(1.)));
 627   for(Int_t ipart=0;ipart<=2;ipart++){
 628     if(NormalizeData){
 629       hITSsaPos[ipart]->Scale(norm,"width");
 630       hITSsaNeg[ipart]->Scale(norm,"width");
 631       hITSsaPosMean[ipart]->Scale(norm,"width");
 632       hITSsaNegMean[ipart]->Scale(norm,"width");
 633       hITSsaRawPos[ipart]->Scale(norm,"width");
 634       hITSsaRawNeg[ipart]->Scale(norm,"width");
 635       hITSsaRawPosMean[ipart]->Scale(norm,"width");
 636       hITSsaRawNegMean[ipart]->Scale(norm,"width");
 637     }
 638     else{
 639       hITSsaPos[ipart]->Scale(1,"width");
 640       hITSsaNeg[ipart]->Scale(1,"width");
 641       hITSsaPosMean[ipart]->Scale(1,"width");
 642       hITSsaNegMean[ipart]->Scale(1,"width");
 643       hITSsaRawPos[ipart]->Scale(1,"width");
 644       hITSsaRawNeg[ipart]->Scale(1,"width");
 645       hITSsaRawPosMean[ipart]->Scale(1,"width");
 646       hITSsaRawNegMean[ipart]->Scale(1,"width");
 647     }
 648     if(NormalizeToINEL){
 649       Float_t INEL=0.86;
 650       printf("Normalization to INEL : %.2f \n",INEL);
 651       hITSsaPos[ipart]->Scale(INEL,"");
 652       hITSsaNeg[ipart]->Scale(INEL,"");
 653       hITSsaPosMean[ipart]->Scale(INEL,"");
 654       hITSsaNegMean[ipart]->Scale(INEL,"");
 655       hITSsaRawPos[ipart]->Scale(INEL,"");
 656       hITSsaRawNeg[ipart]->Scale(INEL,"");
 657       hITSsaRawPosMean[ipart]->Scale(INEL,"");
 658       hITSsaRawNegMean[ipart]->Scale(INEL,"");
 659     }
 660     else printf("Normalization to INEL ignored \n");
 661
 662     hITSsaPos[ipart]->SetMinimum(0.00001);
 663     hITSsaNeg[ipart]->SetMinimum(0.00001);
 664     hITSsaPosMean[ipart]->SetMinimum(0.00001);
 665     hITSsaNegMean[ipart]->SetMinimum(0.00001);
 666     hITSsaRawPos[ipart]->SetMinimum(0.00001);
 667     hITSsaRawNeg[ipart]->SetMinimum(0.00001);
 668     hITSsaRawPosMean[ipart]->SetMinimum(0.00001);
 669     hITSsaRawNegMean[ipart]->SetMinimum(0.00001);
 670     hITSsaPos[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
 671     hITSsaNeg[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
 672     hITSsaPosMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
 673     hITSsaNegMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
 674     hITSsaRawPos[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
 675     hITSsaRawNeg[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
 676     hITSsaRawPosMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
 677     hITSsaRawNegMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
 678     hITSsaPos[ipart]->SetTitle(Form("ITSsa NSigma Symm Pos%i",ipart));
 679     hITSsaNeg[ipart]->SetTitle(Form("ITSsa NSigma Symm Neg%i",ipart));
 680     hITSsaPosMean[ipart]->SetTitle(Form("ITSsa NSigma Asymm Pos%i",ipart));
 681     hITSsaNegMean[ipart]->SetTitle(Form("ITSsa NSigma Asymm Neg%i",ipart));
 682     hITSsaRawPos[ipart]->SetTitle(Form("ITSsaRaw NSigma Symm Pos%i",ipart));
 683     hITSsaRawNeg[ipart]->SetTitle(Form("ITSsaRaw NSigma Symm Neg%i",ipart));
 684     hITSsaRawPosMean[ipart]->SetTitle(Form("ITSsaRaw NSigma Asymm Pos%i",ipart));
 685     hITSsaRawNegMean[ipart]->SetTitle(Form("ITSsaRaw NSigma Asymm Neg%i",ipart));
 686
 687
 688     ResetOutOfRange(hITSsaPos[ipart],ipart,lowRange,upRange);
 689     ResetOutOfRange(hITSsaNeg[ipart],ipart,lowRange,upRange);
 690     ResetOutOfRange(hITSsaPosMean[ipart],ipart,lowRange,upRange);
 691     ResetOutOfRange(hITSsaNegMean[ipart],ipart,lowRange,upRange);
 692
 693   }
 694   TH1F *hSystematicsITSsaPos[3];
 695   TH1F *hSystematicsITSsaNeg[3];
 696   TH1F *hSystematicsITSsaPosMean[3];
 697   TH1F *hSystematicsITSsaNegMean[3];
 698
 699   if(MakeSystErr){
 700     printf("- making Systematic Error\n");
 701     //preliminary systematic error
 702     TFile *fsys = new TFile("RootFilesSystError/ITSsa-systematics-20101014.root");
 703     //fsys->Print("all");
 704     //cout<<fsys->GetSize()<<endl;
 705     TH1F* hSystPos[3];
 706     TH1F* hSystNeg[3];
 707     hSystPos[0]=(TH1F*)fsys->Get("hSystTotPosPion");
 708     hSystNeg[0]=(TH1F*)fsys->Get("hSystTotNegPion");
 709     hSystPos[1]=(TH1F*)fsys->Get("hSystTotPosKaon");
 710     hSystNeg[1]=(TH1F*)fsys->Get("hSystTotNegKaon");
 711     hSystPos[2]=(TH1F*)fsys->Get("hSystTotPosProton");
 712     hSystNeg[2]=(TH1F*)fsys->Get("hSystTotNegProton");
 713
 714     for(Int_t ipart=0;ipart<=2;ipart++){
 715       hSystematicsITSsaPos[ipart]=(TH1F*)hITSsaPos[ipart]->Clone(Form("hSystematicsITSsaPos%i",ipart));
 716       hSystematicsITSsaNeg[ipart]=(TH1F*)hITSsaNeg[ipart]->Clone(Form("hSystematicsITSsaNeg%i",ipart));
 717       hSystematicsITSsaPosMean[ipart]=(TH1F*)hITSsaPosMean[ipart]->Clone(Form("hSystematicsITSsaPosAsymm%i",ipart));
 718       hSystematicsITSsaNegMean[ipart]=(TH1F*)hITSsaNegMean[ipart]->Clone(Form("hSystematicsITSsaNegAsymm%i",ipart));
 719
 720       for(Int_t ibin=0;ibin<=hSystematicsITSsaPos[ipart]->GetNbinsX();ibin++){
 721         Float_t syserr=hSystPos[ipart]->GetBinContent(hSystPos[ipart]->FindBin(hSystematicsITSsaPos[ipart]->GetBinCenter(ibin)));
 722         hSystematicsITSsaPos[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaPos[ipart]->GetBinContent(ibin));
 723         hSystematicsITSsaPosMean[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaPosMean[ipart]->GetBinContent(ibin));
 724       }
 725       for(Int_t ibin=0;ibin<=hSystematicsITSsaNeg[ipart]->GetNbinsX();ibin++){
 726         Float_t syserr=hSystNeg[ipart]->GetBinContent(hSystNeg[ipart]->FindBin(hSystematicsITSsaNeg[ipart]->GetBinCenter(ibin)));
 727         hSystematicsITSsaNeg[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaNeg[ipart]->GetBinContent(ibin));
 728         hSystematicsITSsaNegMean[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaNegMean[ipart]->GetBinContent(ibin));
 729       }
 730     }
 731   }
 732
 733   TCanvas *cSpectraSymm=new TCanvas("cSpectraSymm","cSpectraSymm");
 734   cSpectraSymm->SetGridy();
 735   for(Int_t ipart=0;ipart<=2;ipart++){
 736     if(ipart==0)hITSsaPos[ipart]->DrawCopy("P");
 737     else hITSsaPos[ipart]->DrawCopy("PSAME");
 738     hITSsaNeg[ipart]->DrawCopy("PSAME");
 739   }
 740   cSpectraSymm->BuildLegend();
 741
 742   TCanvas *cSpectraAsymm=new TCanvas("cSpectraAsymm","cSpectraAsymm");
 743   cSpectraAsymm->SetGridy();
 744   for(Int_t ipart=0;ipart<=2;ipart++){
 745     if(ipart==0)hITSsaPosMean[ipart]->DrawCopy("P");
 746     else hITSsaPosMean[ipart]->DrawCopy("PSAME");
 747     hITSsaNegMean[ipart]->DrawCopy("PSAME");
 748   }
 749   cSpectraAsymm->BuildLegend();
 750
 751
 752
 753   ///////////////////////////  SCALING SPECTRA FOR TRACKING EFFICIENCY   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 754   Float_t ScalingFactor=1.015;
 755   Printf("SCALING FACTOR: %f",ScalingFactor);
 756   for(Int_t ipart=0;ipart<=2;ipart++){
 757     hITSsaPosMean[ipart]->Scale(ScalingFactor);
 758     Printf("Scaling %s by %f",hITSsaPosMean[ipart]->GetName(),ScalingFactor);
 759     hITSsaNegMean[ipart]->Scale(ScalingFactor);
 760     Printf("Scaling %s by %f",hITSsaNegMean[ipart]->GetName(),ScalingFactor);
 761     hITSsaPos[ipart]->Scale(ScalingFactor);
 762     Printf("Scaling %s by %f",hITSsaPos[ipart]->GetName(),ScalingFactor);
 763     hITSsaNeg[ipart]->Scale(ScalingFactor);
 764     Printf("Scaling %s by %f",hITSsaNeg[ipart]->GetName(),ScalingFactor);
 765   }
 766
 767   TFile *outDATA2=new TFile(foutDATA.Data(),"update");
 768   TList *lDATA=new TList();
 769   lDATA->SetOwner();
 770   lDATA->SetName(Form("DATA_Mult%ito%i",lowmult,upmult));
 771   if(SaveOnlyAsymm){
 772     for(Int_t ipart=0;ipart<=2;ipart++)  lDATA->Add(hITSsaPosMean[ipart]);
 773     for(Int_t ipart=0;ipart<=2;ipart++)  lDATA->Add(hITSsaNegMean[ipart]);
 774   }else{
 775     for(Int_t ipart=0;ipart<=2;ipart++){
 776       lDATA->Add(fHistPrimMCposBefEvSel[ipart]);
 777       lDATA->Add(fHistPrimMCnegBefEvSel[ipart]);
 778       lDATA->Add(hHistPosNSigmaPrimMean[ipart]);
 779       lDATA->Add(hHistNegNSigmaPrimMean[ipart]);
 780       lDATA->Add(hHistPosNSigmaPrimMCMean[ipart]);
 781       lDATA->Add(hHistNegNSigmaPrimMCMean[ipart]);
 782       lDATA->Add(hITSsaPos[ipart]);
 783       lDATA->Add(hITSsaNeg[ipart]);
 784       lDATA->Add(hITSsaPosMean[ipart]);
 785       lDATA->Add(hITSsaNegMean[ipart]);
 786       lDATA->Add(hITSsaRawPos[ipart]);
 787       lDATA->Add(hITSsaRawNeg[ipart]);
 788       lDATA->Add(hITSsaRawPosMean[ipart]);
 789       lDATA->Add(hITSsaRawNegMean[ipart]);
 790     }
 791     if(MakeSystErr){
 792       for(Int_t ipart=0;ipart<=2;ipart++){
 793         lDATA->Add(hSystematicsITSsaPos[ipart]);
 794         lDATA->Add(hSystematicsITSsaNeg[ipart]);
 795         lDATA->Add( hSystematicsITSsaPosMean[ipart]);
 796         lDATA->Add( hSystematicsITSsaNegMean[ipart]);
 797       }
 798     }
 799   }
 800   lDATA->Write(Form("DATA_Mult%ito%i",lowmult,upmult),1);
 801   outDATA2->Close();
 802   delete outDATA2;
 803
 804
 805   RatioPlot(hITSsaNeg,hITSsaPos,"Neg","Pos","Symm NSigma",Low,Up);
 806   RatioPlot(hITSsaNegMean,hITSsaPosMean,"Neg","Pos","Assymm NSigma",Low,Up);
 807   RatioPlot(hITSsaPos,hITSsaPosMean,"Symm","Asymm","Positive",Low,Up);
 808   RatioPlot(hITSsaNeg,hITSsaNegMean,"Symm","Asymm","Negative",Low,Up);
 809
 810
 811 } //end main
 812
 813
 814 void ResetOutOfRange(TH1F *histo, Int_t ipart, Double_t lowRange[3], Double_t upRange[3]){
 815   // set to -1 the bin contents out of the selcted pt range
 816   for(Int_t ibin=0;ibin<=histo->GetNbinsX();ibin++){
 817     if(histo->GetBinCenter(ibin)<lowRange[ipart] || histo->GetBinCenter(ibin)>upRange[ipart]) histo->SetBinContent(ibin,-1);
 818   }
 819 }
 820
 821 void RatioPlot(TH1F **num,TH1F **den,TString t1,TString t2,TString opt,Double_t Low[],Double_t Up[])
 822 {
 823   gROOT->SetStyle("Plain");
 824   gStyle->SetOptTitle(0);
 825   TH1F *Num[3];
 826   TString title=Form("%s/%s - %s",t1.Data(),t2.Data(),opt.Data());
 827   TCanvas *c=new TCanvas(title.Data(),title.Data());
 828   c->Divide(3,2);
 829   for(Int_t ipart=0;ipart<=2;ipart++){
 830     num[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
 831     num[ipart]->SetMarkerColor(2);
 832     num[ipart]->SetLineColor(2);
 833     num[ipart]->SetStats(kFALSE);
 834     den[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
 835     den[ipart]->SetMarkerColor(1);
 836     den[ipart]->SetLineColor(1);
 837     den[ipart]->SetStats(kFALSE);
 838     c->cd(ipart+1);
 839     gPad->SetGridy();
 840     //gPad->SetLogy();
 841     num[ipart]->Draw("");
 842     den[ipart]->Draw("same");
 843     TPaveText *tpave=new TPaveText(0.1,0.2,0.5,0.29,"brNDC");
 844     tpave->SetBorderSize(0);
 845     tpave->SetFillStyle(0);
 846     tpave->SetFillColor(0);
 847     tpave->SetTextColor(2);
 848     tpave->SetTextSize(.04);
 849     TText *txt1=tpave->AddText(t1.Data());
 850     txt1->SetTextFont(62);
 851     txt1->SetTextColor(2);
 852     TText *txt2=tpave->AddText(t2.Data());
 853     txt2->SetTextFont(62);
 854     txt2->SetTextColor(1);
 855     tpave->Draw();
 856     Num[ipart]=(TH1F*)num[ipart]->Clone(title.Data());
 857     Num[ipart]->SetTitle(title.Data());
 858     Num[ipart]->Divide(den[ipart]);
 859     Num[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
 860     Num[ipart]->SetMarkerColor(4);
 861     Num[ipart]->SetLineColor(4);
 862     Num[ipart]->SetMaximum(1.2);
 863     Num[ipart]->SetMinimum(.8);
 864     c->cd(ipart+4);
 865     gPad->SetGridy();
 866     Num[ipart]->SetStats(kFALSE);
 867     Num[ipart]->Draw();
 868     TPaveText *tpave_1=new TPaveText(0.3,0.7,0.7,0.99,"brNDC");
 869     tpave_1->SetBorderSize(0);
 870     tpave_1->SetFillStyle(0);
 871     tpave_1->SetFillColor(0);
 872     tpave_1->SetTextColor(2);
 873     tpave_1->SetTextSize(.04);
 874     TText *txt1_1=tpave_1->AddText(title.Data());
 875     txt1_1->SetTextFont(62);
 876     tpave_1->Draw();
 877     // TLine line = TLine(Low[ipart],1,Up[ipart],1);
 878     // line.SetLineColor(2);
 879     // line.SetLineStyle(2);
 880     // line.SetLineWidth(2);
 881     // line.DrawClone("same");
 882   }
 883 }
 884
 885
 886
 887
 888  void SetDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1)
 889 {
 890   h1->SetMarkerStyle(markerstyle);
 891   h1->SetMarkerColor(markercolor);
 892   h1->SetMarkerSize(markersize);
 893   h1->SetLineColor(linecolor);
 894   h1->SetLineWidth(linewidth);
 895 }
 896
 897