PWGLF/SPECTRA/PiKaPr/ITSsa/MakeRawITSsaSpectraMultiBin.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include <TLatex.h>
   4 #include <TImage.h>
   5 #include <TSystem.h>
   6 #include <TPaveText.h>
   7 #include <TH1F.h>
   8 #include <TF1.h>
   9 #include <TH1D.h>
  10 #include <TH2F.h>
  11 #include <TMath.h>
  12 #include <TNtuple.h>
  13 #include <TGraphErrors.h>
  14 #include <TList.h>
  15 #include <TLegend.h>
  16 #include <TLegendEntry.h>
  17 #include <TCanvas.h>
  18 #include <TFile.h>
  19 #include <TStyle.h>
  20 #include <Rtypes.h>
  21 #include "AliITSsadEdxFitter.h"
  22 #endif
  23
  24
  25 Double_t BetheBloch(Double_t *mom, Double_t *mass);
  26 void Logo();
  27 //
  28
  29 //______________________________________________________________________
  30 void MakeRawITSsaSpectraMultiBin(Bool_t optMC=kTRUE, Int_t multibin=0){
  31
  32
  33   AliITSsadEdxFitter *ITSsa=new AliITSsadEdxFitter(optMC);
  34
  35   TString filename, dirName, ps0, ps1, ps2;
  36   Char_t listname[50];
  37   if(optMC) dirName=(Form("../gridmultiplicitybins/LHC10d1_1.5sigma_7DCA_negmag/Spectra_MC_negmag_MultBin%d",multibin));
  38   else dirName=(Form("../gridmultiplicitybins/data_1.5sigma_7DCA_negmag/Spectra_data_negmag_MultBin%d",multibin));
  39   switch(multibin){
  40   case 0:
  41     sprintf(listname,"clistITSsaMult0to9999");
  42     break;
  43   case 1:
  44     sprintf(listname,"clistITSsaMult0to5");
  45     break;
  46   case 2:
  47     sprintf(listname,"clistITSsaMult6to9");
  48     break;
  49   case 3:
  50     sprintf(listname,"clistITSsaMult10to14");
  51     break;
  52   case 4:
  53     sprintf(listname,"clistITSsaMult15to22");
  54     break;
  55   case 5:
  56     sprintf(listname,"clistITSsaMult23to9999");
  57     break;
  58   }
  59   if(optMC){
  60     filename=Form("%s/AnalysisResults.root",dirName.Data());
  61     ps0 = Form("%s/outSIM.ps[",dirName.Data());
  62     ps1 = Form("%s/outSIM.ps",dirName.Data());
  63     ps2 = Form("%s/outSIM.ps]",dirName.Data());
  64   }
  65   else{
  66     filename=Form("%s/AnalysisResults.root",dirName.Data());
  67     ps0 = Form("%s/outDATA.ps[",dirName.Data());
  68     ps1 = Form("%s/outDATA.ps",dirName.Data());
  69     ps2 = Form("%s/outDATA.ps]",dirName.Data());
  70   }
  71   TString openfilename="./";
  72   openfilename+=filename;
  73   cout<<openfilename<<endl;
  74   TFile *fi=new TFile(openfilename.Data());
  75   if(!fi){
  76     cout<<"TFile loading failed"<<endl;
  77     return;
  78   }
  79   TDirectoryFile *di=(TDirectoryFile*) fi->Get("PWG2SpectraITSsa");
  80   if(!di){
  81     cout<<"TDirectory loading failed!"<<endl;
  82     return;
  83   }
  84   TList *li=(TList*)di->Get(listname);
  85   if(!li){
  86     cout<<"TList loading failed"<<endl;
  87     return;
  88   }
  89   cout<<"File loaded"<<endl;
  90
  91   TCanvas *cdummy=new TCanvas("dummy","IntegralMethod",1000,800);
  92   cdummy->Print(ps0.Data());
  93
  94   //binning
  95   const Int_t nbins = 22;
  96   Double_t xbins[nbins+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};
  97
  98   //histograms
  99   TH1F *fHistMCPosPi[nbins];
 100   TH1F *fHistMCPosK[nbins];
 101   TH1F *fHistMCPosP[nbins];
 102   TH1F *fHistMCNegPi[nbins];
 103   TH1F *fHistMCNegK[nbins];
 104   TH1F *fHistMCNegP[nbins];
 105   TH1F *fHistPosPi[nbins];
 106   TH1F *fHistPosK[nbins];
 107   TH1F *fHistPosP[nbins];
 108   TH1F *fHistNegPi[nbins];
 109   TH1F *fHistNegK[nbins];
 110   TH1F *fHistNegP[nbins];
 111   for(Int_t m=0;m<nbins;m++){
 112     fHistMCNegPi[m] = (TH1F*)li->FindObject(Form("fHistMCNegPi%d",m));
 113     fHistMCNegK[m] = (TH1F*)li->FindObject(Form("fHistMCNegK%d",m));
 114     fHistMCNegP[m] = (TH1F*)li->FindObject(Form("fHistMCNegP%d",m));
 115     fHistMCPosPi[m] = (TH1F*)li->FindObject(Form("fHistMCPosPi%d",m));
 116     fHistMCPosK[m] = (TH1F*)li->FindObject(Form("fHistMCPosK%d",m));
 117     fHistMCPosP[m] = (TH1F*)li->FindObject(Form("fHistMCPosP%d",m));
 118     fHistPosPi[m] = (TH1F*)li->FindObject(Form("fHistPosPi%d",m));
 119     fHistPosK[m] = (TH1F*)li->FindObject(Form("fHistPosK%d",m));
 120     fHistPosP[m] = (TH1F*)li->FindObject(Form("fHistPosP%d",m));
 121     fHistNegPi[m] = (TH1F*)li->FindObject(Form("fHistNegPi%d",m));
 122     fHistNegK[m] = (TH1F*)li->FindObject(Form("fHistNegK%d",m));
 123     fHistNegP[m] = (TH1F*)li->FindObject(Form("fHistNegP%d",m));
 124   }
 125
 126   TH1F *fHistNEvents = (TH1F*)li->FindObject("fHistNEvents");
 127   TH2F *fHistDEDX = (TH2F*)li->FindObject("fHistDEDX");
 128   TH2F *fHistDEDXdouble = (TH2F*)li->FindObject("fHistDEDXdouble");
 129   TH1F *fHistCharge[4];
 130   for(Int_t j=0;j<4;j++) fHistCharge[j] = (TH1F*)li->FindObject((Form("fHistChargeLay%d",j)));
 131
 132   TH1F *hEffPos[3];
 133   TH1F *hEffNeg[3];
 134   TH1F *hCorrFacPos[3];
 135   TH1F *hCorrFacNeg[3];
 136   TH1F *hEffMCPIDPos[3];
 137   TH1F *hEffMCPIDNeg[3];
 138   TH1F *hCorrFacMCPIDNeg[3];
 139   TH1F *hCorrFacMCPIDPos[3];
 140   TH1F *hSpectraPrimPosMC[3];
 141   TH1F *hSpectraPrimNegMC[3];
 142   TH1F *hSpectraPrimPosMCBefEvSel[3];
 143   TH1F *hSpectraPrimNegMCBefEvSel[3];
 144   TH1F *hSpectraPos[3];
 145   TH1F *hSpectraNeg[3];
 146   TH1F *hSpectraMCPIDPos[3];
 147   TH1F *hSpectraMCPIDNeg[3];
 148   TH1F* hMeanPos[3];
 149   TH1F* hMeanNeg[3];
 150   TH1F* hSigmaPos[3];
 151   TH1F* hSigmaNeg[3];
 152   TGraph *gres[6][22];
 153
 154   for(Int_t i=0; i<3; i++){
 155     hSpectraPrimPosMC[i]=(TH1F*)li->FindObject(Form("fHistPrimMCpos%d",i));
 156     hSpectraPrimNegMC[i]=(TH1F*)li->FindObject(Form("fHistPrimMCneg%d",i));
 157     hSpectraPrimPosMCBefEvSel[i]=(TH1F*)li->FindObject(Form("fHistPrimMCposBefEvSel%d",i));
 158     hSpectraPrimNegMCBefEvSel[i]=(TH1F*)li->FindObject(Form("fHistPrimMCnegBefEvSel%d",i));
 159     hSpectraMCPIDPos[i]=new TH1F(Form("hSpectraMCPIDPos%d",i),Form("hSpectraMCPIDPos%d",i),nbins,xbins);
 160     hSpectraMCPIDNeg[i]=new TH1F(Form("hSpectraMCPIDNeg%d",i),Form("hSpectraMCPIDNeg%d",i),nbins,xbins);
 161     hSpectraPos[i]=new TH1F(Form("hSpectraPos%d",i),Form("hSpectraPos%d",i),nbins,xbins);
 162     hSpectraNeg[i]=new TH1F(Form("hSpectraNeg%d",i),Form("hSpectraNeg%d",i),nbins,xbins);
 163     hMeanPos[i]=new TH1F(Form("hMeanPos%d",i),Form("hMeanPos%d",i),nbins,xbins);
 164     hMeanNeg[i]=new TH1F(Form("hMeanNeg%d",i),Form("hMeanNeg%d",i),nbins,xbins);
 165     hSigmaPos[i]=new TH1F(Form("hSigmaPos%d",i),Form("hSigmaPos%d",i),nbins,xbins);
 166     hSigmaNeg[i]=new TH1F(Form("hSigmaNeg%d",i),Form("hSigmaNeg%d",i),nbins,xbins);
 167   }
 168
 169   //division for DeltaPt (for MC spectra)
 170   for(Int_t ipt=0;ipt<3;ipt++){
 171     for(Int_t bin=1; bin <= hSpectraPrimPosMC[ipt]->GetNbinsX(); bin++){
 172       Float_t binSize=hSpectraPrimPosMC[ipt]->GetBinLowEdge(bin+1) - hSpectraPrimPosMC[ipt]->GetBinLowEdge(bin);
 173       hSpectraPrimPosMC[ipt]->SetBinContent(bin, hSpectraPrimPosMC[ipt]->GetBinContent(bin) / binSize);
 174       hSpectraPrimPosMC[ipt]->SetBinError(bin, 0);
 175       hSpectraPrimNegMC[ipt]->SetBinContent(bin, hSpectraPrimNegMC[ipt]->GetBinContent(bin) / binSize);
 176       hSpectraPrimNegMC[ipt]->SetBinError(bin, 0);
 177       if(hSpectraPrimPosMCBefEvSel[ipt]){
 178         hSpectraPrimPosMCBefEvSel[ipt]->SetBinContent(bin, hSpectraPrimPosMCBefEvSel[ipt]->GetBinContent(bin) / binSize);
 179         hSpectraPrimPosMCBefEvSel[ipt]->SetBinError(bin, 0);
 180       }
 181       if(hSpectraPrimNegMCBefEvSel[ipt]){
 182         hSpectraPrimNegMCBefEvSel[ipt]->SetBinContent(bin, hSpectraPrimNegMCBefEvSel[ipt]->GetBinContent(bin) / binSize);
 183         hSpectraPrimNegMCBefEvSel[ipt]->SetBinError(bin, 0);
 184       }
 185     }
 186   }
 187   cout<<"All plots loaded"<<endl;
 188
 189   //open output file to store the dedx distribution and the spectra
 190   TString savename=filename.Data();
 191   savename.ReplaceAll("AnalysisResults","SpectraReco");
 192   TFile *fout=new TFile(savename.Data(),"recreate");
 193   fout->cd();
 194
 195   gStyle->SetOptStat(0);
 196   gStyle->SetOptFit(111);
 197   gStyle->SetPalette(1);
 198
 199   //propaganda plot
 200   Float_t pdgmass[5]={0.13957,0.493677,0.938272,1.875612762,0.00596}; //mass for pi, K, P, d (Gev/c^2)
 201   TF1 *funPos[5];
 202   TF1 *funNeg[5];
 203   for(Int_t m=0;m<5;m++){
 204     funPos[m] = new TF1(Form("funPos%d",m),BetheBloch,0.02,5,1);
 205     funPos[m]->SetParameter(0,pdgmass[m]);
 206     funPos[m]->SetLineWidth(2);
 207     funNeg[m] = new TF1(Form("funNeg%d",m),BetheBloch,-5,0.02,1);
 208     funNeg[m]->SetParameter(0,-pdgmass[m]);
 209     funNeg[m]->SetLineWidth(2);
 210   }
 211
 212   TCanvas *cEvents=new TCanvas("cEvents","cEvents",500,400);
 213   cEvents->cd();
 214   fHistNEvents->Draw("text");
 215   fHistNEvents->Write();
 216
 217   TCanvas *cDEDX=new TCanvas("cDEDX","DEDX",1000,800);
 218   cDEDX->cd();
 219   gPad->SetLogx();
 220   gPad->SetLogz();
 221   fHistDEDX->GetXaxis()->SetRangeUser(0.08,5);
 222   fHistDEDX->GetYaxis()->SetRangeUser(0.,700);
 223   fHistDEDX->GetXaxis()->SetTitle("momentum [GeV/c]");
 224   fHistDEDX->GetYaxis()->SetTitle("dE [keV/300#mum]");
 225   fHistDEDX->Draw("colz");
 226   fHistDEDX->Write();
 227   for(Int_t m=0;m<3;m++) funPos[m]->Draw("same");
 228   Logo();
 229
 230   TCanvas *cDEDXdouble=new TCanvas("cDEDXdouble","DEDXdouble",1000,800);
 231   cDEDXdouble->cd();
 232   gPad->SetLogz();
 233   fHistDEDXdouble->GetXaxis()->SetRangeUser(-3,3);
 234   fHistDEDXdouble->GetXaxis()->SetTitle("momentum * sign [GeV/c]");
 235   fHistDEDXdouble->GetYaxis()->SetTitle("dE [keV/300#mum]");
 236   fHistDEDXdouble->Draw("colz");
 237   fHistDEDXdouble->Write();
 238   for(Int_t m=0;m<3;m++) {
 239     funPos[m]->Draw("same");
 240     funNeg[m]->Draw("same");
 241   }
 242   Logo();
 243
 244   //calibration check histo
 245   TCanvas *cs=new TCanvas("cs","cs",1000,800);
 246   cs->Divide(2,2);
 247   for(Int_t j=0;j<4;j++){
 248     cs->cd(j+1);
 249     fHistCharge[j]->GetXaxis()->SetTitle("Charge [keV]");
 250     if(j==0) fHistCharge[j]->SetTitle("Drift inner layer");
 251     if(j==1) fHistCharge[j]->SetTitle("Drift outer layer");
 252     if(j==2) fHistCharge[j]->SetTitle("Strip inner layer");
 253     if(j==3) fHistCharge[j]->SetTitle("Strip inner layer");
 254     fHistCharge[j]->Draw();
 255     fHistCharge[j]->SetFillColor(7);
 256     fHistCharge[j]->Fit("gaus","QR","",70,100);
 257     fHistCharge[j]->Write();
 258   }
 259
 260   //canvas MC spectra
 261   if(optMC){
 262     TCanvas *cspectraMC=new TCanvas("cspectraMC","cspectraMC",1000,800);
 263     cspectraMC->Divide(2,1);
 264     cspectraMC->cd(1);
 265     gPad->SetLogy();
 266     gPad->SetGridy();
 267     for(Int_t i=0;i<3;i++){
 268       if(i==0){
 269         hSpectraPrimPosMC[i]->Draw("");
 270         hSpectraPrimPosMC[i]->SetTitle("ITS EvMC truth positive");
 271         hSpectraPrimPosMC[i]->SetMinimum(100);
 272         hSpectraPrimPosMC[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
 273         hSpectraPrimPosMC[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 274       }
 275       hSpectraPrimPosMC[i]->Draw("same");
 276       hSpectraPrimPosMC[i]->SetLineColor(i+2);
 277       hSpectraPrimPosMC[i]->SetMarkerColor(i+2);
 278       hSpectraPrimPosMC[i]->SetMarkerStyle(21+i);
 279       hSpectraPrimPosMC[i]->Write();
 280     }
 281     cspectraMC->cd(2);
 282     gPad->SetLogy();
 283     gPad->SetGridy();
 284     for(Int_t i=0;i<3;i++){
 285       if(i==0){
 286         hSpectraPrimNegMC[i]->Draw("");
 287         hSpectraPrimNegMC[i]->SetTitle("ITS MC truth negative");
 288         hSpectraPrimNegMC[i]->SetMinimum(100);
 289         hSpectraPrimNegMC[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
 290         hSpectraPrimNegMC[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 291       }
 292       hSpectraPrimNegMC[i]->Draw("same");
 293       hSpectraPrimNegMC[i]->SetLineColor(i+2);
 294       hSpectraPrimNegMC[i]->SetMarkerColor(i+2);
 295       hSpectraPrimNegMC[i]->SetMarkerStyle(21+i);
 296       hSpectraPrimNegMC[i]->Write();
 297     }
 298     cspectraMC->Print(ps0.Data());
 299   }
 300
 301   //dedx distribution to be fitted
 302   TCanvas *cgausPipos=new TCanvas("cgausPipos","PIONS pos",1000,800);
 303   cgausPipos->Divide(5,4,0.001,0.001);
 304   TCanvas *cgausKpos=new TCanvas("cgausKpos","KAONS pos",1000,800);
 305   cgausKpos->Divide(5,4,0.001,0.001);
 306   TCanvas *cgausPpos=new TCanvas("cgausPpos","PROTONS pos",1000,800);
 307   cgausPpos->Divide(5,4,0.001,0.001);
 308   TCanvas *cgausPineg=new TCanvas("cgausPineg","PIONS neg",1000,800);
 309   cgausPineg->Divide(5,4,0.001,0.001);
 310   TCanvas *cgausKneg=new TCanvas("cgausKneg","KAONS neg",1000,800);
 311   cgausKneg->Divide(5,4,0.001,0.001);
 312   TCanvas *cgausPneg=new TCanvas("cgausPneg","PROTONS neg",1000,800);
 313   cgausPneg->Divide(5,4,0.001,0.001);
 314
 315   //binsize for one histo only because are all equal
 316   Float_t binsize = fHistPosPi[0]->GetBinWidth(1);
 317   Double_t fpar[5],efpar[5];
 318   for(Int_t i=0; i<nbins-2; i++){
 319     //------------------- positive particles
 320     //pions
 321     cgausPipos->cd(i+1);
 322     gPad->SetLogy();
 323     fHistPosPi[i]->Write();
 324     gres[0][i]=new TGraph();
 325     ITSsa->DoFit(fHistPosPi[i],i,211,gres[0][i]);
 326     ITSsa->FillHisto(hSpectraPos[0],i+1,binsize,211);
 327     if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDPos[0],i+1,211,fHistMCPosPi[i]);
 328     if(ITSsa->IsGoodBin(i,211)){
 329       ITSsa->GetFitPar(fpar,efpar);
 330       hMeanPos[0]->SetBinContent(i+1,fpar[1]);
 331       hMeanPos[0]->SetBinError(i+1,efpar[1]);
 332       hSigmaPos[0]->SetBinContent(i+1,fpar[2]);
 333       hSigmaPos[0]->SetBinError(i+1,efpar[2]);
 334     }
 335     cgausPipos->Update();
 336     //kaons
 337     cgausKpos->cd(i+1);
 338     gPad->SetLogy();
 339     fHistPosK[i]->Write();
 340     gres[1][i]=new TGraph();
 341     ITSsa->DoFit(fHistPosK[i],i,321,gres[1][i]);
 342     ITSsa->FillHisto(hSpectraPos[1],i+1,binsize,321);
 343     if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDPos[1],i+1,321,fHistMCPosK[i]);
 344     if(ITSsa->IsGoodBin(i,321)){
 345       ITSsa->GetFitPar(fpar,efpar);
 346       hMeanPos[1]->SetBinContent(i+1,fpar[1]);
 347       hMeanPos[1]->SetBinError(i+1,efpar[1]);
 348       hSigmaPos[1]->SetBinContent(i+1,fpar[2]);
 349       hSigmaPos[1]->SetBinError(i+1,efpar[2]);
 350     }
 351     cgausKpos->Update();
 352     //protons
 353     cgausPpos->cd(i+1);
 354     gPad->SetLogy();
 355     fHistPosP[i]->Write();
 356     fHistPosP[i]->SetFillColor(16);
 357     gres[2][i]=new TGraph();
 358     ITSsa->DoFitProton(fHistPosP[i],i,2212,gres[2][i]);
 359     ITSsa->FillHisto(hSpectraPos[2],i+1,binsize,2212);
 360     if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDPos[2],i+1,2212,fHistMCPosP[i]);
 361     if(ITSsa->IsGoodBin(i,2212)){
 362       ITSsa->GetFitPar(fpar,efpar);
 363       hMeanPos[2]->SetBinContent(i+1,fpar[1]);
 364       hMeanPos[2]->SetBinError(i+1,efpar[1]);
 365       hSigmaPos[2]->SetBinContent(i+1,fpar[2]);
 366       hSigmaPos[2]->SetBinError(i+1,efpar[2]);
 367     }
 368     cgausPpos->Update();
 369
 370     //------------------- negative particles
 371     //pions
 372     cgausPineg->cd(i+1);
 373     gPad->SetLogy();
 374     fHistNegPi[i]->Write();
 375     gres[3][i]=new TGraph();
 376     ITSsa->DoFit(fHistNegPi[i],i,211,gres[3][i]);
 377     ITSsa->FillHisto(hSpectraNeg[0],i+1,binsize,211);
 378     if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDNeg[0],i+1,211,fHistMCNegPi[i]);
 379     if(ITSsa->IsGoodBin(i,211)){
 380       ITSsa->GetFitPar(fpar,efpar);
 381       hMeanNeg[0]->SetBinContent(i+1,fpar[1]);
 382       hMeanNeg[0]->SetBinError(i+1,efpar[1]);
 383       hSigmaNeg[0]->SetBinContent(i+1,fpar[2]);
 384       hSigmaNeg[0]->SetBinError(i+1,efpar[2]);
 385     }
 386     cgausPineg->Update();
 387     //kaons
 388     cgausKneg->cd(i+1);
 389     gPad->SetLogy();
 390     fHistNegK[i]->Write();
 391     gres[4][i]=new TGraph();
 392     ITSsa->DoFit(fHistNegK[i],i,321,gres[4][i]);
 393     ITSsa->FillHisto(hSpectraNeg[1],i+1,binsize,321);
 394     if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDNeg[1],i+1,321,fHistMCNegK[i]);
 395     if(ITSsa->IsGoodBin(i,321)){
 396       ITSsa->GetFitPar(fpar,efpar);
 397       hMeanNeg[1]->SetBinContent(i+1,fpar[1]);
 398       hMeanNeg[1]->SetBinError(i+1,efpar[1]);
 399       hSigmaNeg[1]->SetBinContent(i+1,fpar[2]);
 400       hSigmaNeg[1]->SetBinError(i+1,efpar[2]);
 401     }
 402     cgausKneg->Update();
 403     //protons
 404     cgausPneg->cd(i+1);
 405     gPad->SetLogy();
 406     fHistNegP[i]->Write();
 407     gres[5][i]=new TGraph();
 408     ITSsa->DoFitProton(fHistNegP[i],i,2212,gres[5][i]);
 409     ITSsa->FillHisto(hSpectraNeg[2],i+1,binsize,2212);
 410     if(optMC) ITSsa->FillHistoMC(hSpectraMCPIDNeg[2],i+1,2212,fHistMCNegP[i]);
 411     if(ITSsa->IsGoodBin(i,2212)){
 412       ITSsa->GetFitPar(fpar,efpar);
 413       hMeanNeg[2]->SetBinContent(i+1,fpar[1]);
 414       hMeanNeg[2]->SetBinError(i+1,efpar[1]);
 415       hSigmaNeg[2]->SetBinContent(i+1,fpar[2]);
 416       hSigmaNeg[2]->SetBinError(i+1,efpar[2]);
 417     }
 418     cgausPneg->Update();
 419   }
 420
 421   //save histograms in the ps file
 422   cgausPipos->Print(ps1.Data());
 423   cgausKpos->Print(ps1.Data());
 424   cgausPpos->Print(ps1.Data());
 425   cgausPineg->Print(ps1.Data());
 426   cgausKneg->Print(ps1.Data());
 427   cgausPneg->Print(ps1.Data());
 428
 429   //spectra REC
 430   TCanvas *cspecREC=new TCanvas("cspecREC","SPECTRA rec",1000,800);
 431   cspecREC->Divide(2,1);
 432   cspecREC->cd(1);
 433   gPad->SetLogy();
 434   gPad->SetGridy();
 435   for(Int_t i=0;i<3;i++){
 436     if(i==0){
 437       hSpectraPos[i]->Draw("text");
 438       hSpectraPos[i]->SetTitle("ITSsa RAW positive");
 439       hSpectraPos[i]->SetMinimum(10);
 440       hSpectraPos[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
 441       hSpectraPos[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 442     }
 443     hSpectraPos[i]->Draw("esame");
 444     hSpectraPos[i]->SetLineColor(i+2);
 445     hSpectraPos[i]->SetMarkerColor(i+2);
 446     hSpectraPos[i]->SetMarkerStyle(21+i);
 447     hSpectraPos[i]->Write();
 448   }
 449   TLegend *leg=new TLegend(0.51,0.11,0.84,0.35,NULL,"brNDC");
 450   leg->SetFillColor(0);
 451   leg->SetBorderSize(0);
 452   TLegendEntry *entry0=leg->AddEntry(hSpectraPos[0],"pions","p");
 453   entry0->SetTextColor(2);
 454   TLegendEntry *entry2=leg->AddEntry(hSpectraPos[1],"kaons","p");
 455   entry2->SetTextColor(3);
 456   TLegendEntry *entry4=leg->AddEntry(hSpectraPos[2],"protons","p");
 457   entry4->SetTextColor(4);
 458   leg->Draw("same");
 459
 460   cspecREC->cd(2);
 461   gPad->SetLogy();
 462   gPad->SetGridy();
 463   for(Int_t i=0;i<3;i++){
 464     if(i==0){
 465       hSpectraNeg[i]->Draw("e");
 466       hSpectraNeg[i]->SetTitle("ITSsa RAW negative");
 467       hSpectraNeg[i]->SetMinimum(10);
 468       hSpectraNeg[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
 469       hSpectraNeg[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 470     }
 471     hSpectraNeg[i]->Draw("esame");
 472     hSpectraNeg[i]->SetLineColor(i+2);
 473     hSpectraNeg[i]->SetMarkerColor(i+2);
 474     hSpectraNeg[i]->SetMarkerStyle(21+i);
 475     hSpectraNeg[i]->Write();
 476   }
 477
 478
 479   //FitParameters
 480   TCanvas *cFitPar=new TCanvas("cFitPar","FitParameters",1000,800);
 481   cFitPar->Divide(2,3);
 482   TLatex** tplus=new TLatex*[3];
 483   TLatex** tminus=new TLatex*[3];
 484   tplus[0]=new TLatex(0.7,0.75,"#pi^{+}");
 485   tplus[1]=new TLatex(0.7,0.75,"K^{+}");
 486   tplus[2]=new TLatex(0.7,0.75,"p");
 487   tminus[0]=new TLatex(0.7,0.65,"#pi^{-}");
 488   tminus[1]=new TLatex(0.7,0.65,"K^{-}");
 489   tminus[2]=new TLatex(0.7,0.65,"#bar{p}");
 490
 491   for(Int_t ipart=0; ipart<3; ipart++){
 492     tplus[ipart]->SetNDC();
 493     tplus[ipart]->SetTextColor(1);
 494     tplus[ipart]->SetTextSize(0.05);
 495     tminus[ipart]->SetNDC();
 496     tminus[ipart]->SetTextColor(4);
 497     tminus[ipart]->SetTextSize(0.05);
 498     cFitPar->cd(1+2*ipart);
 499     gPad->SetGridy();
 500     hMeanPos[ipart]->SetMarkerStyle(20);
 501     hMeanPos[ipart]->Draw("e");
 502     hMeanPos[ipart]->GetYaxis()->SetRangeUser(-1,1);
 503     hMeanNeg[ipart]->SetMarkerStyle(24);
 504     hMeanNeg[ipart]->SetMarkerColor(4);
 505     hMeanNeg[ipart]->SetLineColor(4);
 506     hMeanNeg[ipart]->Draw("esame");
 507     tplus[ipart]->Draw();
 508     tminus[ipart]->Draw();
 509     cFitPar->cd(2+2*ipart);
 510     gPad->SetGridy();
 511     hSigmaPos[ipart]->SetMarkerStyle(20);
 512     hSigmaPos[ipart]->Draw("e");
 513     hSigmaPos[ipart]->GetYaxis()->SetRangeUser(0.1,0.3);
 514     hSigmaNeg[ipart]->SetMarkerStyle(24);
 515     hSigmaNeg[ipart]->SetMarkerColor(4);
 516     hSigmaNeg[ipart]->SetLineColor(4);
 517     hSigmaNeg[ipart]->Draw("esame");
 518     hSigmaPos[ipart]->Write();
 519     hSigmaNeg[ipart]->Write();
 520     tplus[ipart]->Draw();
 521     tminus[ipart]->Draw();
 522   }
 523
 524   //plus/minus ratio plots
 525   TH1F* hRatioPions=new TH1F(*hSpectraPos[0]);
 526   hRatioPions->Divide(hSpectraNeg[0]);
 527   TH1F* hRatioKaons=new TH1F(*hSpectraPos[1]);
 528   hRatioKaons->Divide(hSpectraNeg[1]);
 529   TH1F* hRatioProtons=new TH1F(*hSpectraPos[2]);
 530   hRatioProtons->Divide(hSpectraNeg[2]);
 531
 532   TCanvas *cratios=new TCanvas("cratios","Ratios +/-",1000,800);
 533   cratios->Divide(1,3);
 534   cratios->cd(1);
 535   hRatioPions->SetMinimum(0.7);
 536   hRatioPions->SetMaximum(1.3);
 537   hRatioPions->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 538   hRatioPions->GetYaxis()->SetTitle("#pi^{+}/#pi^{-}");
 539   hRatioPions->Draw("e");
 540   cratios->cd(2);
 541   hRatioKaons->SetMinimum(0.7);
 542   hRatioKaons->SetMaximum(1.3);
 543   hRatioKaons->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 544   hRatioKaons->GetYaxis()->SetTitle("K^{+}/K^{-}");
 545   hRatioKaons->Draw("e");
 546   cratios->cd(3);
 547   hRatioProtons->SetMinimum(0.7);
 548   hRatioProtons->SetMaximum(1.3);
 549   hRatioProtons->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 550   hRatioProtons->GetYaxis()->SetTitle("p/#bar{p}");
 551   hRatioProtons->Draw("e");
 552
 553   //efficiency and correction factor histograms
 554   if(optMC){
 555     for(Int_t i=0;i<3;i++){
 556       hEffPos[i] = (TH1F*)hSpectraPos[i]->Clone(Form("hEffPos%d",i));
 557       hEffPos[i]->SetTitle(Form("hEffPos%d",i));
 558       hEffPos[i]->Divide(hEffPos[i], hSpectraPrimPosMC[i], 1.0, 1.0, "B");//binomial errors
 559       hEffPos[i]->SetLineColor(i+2);
 560       hEffPos[i]->SetMarkerColor(i+2);
 561       hEffPos[i]->SetMarkerStyle(21+i);
 562       hEffPos[i]->Write();
 563
 564       if(hSpectraPrimPosMCBefEvSel[i]){
 565         hCorrFacPos[i] = (TH1F*)hSpectraPos[i]->Clone(Form("hCorrFacPos%d",i));
 566         hCorrFacPos[i]->SetTitle(Form("hCorrFacPos%d",i));
 567         hCorrFacPos[i]->Divide(hCorrFacPos[i], hSpectraPrimPosMCBefEvSel[i], 1.0, 1.0, "B");//binomial errors
 568         hCorrFacPos[i]->SetLineColor(i+2);
 569         hCorrFacPos[i]->SetMarkerColor(i+2);
 570         hCorrFacPos[i]->SetMarkerStyle(25+i);
 571         hCorrFacPos[i]->Write();
 572       }
 573       hEffMCPIDPos[i] = (TH1F*)hSpectraMCPIDPos[i]->Clone(Form("hEffMCPIDPos%d",i));
 574       hEffMCPIDPos[i]->SetTitle(Form("hEffMCPIDPos%d",i));
 575       hEffMCPIDPos[i]->Divide(hEffMCPIDPos[i], hSpectraPrimPosMC[i], 1.0, 1.0, "B");//binomial errors
 576       hEffMCPIDPos[i]->SetLineColor(i+2);
 577       hEffMCPIDPos[i]->SetMarkerColor(i+2);
 578       hEffMCPIDPos[i]->SetMarkerStyle(25+i);
 579       hEffMCPIDPos[i]->Write();
 580
 581       if(hSpectraPrimPosMCBefEvSel[i]){
 582         hCorrFacMCPIDPos[i] = (TH1F*)hSpectraMCPIDPos[i]->Clone(Form("hCorrFacMCPIDPos%d",i));
 583         hCorrFacMCPIDPos[i]->SetTitle(Form("hCorrFacMCPIDPos%d",i));
 584         hCorrFacMCPIDPos[i]->Divide(hCorrFacMCPIDPos[i], hSpectraPrimPosMCBefEvSel[i], 1.0, 1.0, "B");//binomial errors
 585         hCorrFacMCPIDPos[i]->SetLineColor(i+2);
 586         hCorrFacMCPIDPos[i]->SetMarkerColor(i+2);
 587         hCorrFacMCPIDPos[i]->SetMarkerStyle(25+i);
 588         hCorrFacMCPIDPos[i]->Write();
 589       }
 590
 591       hEffNeg[i] = (TH1F*)hSpectraNeg[i]->Clone(Form("hEffNeg%d",i));
 592       hEffNeg[i]->SetTitle(Form("hEffNeg%d",i));
 593       hEffNeg[i]->Divide(hEffNeg[i], hSpectraPrimNegMC[i], 1.0, 1.0, "B");//binomial errors
 594       hEffNeg[i]->SetLineColor(i+2);
 595       hEffNeg[i]->SetMarkerColor(i+2);
 596       hEffNeg[i]->SetMarkerStyle(21+i);
 597       hEffNeg[i]->Write();
 598
 599       if(hSpectraPrimNegMCBefEvSel[i]){
 600         hCorrFacNeg[i] = (TH1F*)hSpectraNeg[i]->Clone(Form("hCorrFacNeg%d",i));
 601         hCorrFacNeg[i]->SetTitle(Form("hCorrFacNeg%d",i));
 602         hCorrFacNeg[i]->Divide(hCorrFacNeg[i], hSpectraPrimNegMCBefEvSel[i], 1.0, 1.0, "B");//binomial errors
 603         hCorrFacNeg[i]->SetLineColor(i+2);
 604         hCorrFacNeg[i]->SetMarkerColor(i+2);
 605         hCorrFacNeg[i]->SetMarkerStyle(25+i);
 606         hCorrFacNeg[i]->Write();
 607       }
 608
 609       hEffMCPIDNeg[i] = (TH1F*)hSpectraMCPIDNeg[i]->Clone(Form("hEffMCPIDNeg%d",i));
 610       hEffMCPIDNeg[i]->SetTitle(Form("hEffMCPIDNeg%d",i));
 611       hEffMCPIDNeg[i]->Divide(hEffMCPIDNeg[i], hSpectraPrimNegMC[i], 1.0, 1.0, "B");//binomial errors
 612       hEffMCPIDNeg[i]->SetLineColor(i+2);
 613       hEffMCPIDNeg[i]->SetMarkerColor(i+2);
 614       hEffMCPIDNeg[i]->SetMarkerStyle(25+i);
 615       hEffMCPIDNeg[i]->Write();
 616
 617       if(hSpectraPrimPosMCBefEvSel[i]){
 618         hCorrFacMCPIDNeg[i] = (TH1F*)hSpectraMCPIDPos[i]->Clone(Form("hCorrFacMCPIDNeg%d",i));
 619         hCorrFacMCPIDNeg[i]->SetTitle(Form("hCorrFacMCPIDNeg%d",i));
 620         hCorrFacMCPIDNeg[i]->Divide(hCorrFacMCPIDNeg[i], hSpectraPrimNegMCBefEvSel[i], 1.0, 1.0, "B");//binomial errors
 621         hCorrFacMCPIDNeg[i]->SetLineColor(i+2);
 622         hCorrFacMCPIDNeg[i]->SetMarkerColor(i+2);
 623         hCorrFacMCPIDNeg[i]->SetMarkerStyle(25+i);
 624         hCorrFacMCPIDNeg[i]->Write();
 625       }
 626
 627       hSpectraMCPIDPos[i]->SetLineColor(i+2);
 628       hSpectraMCPIDPos[i]->SetMarkerColor(i+2);
 629       hSpectraMCPIDPos[i]->SetMarkerStyle(25+i);
 630       hSpectraMCPIDPos[i]->Write();
 631
 632       hSpectraMCPIDNeg[i]->SetLineColor(i+2);
 633       hSpectraMCPIDNeg[i]->SetMarkerColor(i+2);
 634       hSpectraMCPIDNeg[i]->SetMarkerStyle(25+i);
 635       hSpectraMCPIDNeg[i]->Write();
 636     }
 637
 638     //spectra REC Ideal PID
 639     TCanvas *cspecREC2=new TCanvas("cspecIdPid","SPECTRA rec Ideal PID",1000,800);
 640     cspecREC2->Divide(2,1);
 641     cspecREC2->cd(1);
 642     gPad->SetLogy();
 643     gPad->SetGridy();
 644     for(Int_t i=0;i<3;i++){
 645       if(i==0){
 646         hSpectraMCPIDPos[i]->Draw("text");
 647         hSpectraMCPIDPos[i]->SetTitle("ITSsa RAW positive - Ideal PID");
 648         hSpectraMCPIDPos[i]->SetMinimum(10);
 649         hSpectraMCPIDPos[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
 650         hSpectraMCPIDPos[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 651       }
 652       hSpectraMCPIDPos[i]->Draw("esame");
 653     }
 654     leg->Draw("same");
 655     cspecREC2->cd(2);
 656     gPad->SetLogy();
 657     gPad->SetGridy();
 658     for(Int_t i=0;i<3;i++){
 659       if(i==0){
 660         hSpectraMCPIDNeg[i]->Draw("e");
 661         hSpectraMCPIDNeg[i]->SetTitle("ITSsa RAW negative -Ideal PID");
 662         hSpectraMCPIDNeg[i]->SetMinimum(10);
 663         hSpectraMCPIDNeg[i]->GetYaxis()->SetTitle("d^{2}N/dp_{t}dy");
 664         hSpectraMCPIDNeg[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 665       }
 666       hSpectraMCPIDNeg[i]->Draw("esame");
 667     }
 668
 669     // Efficiency
 670     TCanvas *ceff=new TCanvas("ceff","ceff",1000,800);
 671     ceff->Divide(2,1);
 672     ceff->cd(1);
 673     gPad->SetGridy();
 674     for(Int_t i=0;i<3;i++){
 675       if(i==0){
 676         hEffPos[i]->SetTitle("ITSsa efficiency (positive)");
 677         hEffPos[i]->Draw("e");
 678         hEffPos[i]->GetYaxis()->SetRangeUser(0.1,0.9);
 679         hEffPos[i]->GetYaxis()->SetTitle("#epsilon");
 680         hEffPos[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 681       }
 682       hEffPos[i]->Draw("esame");
 683       hEffMCPIDPos[i]->Draw("esame");
 684     }
 685     TLegend *legeff=new TLegend(0.51,0.11,0.89,0.5,NULL,"brNDC");
 686     legeff->SetBorderSize(0);
 687     legeff->SetFillColor(0);
 688     TLegendEntry *entryeff0=legeff->AddEntry(hEffPos[0],"pions from fit","p");
 689     entryeff0->SetTextColor(2);
 690     TLegendEntry *entryeff2=legeff->AddEntry(hEffPos[1],"kaons from fit","p");
 691     entryeff2->SetTextColor(3);
 692     TLegendEntry *entryeff4=legeff->AddEntry(hEffPos[2],"protons from fit","p");
 693     entryeff4->SetTextColor(4);
 694     TLegendEntry *entryeff1=legeff->AddEntry(hEffMCPIDPos[0],"pions ideal PID","p");
 695     entryeff1->SetTextColor(2);
 696     TLegendEntry *entryeff3=legeff->AddEntry(hEffMCPIDPos[1],"kaons ideal PID","p");
 697     entryeff3->SetTextColor(3);
 698     TLegendEntry *entryeff5=legeff->AddEntry(hEffMCPIDPos[2],"protons ideal PID","p");
 699     entryeff5->SetTextColor(4);
 700     legeff->Draw("same");
 701
 702     ceff->cd(2);
 703     gPad->SetGridy();
 704     for(Int_t i=0;i<3;i++){
 705       if(i==0){
 706         hEffNeg[i]->SetTitle("ITSsa efficiency (negative)");
 707         hEffNeg[i]->Draw("e");
 708         hEffNeg[i]->GetYaxis()->SetRangeUser(0.1,0.9);
 709         hEffNeg[i]->GetYaxis()->SetTitle("#epsilon");
 710         hEffNeg[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 711       }
 712       hEffNeg[i]->Draw("esame");
 713       hEffMCPIDNeg[i]->Draw("esame");
 714     }
 715     ceff->Update();
 716
 717     // Correction factor
 718     TCanvas *ccf=new TCanvas("ccf","Corr Factor",1000,800);
 719     ccf->Divide(2,1);
 720     ccf->cd(1);
 721     gPad->SetGridy();
 722     for(Int_t i=0;i<3;i++){
 723       if(i==0){
 724         hEffPos[i]->SetTitle("ITSsa efficiency (positive)");
 725         hEffPos[i]->Draw("e");
 726         hEffPos[i]->GetYaxis()->SetRangeUser(0.1,0.9);
 727         hEffPos[i]->GetYaxis()->SetTitle("#epsilon");
 728         hEffPos[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 729       }
 730       hEffPos[i]->Draw("esame");
 731       if(hSpectraPrimPosMCBefEvSel[i]) hCorrFacPos[i]->Draw("esame");
 732     }
 733     TLegend *legcf=new TLegend(0.51,0.11,0.89,0.5,NULL,"brNDC");
 734     legcf->SetBorderSize(0);
 735     legcf->SetFillColor(0);
 736     TLegendEntry *entrycf0=legcf->AddEntry(hEffPos[0],"pions - eff.","p");
 737     entrycf0->SetTextColor(2);
 738     TLegendEntry *entrycf2=legcf->AddEntry(hEffPos[1],"kaons  - eff.","p");
 739     entrycf2->SetTextColor(3);
 740     TLegendEntry *entrycf4=legcf->AddEntry(hEffPos[2],"protons  - eff.","p");
 741     entrycf4->SetTextColor(4);
 742     TLegendEntry *entrycf1=legcf->AddEntry(hCorrFacPos[0],"pions - Corr. Fac.","p");
 743     entrycf1->SetTextColor(2);
 744     TLegendEntry *entrycf3=legcf->AddEntry(hCorrFacPos[1],"kaons - Corr. Fac.","p");
 745     entrycf3->SetTextColor(3);
 746     TLegendEntry *entrycf5=legcf->AddEntry(hCorrFacPos[2],"protons - Corr. Fac.","p");
 747     entrycf5->SetTextColor(4);
 748     legcf->Draw("same");
 749
 750     ccf->cd(2);
 751     gPad->SetGridy();
 752     for(Int_t i=0;i<3;i++){
 753       if(i==0){
 754         hEffNeg[i]->SetTitle("ITSsa efficiency (negative)");
 755         hEffNeg[i]->Draw("e");
 756         hEffNeg[i]->GetYaxis()->SetRangeUser(0.1,0.9);
 757         hEffNeg[i]->GetYaxis()->SetTitle("#epsilon");
 758         hEffNeg[i]->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 759       }
 760       hEffNeg[i]->Draw("esame");
 761       if(hSpectraPrimNegMCBefEvSel[i]) hCorrFacNeg[i]->Draw("esame");
 762     }
 763     ccf->Update();
 764
 765     cspecREC2->Print(ps1.Data());
 766     ceff->Print(ps1.Data());
 767     ccf->Print(ps1.Data());
 768   }
 769
 770   //save canvas in the ps file
 771   cratios->Print(ps1.Data());
 772   cspecREC->Print(ps1.Data());
 773   cFitPar->Print(ps1.Data());
 774   cs->Print(ps1.Data());
 775   cDEDX->Print(ps1.Data());
 776   cDEDXdouble->Print(ps1.Data());
 777   cdummy->Print(ps2.Data());
 778   //    fout->Close();
 779   return;
 780
 781 }//end of the main
 782
 783
 784 //______________________________________________________________________
 785 Double_t BetheBloch(Double_t *mom, Double_t *mass) {
 786   // BB PHobos parametrization fine tuned by G. Ortona
 787   // on data and MC in June 2010
 788   Double_t bg=mom[0]/mass[0];
 789   Double_t par[5];
 790   Bool_t MC=kFALSE;
 791   if(MC){
 792     par[0]=1.39126e+02;//tuned on LHC10d (sim pass4)
 793     par[1]=2.33589e+01;
 794     par[2]=6.05219e-02;
 795     par[3]=2.04336e-01;
 796     par[4]=-4.99854e-04;
 797   }
 798   else{
 799     par[0]=5.33458e+04;//tuned on data pass6
 800     par[1]=1.65303e+01;
 801     par[2]=2.60065e-03;
 802     par[3]=3.59533e-04;
 803     par[4]=7.51168e-05;
 804   }
 805   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 806   Double_t gamma=bg/beta;
 807   Double_t eff=1.0;
 808   if(bg<par[2]) eff=(bg-par[3])*(bg-par[3])+par[4];
 809   else eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 810   return (par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 811 }
 812
 813 //______________________________________________________________________
 814 void Logo(){
 815   //method to plot the ALICE logo in the propaganda plot
 816   TPaveText *tpave=new TPaveText(0.3,0.7,0.7,0.89,"brNDC");
 817   tpave->SetBorderSize(0);
 818   tpave->SetFillStyle(0);
 819   tpave->SetFillColor(0);
 820   TText *txt1=tpave->AddText("ALICE performance");
 821   TText *txt2=tpave->AddText("ITS stand-alone tracks");
 822   TText *txt3=tpave->AddText("pp @ #sqrt{s} = 7 TeV");
 823   txt1->SetTextFont(62);
 824   txt2->SetTextFont(62);
 825   txt3->SetTextFont(62);
 826   tpave->Draw();
 827   TPad *pad1=new TPad("pad1", "pad1", 0.75, 0.7, 0.89, 0.89);
 828   pad1->Draw();
 829   pad1->cd();
 830   TImage *img1 = TImage::Open("ALICElogo.gif");
 831   img1->FillRectangle(0);
 832   img1->Draw();
 833 }
 834
 835 //EOF
 836