PWG2/EVCHAR/correctSPDdNdEtapp.C

   1 /*************************************************************************
   2 * Macro correctSPDdNdEta                                                 *
   3 * To read data and correction histograms produced with                   *
   4 * AliAnalysisTaskSPDdNdEta and apply corrections to data                 *
   5 *                                                                        *
   6 * Author:  M. Nicassio (INFN Bari)                                       *
   7 * Contact: Maria.Nicassio@ba.infn.it, Domenico.Elia@ba.infn.it           *
   8 **************************************************************************/
   9
  10 #include "Riostream.h"
  11 #include "TROOT.h"
  12 #include "TStyle.h"
  13 #include "TMath.h"
  14 #include "TH1.h"
  15 #include "TH2.h"
  16 #include "TFile.h"
  17 #include "TList.h"
  18 #include "TTree.h"
  19 #include "TLegend.h"
  20 #include "TLegendEntry.h"
  21 #include "TCanvas.h"
  22
  23 void correctSPDdNdEtapp (Char_t* fileRaw, Char_t* fileCorr, Char_t* fileMC, Char_t* fileAccPb,
  24                        Bool_t kPrimariesTest, Bool_t kAccPb,
  25                        Int_t binVtxStart, Int_t binVtxStop,
  26                        Double_t cutAccPb, Double_t cutAcc, Double_t cutBkg,
  27                        Double_t cutAlgEff, Double_t cutDec, Double_t cutVtx, Double_t cutTrig)  {
  28
  29 //lim vtx: 5-15 ->[-10cm,10cm[; 0-20 ->[-20cm,20cm[
  30
  31   gStyle->SetOptLogy(kFALSE);
  32   gStyle->SetFrameLineWidth(2);
  33   gStyle->SetFrameFillColor(0);
  34   gStyle->SetCanvasColor(0);
  35   gStyle->SetPadColor(0);
  36   gStyle->SetHistLineWidth(2);
  37   gStyle->SetLabelSize(0.05, "x");
  38   gStyle->SetLabelSize(0.05, "y");
  39   gStyle->SetTitleSize(0.05, "x");
  40   gStyle->SetTitleSize(0.05, "y");
  41   gStyle->SetTitleOffset(1.1, "x");
  42   gStyle->SetTitleOffset(0.9, "y");
  43   gStyle->SetPadBottomMargin(0.14);
  44   gStyle->SetPalette(1,0);
  45   gStyle->SetTitleFillColor(0);
  46   gStyle->SetStatColor(0);
  47
  48   //Input files
  49   TFile *f_accPb = new TFile(fileAccPb);
  50   TFile *f_corr = new TFile(fileCorr);
  51   TFile *f_MC = new TFile(fileMC);
  52   TFile *f_raw = new TFile(fileRaw);
  53
  54   TList *l_corr = (TList*)f_corr->Get("cOutput");
  55   TList *l_MC = (TList*)f_MC->Get("cOutput");
  56   TList *l_raw = (TList*)f_raw->Get("cOutput");
  57
  58   //Corrections at track level
  59   //Background
  60   TH2F *hBackgroundCorr = new TH2F("backgroundCorr","Background correction",60,-3.,3.,20,-20.,20.);
  61   hBackgroundCorr->Sumw2();
  62   TH2F *hBackgroundCorrNum = (TH2F*)l_corr->At(32);
  63
  64   Int_t corrHistoName;
  65   if (kPrimariesTest) {
  66     corrHistoName = 32;
  67   } else {
  68     corrHistoName = 33;
  69   }
  70
  71   TH2F *hBackgroundCorrDen = (TH2F*)l_corr->At(corrHistoName);
  72   hBackgroundCorr->Divide(hBackgroundCorrNum,hBackgroundCorrDen,1.,1.);
  73
  74   for (Int_t ieta=0; ieta<60; ieta++) {
  75     for (Int_t jz=0; jz<20; jz++) {
  76       Float_t den = hBackgroundCorrNum->GetBinContent(ieta+1,jz+1);
  77       if (den!=0) {
  78         Float_t num = hBackgroundCorrDen->GetBinContent(ieta+1,jz+1)-den;
  79         hBackgroundCorr->SetBinError(ieta+1,jz+1,(TMath::Sqrt(num)/den)/TMath::Power(hBackgroundCorr->GetBinContent(ieta+1,jz+1),2));
  80       } else {
  81         hBackgroundCorr->SetBinError(ieta+1,jz+1,0);
  82       }
  83     }
  84   }
  85
  86   //Algorithm efficiency
  87   TH2F *hTrackletRecEff = new TH2F("efficiencyCorr","Tracklet algorithm + SPD efficiency",60,-3,3,20,-20,20);
  88   hTrackletRecEff->Sumw2();
  89   TH2F *hTrackletEffCorr_num = new TH2F(*((TH2F*)l_corr->At(34)));
  90   hTrackletEffCorr_num->Rebin2D(2,2,"");
  91   hTrackletRecEff->Divide(hBackgroundCorrNum,hTrackletEffCorr_num,1.,1.,"b");   // 1/Efficiency=correction_weight
  92
  93   //Correction for particles that do not reach the sensitive layer
  94   TH2F *hDecPartLay2Corr = new TH2F("decPartLay2Corr","Correction for particles not reaching the detector",60,-3.,3.,20,-20.,20.);
  95   hDecPartLay2Corr->Sumw2();
  96   TH2F *hPrimaries_evtTrVtx = new TH2F(*((TH2F*) l_corr->At(35)));
  97   TH2F *hDetectablePrimariesLay2_evtTrVtx = new TH2F(*((TH2F*) l_corr->At(36)));
  98   hDetectablePrimariesLay2_evtTrVtx->Rebin2D(2,2,"");
  99   hDecPartLay2Corr->Divide(hDetectablePrimariesLay2_evtTrVtx,hPrimaries_evtTrVtx,1.,1.,"b");
 100
 101   //Acceptance
 102   TH2F *hAccCorr = new TH2F("accCorr_test","SPD acceptance correction",120,-3,3,40,-20,20);
 103   hAccCorr->Sumw2();
 104   hAccCorr->Rebin2D();
 105   TH2F *hAccNum = new TH2F(*((TH2F*)l_corr->At(34)));
 106   hAccNum->Sumw2();
 107   TH2F *hAccDen = new TH2F(*((TH2F*)l_corr->At(36)));
 108   hAccDen->Sumw2();
 109   hAccNum->Rebin2D();
 110   hAccDen->Rebin2D();
 111   hAccCorr->Divide(hAccNum,hAccDen,1.,1.,"b");
 112   TH2F *hAccCorrRebin = new TH2F(*hAccCorr);
 113   hAccCorrRebin->Sumw2();
 114 //  hAccCorrRebin->Rebin2D();
 115
 116   TH2F *hAccCorrPb = (TH2F*) f_accPb->Get("AccTrac");
 117
 118   //Vertex-trigger correction
 119   TH2F *hVertexTrigCorrEtaNum = (TH2F*)l_corr->At(37);
 120   TH2F *hTriggerCorrEta_norm  = (TH2F*)l_corr->At(38);
 121   TH2F *hVertexCorrEta_norm   = (TH2F*)l_corr->At(35);
 122   TH2F *hTrigCorrEtaNumNSD = (TH2F*)l_corr->At(39);
 123
 124   TH2F *hVertexCorrEta = new TH2F("vertexCorrEta","",60,-3.,3.,20,-20.,20.);
 125   hVertexCorrEta->Sumw2();
 126   hVertexCorrEta->Divide(hVertexCorrEta_norm,hTriggerCorrEta_norm,1.,1.,"b");
 127
 128   TH2F *hTriggerCorrEta = new TH2F("triggerCorrEta","",60,-3.,3.,20,-20.,20.);
 129   hTriggerCorrEta->Sumw2();
 130   hTriggerCorrEta->Divide(hTriggerCorrEta_norm,hVertexTrigCorrEtaNum,1.,1.,"b");
 131
 132   TH2F *hVertexTrigCorrEta = new TH2F("vertexTrigEtaCorr","Correction for MB trigger and vertexer inefficiency",60,-3.,3.,20,-20.,20.);
 133   hVertexTrigCorrEta->Sumw2();
 134   hVertexTrigCorrEta->Divide(hVertexCorrEta_norm,hVertexTrigCorrEtaNum,1.,1.,"b");
 135
 136   TH2F *hVertexTrigCorrEtaNSD = new TH2F("vertexTrigEtaCorrNSD","",60,-3.,3.,20,-20.,20.);
 137   hVertexTrigCorrEtaNSD->Sumw2();
 138   hVertexTrigCorrEtaNSD->Divide(hTrigCorrEtaNumNSD,hVertexCorrEta_norm,1.,1.);
 139
 140   TH2F *hTrigCorrEtaNSD = new TH2F("triggerCorrEtaNSD","",60,-3.,3.,20,-20.,20.);
 141   hTrigCorrEtaNSD->Divide(hTriggerCorrEta_norm,hTrigCorrEtaNumNSD,1.,1.);
 142   TH2F *hTrigCorrEtaNSDInv = new TH2F("triggerCorrEtaNSDInv","",60,-3.,3.,20,-20.,20.);
 143   hTrigCorrEtaNSDInv->Divide(hTrigCorrEtaNumNSD,hTriggerCorrEta_norm,1.,1.);
 144
 145   TH2F *hTracksinTriggeredEvts_NSD = new TH2F(*(TH2F*)l_corr->At(40));
 146   for (Int_t ieta=0; ieta<60; ieta++) {
 147     for (Int_t jz=0; jz<20; jz++) {
 148       Float_t den = hTrigCorrEtaNumNSD->GetBinContent(ieta+1,jz+1);
 149       if (den!=0) {
 150         Float_t acc=hTracksinTriggeredEvts_NSD->GetBinContent(ieta+1,jz+1)/hTrigCorrEtaNumNSD->GetBinContent(ieta+1,jz+1);
 151         Float_t binerr = TMath::Sqrt(acc*(1-acc)/den);
 152         Float_t nonbinerr = TMath::Sqrt(hTriggerCorrEta_norm->GetBinContent(ieta+1,jz+1)-hTracksinTriggeredEvts_NSD->GetBinContent(ieta+1,jz+1))/den;
 153         hTrigCorrEtaNSD->SetBinError(ieta+1,jz+1,TMath::Sqrt(TMath::Power(binerr,2)+TMath::Power(nonbinerr,2)));
 154       } else {
 155         hTrigCorrEtaNSD->SetBinError(ieta+1,jz+1,0);
 156       }
 157     }
 158   }
 159
 160   //Corrections at event level
 161   TH2F *hVertexTrigCorrNum = (TH2F*)l_corr->At(41);
 162   TH2F *hTriggerCorr_norm  = (TH2F*)l_corr->At(43);
 163   TH2F *hVertexCorr_norm   = (TH2F*)l_corr->At(42);
 164   TH2F *hTriggerCorrNumNSD = (TH2F*)l_corr->At(44);
 165
 166   Int_t nBinMultCorr = 28;
 167   Double_t binLimMultCorr[nBinMultCorr+1];
 168   for (Int_t i = 0; i<nBinMultCorr+1; ++i) {
 169     if (i<21) binLimMultCorr[i] = (Double_t) i;
 170     else if (i<27) binLimMultCorr[i] = (Double_t) 20 +5*(i-20);
 171     else if (i==27) binLimMultCorr[i] = 100;
 172     else if (i==28) binLimMultCorr[i] = 200;
 173   }
 174
 175
 176   TH2F *hVertexCorr = new TH2F("vertexCorr","Vertex inefficiency correction",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 177   hVertexCorr->Sumw2();
 178   hVertexCorr->Divide(hVertexCorr_norm,hTriggerCorr_norm,1.,1.,"b");
 179
 180   TH2F *hTriggerCorr = new TH2F("triggerCorr","MB trigger correction - inelastic events",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 181   hTriggerCorr->Sumw2();
 182   hTriggerCorr->Divide(hTriggerCorr_norm,hVertexTrigCorrNum,1.,1.,"b");
 183
 184   TH2F *hTriggerCorrNSD = new TH2F("triggerCorrNSD","MB trigger correction - NSD events",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 185   hTriggerCorrNSD->Sumw2();
 186   hTriggerCorrNSD->Divide(hTriggerCorr_norm,hTriggerCorrNumNSD,1.,1.);
 187   TH2F *hTriggerCorrNSD2 = new TH2F("triggerCorrNSD2","",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 188   hTriggerCorrNSD2->Sumw2();
 189   hTriggerCorrNSD2->Divide(hTriggerCorrNumNSD,hTriggerCorr_norm,1.,1.);
 190
 191
 192   TH2F *hTriggeredEvts_NSD = new TH2F(*(TH2F*)l_corr->At(45));
 193   for (Int_t im=0; im<28; im++) {
 194     for (Int_t jz=0; jz<20; jz++) {
 195       Float_t den = hTriggerCorrNumNSD->GetBinContent(im+1,jz+1);
 196       if (den!=0) {
 197         Float_t acc=hTriggeredEvts_NSD->GetBinContent(im+1,jz+1)/hTriggerCorrNumNSD->GetBinContent(im+1,jz+1);
 198         Float_t binerr = TMath::Sqrt(acc*(1-acc)/den);
 199         Float_t nonbinerr = TMath::Sqrt(hTriggerCorr_norm->GetBinContent(im+1,jz+1)-hTriggeredEvts_NSD->GetBinContent(im+1,jz+1))/den;
 200         hTriggerCorrNSD->SetBinError(im+1,jz+1,TMath::Sqrt(TMath::Power(binerr,2)+TMath::Power(nonbinerr,2)));
 201       } else {
 202         hTriggerCorrNSD->SetBinError(im+1,jz+1,0);
 203       }
 204     }
 205   }
 206
 207   TH2F *hVertexMC2D = (TH2F*)l_MC->At(41);
 208   TH1D *hMCvertex = new TH1D(*(hVertexMC2D->ProjectionY("MCvertex",0,-1,"e")));  //MC vertex distrib all events
 209 //  hMCvertex->Draw();
 210
 211   //Raw distributions
 212
 213   //...track level
 214   TH2F *hSPDEta_2Draw_accBin = new TH2F(*((TH2F*)l_raw->At(2)));
 215   hSPDEta_2Draw_accBin->SetNameTitle("SPDEta_2Draw_accBin","");
 216   TH2F *hSPDEta_2Draw = new TH2F(*((TH2F*)l_raw->At(2)));
 217   hSPDEta_2Draw->SetNameTitle("SPDEta_2Draw","Reconstructed tracklets");
 218
 219   //...event level
 220   TH2F *hESDVertexvsNtracklets = new TH2F(*((TH2F*)l_raw->At(0)));
 221   hESDVertexvsNtracklets->SetNameTitle("ESDVertexvsNtracklets","Triggered events with vertex");
 222   TH2F *hESDVertexvsNtracklets_trigEvts = new TH2F(*((TH2F*)l_raw->At(1)));
 223   hESDVertexvsNtracklets_trigEvts->SetNameTitle("ESDVertexvsNtracklets_trigEvts","Triggered events");
 224
 225   TH2F *hESDVertexvsNtracklets_Corr = new TH2F("ESDVertexvsNtracklets_Corr","",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 226   TH2F *hESDVertexvsNtracklets_CorrNSD = new TH2F("ESDVertexvsNtracklets_CorrNSD","",nBinMultCorr,binLimMultCorr,20,-20.,20.);
 227
 228   TH1D *hSPDvertex = (TH1D*)l_raw->At(19);
 229   if (kPrimariesTest) {
 230     hSPDEta_2Draw = new TH2F(*((TH2F*)l_MC->At(32)));
 231     hSPDEta_2Draw->SetNameTitle("SPDEta_2Draw","");
 232     hESDVertexvsNtracklets = new TH2F(*((TH2F*)l_MC->At(42)));
 233     hESDVertexvsNtracklets_trigEvts = new TH2F(*((TH2F*)l_MC->At(43)));
 234
 235     hSPDvertex = hVertexCorr_norm->ProjectionY("MCvertex_RV",0,-1,"e");
 236   }
 237
 238   if (!kPrimariesTest) hSPDEta_2Draw->Rebin2D(2,2,"");
 239
 240   //Vertex correction
 241
 242   //Computing the percentage of triggered events with vertex and tracklet mult!=0 for each rec vertex bin
 243   //to distribute triggered events with tracklet mult = 0 in vertex bins.
 244   TH1F* fHistZdistribTrigVtxEvt = new TH1F("fHistZdistribTrigVtxEvt", "",20,-20.,20.);
 245   fHistZdistribTrigVtxEvt->Sumw2();
 246   for (Int_t ivtx=0; ivtx<20; ivtx++) {
 247     Float_t fracZ = (hESDVertexvsNtracklets->ProjectionY("_y",2,hESDVertexvsNtracklets->GetNbinsX()+1,"e"))->GetBinContent(ivtx+1)/
 248                     Float_t (hESDVertexvsNtracklets->ProjectionY("_y",2,hESDVertexvsNtracklets->GetNbinsX()+1,"e")->Integral(0,hESDVertexvsNtracklets->GetNbinsY()+1));
 249
 250     fHistZdistribTrigVtxEvt->SetBinContent(ivtx+1,fracZ);
 251     //cout<<" Fractions of rec events: Num "<<hESDVertexvsNtracklets->ProjectionY("_y",2,hESDVertexvsNtracklets->GetNbinsX()+1,"e")->GetBinContent(ivtx+1)<<
 252     //      " Den  "<<hESDVertexvsNtracklets->ProjectionY("_y",2,hESDVertexvsNtracklets->GetNbinsX()+1,"e")->Integral(0,hESDVertexvsNtracklets->GetNbinsY()+1)<<
 253     //      " Frac: "<<fracZ<<endl;
 254   }
 255
 256   //Computing the MC correction for the previous %
 257   TH1F* fHistEvtCorrFactor = new TH1F("fHistEvtCorrFactor", "",20,-20.,20.);
 258   fHistEvtCorrFactor->Sumw2();
 259   for (Int_t ivtx=0; ivtx<20; ivtx++) {
 260
 261     Float_t fracZnum = (hTriggerCorr_norm->GetBinContent(1,ivtx+1))/Float_t(hTriggerCorr_norm->ProjectionY("_y",1,1,"e")->Integral(0,hTriggerCorr_norm->GetNbinsY()+1));
 262
 263     //cout<<"Vertex bin "<<ivtx+1<<" Content per bin "<<hTriggerCorr_norm->GetBinContent(1,ivtx+1)<<" Integral "<<hTriggerCorr_norm->ProjectionY("_y",1,1,"e")->Integral(0,hTriggerCorr_norm->GetNbinsY()+1)<<" MC fraction of triggered events in 0-multiplicity bin "<<fracZnum<<endl;
 264
 265     Float_t fracZden = ((hVertexCorr_norm->ProjectionY("_y",2,hVertexCorr_norm->GetNbinsX()+1,"e"))->GetBinContent(ivtx+1))/Float_t(hVertexCorr_norm->ProjectionY("_y",2,hVertexCorr_norm->GetNbinsX()+1,"e")->Integral(0,hVertexCorr_norm->GetNbinsY()+1));
 266
 267     //cout<<" MC fraction of rec events (multSPD>0). Num "<<(hVertexCorr_norm->ProjectionY("_y",2,hVertexCorr_norm->GetNbinsX()+1,"e"))->GetBinContent(ivtx+1)<<" Den "<<hVertexCorr_norm->ProjectionY("_y",2,hVertexCorr_norm->GetNbinsX()+1,"e")->Integral(0,hVertexCorr_norm->GetNbinsY()+1)<<" fractions "<<fracZden<<endl;
 268
 269     Float_t ratio =0.;
 270     Float_t ratioerr =0.;
 271     if (fracZden!=0.) {
 272       ratio = fracZnum/fracZden;
 273     }
 274     if(ratio!=0.) {
 275       fHistEvtCorrFactor->SetBinContent(ivtx+1,ratio);
 276       fHistEvtCorrFactor->SetBinError(ivtx+1,ratioerr);
 277     }
 278   }
 279
 280   TH1F* hprod = new TH1F("prod", "",20,-20.,20.);
 281   hprod->Sumw2();
 282
 283   hprod->Multiply(fHistEvtCorrFactor,fHistZdistribTrigVtxEvt,1.,1.);
 284
 285   //MC distributions
 286
 287   Int_t histoNameRV;
 288   if (kPrimariesTest) {
 289     histoNameRV = 46; //MC pseudorapidity, MC vertex
 290   } else {
 291     histoNameRV = 47; //MC pseudorapidity, rec vertex
 292   }
 293
 294   TH2F *Eta_RV = (TH2F*)l_MC->At(histoNameRV);
 295   TH2F *Eta = (TH2F*)l_MC->At(48);
 296 //  Eta->Draw();
 297
 298   TH1D *hdNdEta_RV = new TH1D("dNdEta_RV","Pseudorapidity ",60,-3,3);
 299   hdNdEta_RV=Eta_RV->ProjectionX("dNdEta_RV",binVtxStart+1,binVtxStop,"e");
 300   Double_t nEvtsRec=hSPDvertex->Integral(binVtxStart+1,binVtxStop);
 301   cout<<"Number of reconstructed events: "<<nEvtsRec <<endl;
 302   hdNdEta_RV->Scale(10./nEvtsRec);
 303
 304   TH1D *hdNdEta_Trigg = new TH1D("dNdEta_Trigg","Pseudorapidity ",60,-3,3);
 305   TH2F *hEta_Trigg = new TH2F(*((TH2F*)l_MC->At(38)));
 306   hdNdEta_Trigg=hEta_Trigg->ProjectionX("dNdEta_Trigg",0,-1,"e");
 307   Double_t nEvtsTrigg=0;
 308   nEvtsTrigg =hTriggerCorr_norm->Integral();
 309   hdNdEta_Trigg->Scale(10./nEvtsTrigg);
 310   cout<<"Number of triggered events: "<<nEvtsTrigg <<endl;
 311   TH1D *hdNdEta = new TH1D("dNdEta","Pseudorapidity ",60,-3.,3.);
 312   TH1D *hdNdEta10 = new TH1D("dNdEta10","Pseudorapidity ",60,-3.,3.);
 313   hdNdEta->Sumw2();
 314   hdNdEta = Eta->ProjectionX("dNdEta",0,-1,"e");
 315   hdNdEta10 = Eta->ProjectionX("dNdEta10",binVtxStart+1,binVtxStop,"e");
 316   Double_t nEvtsTot=hMCvertex->Integral(0,hMCvertex->GetNbinsX()+1);
 317   Double_t nEvtsTot10= hMCvertex->Integral(binVtxStart+1,binVtxStop);
 318
 319   cout<<"Number of generated events in +-20cm: "<<nEvtsTot <<endl;
 320   cout<<"Number of generated events in +-10cm: "<<nEvtsTot10 <<endl;
 321
 322   hdNdEta->Scale(10./nEvtsTot);
 323   hdNdEta10->Scale(10./nEvtsTot10);
 324   hdNdEta->GetXaxis()->SetTitle("Pseudorapidity #eta");
 325   hdNdEta->GetYaxis()->SetTitle("dN_{ch}/d#eta");
 326 //  hdNdEta->Draw();
 327
 328   TH1F *hdNdEtaNSD = (TH1F*)l_MC->At(51);
 329   hdNdEtaNSD->Sumw2();
 330   Double_t nEvtsNSD=0;
 331   TH1F *hProcessType = (TH1F*)l_MC->At(52);
 332   nEvtsNSD =hProcessType->GetBinContent(2);
 333 //   cout<<"NSD"<<endl;
 334   cout<<"Number of generated NSD events in +-20cm: "<<nEvtsNSD <<endl;
 335   hdNdEtaNSD->Scale(10./nEvtsNSD);
 336
 337   for (Int_t ibin=0; ibin<60; ibin++) hdNdEtaNSD->SetBinError(ibin+1,0);
 338
 339   TH1F *hdNdEtaMCInel = (TH1F*)l_MC->At(49);
 340   TH1F *hdNdEtaMCNSD = (TH1F*)l_MC->At(69);
 341   Double_t ntracksInel = hdNdEtaMCInel->GetBinContent(4);
 342   Double_t ntracksNSD = hdNdEtaMCNSD->GetBinContent(4);
 343   cout<<"dNdEtaMCInel: "<<ntracksInel/nEvtsTot<<endl;
 344   cout<<"dNdEtaMCNSD: "<<ntracksNSD/nEvtsNSD<<endl;
 345
 346   //Corrected distributions
 347   TH2F *hSPDEta_2D_bkgCorr =  new TH2F("SPDEta_2D_bkgCorr", "",60,-3.,3.,20,-20.,20.);
 348   hSPDEta_2D_bkgCorr->Sumw2();
 349   TH2F *hSPDEta_2D_bkgEffCorr =  new TH2F("SPDEta_2D_bkgEffCorr", "",60,-3.,3.,20,-20.,20.);
 350   hSPDEta_2D_bkgEffCorr->Sumw2();
 351   TH2F *hSPDEta_2D_bkgEffAccCorr =  new TH2F("SPDEta_2D_bkgEffAccCorr", "",120,-3.,3.,40,-20.,20.);
 352   hSPDEta_2D_bkgEffAccCorr->Sumw2();
 353   TH2F *hSPDEta_2D_fullyCorr =  new TH2F("SPDEta_2D_fullyCorr", "",60,-3.,3.,20,-20.,20.);
 354   hSPDEta_2D_fullyCorr->Sumw2();
 355   TH2F *hSPDEta_2D_triggeredEvts =  new TH2F("SPDEta_2D_triggeredEvts","",60,-3.,3.,20,-20.,20.);
 356   hSPDEta_2D_triggeredEvts->Sumw2();
 357   TH2F *hSPDEta_2D_allEvts =  new TH2F("SPDEta_2D_allEvts","",60,-3.,3.,20,-20.,20.);
 358   hSPDEta_2D_allEvts->Sumw2();
 359   TH2F *hSPDEta_2D_NSDEvts =  new TH2F("SPDEta_2D_NSDEvts","",60,-3.,3.,20,-20.,20.);
 360   hSPDEta_2D_NSDEvts->Sumw2();
 361
 362
 363   TH1F *hnorm_bkgEffAccCorr =  new TH1F("norm_bkgEffAccCorr", "",60,-3.,3.);
 364   hnorm_bkgEffAccCorr->Sumw2();
 365   TH1F *hnorm_fullyCorr =  new TH1F("norm_fullyCorr", "",60,-3.,3.);
 366   hnorm_fullyCorr->Sumw2();
 367   TH1F *hnorm_triggered =  new TH1F("norm_triggered", "",60,-3.,3.);
 368   hnorm_triggered->Sumw2();
 369   TH1F *hnorm =  new TH1F("norm", "",60,-3.,3.);
 370   hnorm->Sumw2();
 371   TH1F *hnormNSD =  new TH1F("normNSD", "",60,-3.,3.);
 372   hnormNSD->Sumw2();
 373
 374   TH1F *hnorm_gen =  new TH1F("norm_gen", "",60,-3.,3.);
 375   hnorm_gen->Sumw2();
 376
 377   //dN/deta
 378   TH1F *hdNdEta_raw = new TH1F("dNdEta_raw","", 60,-3.,3.);
 379   hdNdEta_raw->Sumw2();
 380   hdNdEta_raw->Add(hSPDEta_2Draw->ProjectionX("SPDEta_2Draw_x",binVtxStart+1,binVtxStop,"e"));
 381   hdNdEta_raw->Scale(10./nEvtsRec);
 382
 383   TH1F *hdNdEta_bkgCorr = new TH1F("dNdEta_bkgCorr","", 60,-3.,3.);
 384   hdNdEta_bkgCorr->Sumw2();
 385   TH1F *hdNdEta_bkgEffCorr = new TH1F("dNdEta_bkgEffCorr","", 60,-3.,3.);
 386   hdNdEta_bkgEffCorr->Sumw2();
 387   TH1F *hdNdEta_bkgEffAccCorr = new TH1F("dNdEta_bkgEffAccCorr","", 60,-3.,3.);
 388   hdNdEta_bkgEffAccCorr->Sumw2();
 389   TH1F *hdNdEta_fullyCorr = new TH1F("dNdEta_fullyCorr","", 60,-3.,3.);
 390   hdNdEta_fullyCorr->Sumw2();
 391   TH1F *hdNdEta_triggeredEvts = new TH1F("dNdEta_triggeredEvts","", 60,-3.,3.);
 392   hdNdEta_triggeredEvts->Sumw2();
 393   TH1F *hdNdEta_allEvts = new TH1F("dNdEta_allEvts","", 60,-3.,3.);
 394   hdNdEta_allEvts->Sumw2();
 395   TH1F *hdNdEta_NSDEvts = new TH1F("dNdEta_NSDEvts","", 60,-3.,3.);
 396   hdNdEta_NSDEvts->Sumw2();
 397
 398   //Cutting either on weights or on relative errors of weights
 399   // Acceptance Pb-Pb
 400   for (Int_t jeta=0; jeta<120; jeta++) {
 401     for (Int_t kvtx=0; kvtx<40; kvtx++) {
 402       if (hAccCorrPb->GetBinContent(jeta+1,kvtx+1) < cutAccPb) {
 403         hAccCorrPb->SetBinContent(jeta+1,kvtx+1,0.);
 404         hAccCorrPb->SetBinError(jeta+1,kvtx+1,0.);
 405       }
 406     }
 407   }
 408   for (Int_t jeta=0; jeta<60; jeta++) {
 409     for (Int_t kvtx=0; kvtx<20; kvtx++) {
 410       //Acceptance p-p
 411       Float_t relErr = hAccCorr->GetBinError(jeta+1,kvtx+1)/hAccCorr->GetBinContent(jeta+1,kvtx+1);
 412       if (relErr>cutAcc) {
 413         hAccCorr->SetBinContent(jeta+1,kvtx+1,0.);
 414         hAccCorr->SetBinError(jeta+1,kvtx+1,0.);
 415       }
 416       // Background
 417       Float_t relBinErr = hBackgroundCorr->GetBinError(jeta+1,kvtx+1)/hBackgroundCorr->GetBinContent(jeta+1,kvtx+1);
 418       if (relBinErr>cutBkg) {
 419         hBackgroundCorr->SetBinContent(jeta+1,kvtx+1,0.);
 420         hBackgroundCorr->SetBinError(jeta+1,kvtx+1,0.);
 421       }
 422       //Alg eff
 423       relBinErr = hTrackletRecEff->GetBinError(jeta+1,kvtx+1)/hTrackletRecEff->GetBinContent(jeta+1,kvtx+1);
 424       if (relBinErr>cutAlgEff) {
 425         hTrackletRecEff->SetBinContent(jeta+1,kvtx+1,0.);
 426         hTrackletRecEff->SetBinError(jeta+1,kvtx+1,0.);
 427       }
 428       //Dec particles
 429       relBinErr = hDecPartLay2Corr->GetBinError(jeta+1,kvtx+1)/hDecPartLay2Corr->GetBinContent(jeta+1,kvtx+1);
 430       if (relBinErr>cutDec) {
 431         hDecPartLay2Corr->SetBinContent(jeta+1,kvtx+1,0.);
 432         hDecPartLay2Corr->SetBinError(jeta+1,kvtx+1,0.);
 433       }
 434       //Vtx Corr
 435       relBinErr = hVertexCorrEta->GetBinError(jeta+1,kvtx+1)/hVertexCorrEta->GetBinContent(jeta+1,kvtx+1);
 436       if (relBinErr>cutVtx) {
 437         hVertexCorrEta->SetBinContent(jeta+1,kvtx+1,0.);
 438         hVertexCorrEta->SetBinError(jeta+1,kvtx+1,0.);
 439       }
 440       relBinErr = hTriggerCorrEta->GetBinError(jeta+1,kvtx+1)/hTriggerCorrEta->GetBinContent(jeta+1,kvtx+1);
 441       if (relBinErr>cutTrig) {
 442         hTriggerCorrEta->SetBinContent(jeta+1,kvtx+1,0.);
 443         hTriggerCorrEta->SetBinError(jeta+1,kvtx+1,0.);
 444       }
 445       relBinErr = hTrigCorrEtaNSD->GetBinError(jeta+1,kvtx+1)/hTrigCorrEtaNSD->GetBinContent(jeta+1,kvtx+1);
 446       if (relBinErr>cutTrig) {
 447         hTrigCorrEtaNSD->SetBinContent(jeta+1,kvtx+1,0.);
 448         hTrigCorrEtaNSD->SetBinError(jeta+1,kvtx+1,0.);
 449       }
 450       if (hAccCorr->GetBinContent(jeta+1,kvtx+1)==0) {
 451         hDecPartLay2Corr->SetBinContent(jeta+1,kvtx+1,0.);  //show these correction only in bins inside the SPD acceptance
 452         hDecPartLay2Corr->SetBinError(jeta+1,kvtx+1,0.);
 453         hVertexTrigCorrEta->SetBinContent(jeta+1,kvtx+1,0.);
 454         hVertexTrigCorrEta->SetBinError(jeta+1,kvtx+1,0.);
 455         hVertexTrigCorrEtaNSD->SetBinContent(jeta+1,kvtx+1,0.);
 456         hVertexTrigCorrEtaNSD->SetBinError(jeta+1,kvtx+1,0.);
 457         hTrigCorrEtaNSDInv->SetBinContent(jeta+1,kvtx+1,0.);
 458         hTrigCorrEtaNSDInv->SetBinError(jeta+1,kvtx+1,0.);
 459
 460       }
 461     }
 462   }
 463
 464   //Apply corrections at
 465   //...track level
 466   hSPDEta_2D_bkgCorr->Add(hSPDEta_2Draw);
 467   hSPDEta_2D_bkgCorr->Multiply(hBackgroundCorr);
 468   hSPDEta_2D_bkgEffCorr->Divide(hSPDEta_2D_bkgCorr,hTrackletRecEff,1.,1.);
 469   if (kAccPb) {
 470     hSPDEta_2D_bkgEffAccCorr->Divide(hSPDEta_2Draw_accBin,hAccCorrPb,1.,1.);
 471     hSPDEta_2D_bkgEffAccCorr->Rebin2D(2,2,"");
 472   } else {
 473     hSPDEta_2D_bkgEffAccCorr->Rebin2D(2,2,"");
 474     hSPDEta_2D_bkgEffAccCorr->Divide(hSPDEta_2Draw,hAccCorr,1.,1.);
 475 //    hSPDEta_2D_bkgEffAccCorr->Divide(hSPDEta_2Draw_accBin,hAccCorr,1.,1.);
 476   }
 477 //  hSPDEta_2D_bkgEffAccCorr->Rebin2D(2,2,"");
 478   hSPDEta_2D_bkgEffAccCorr->Multiply(hBackgroundCorr);
 479   hSPDEta_2D_bkgEffAccCorr->Divide(hTrackletRecEff);
 480   hSPDEta_2D_fullyCorr->Add(hSPDEta_2D_bkgEffAccCorr);
 481   hSPDEta_2D_fullyCorr->Divide(hDecPartLay2Corr);
 482   hSPDEta_2D_triggeredEvts->Add(hSPDEta_2D_fullyCorr);
 483   hSPDEta_2D_triggeredEvts->Divide(hVertexCorrEta);
 484   hSPDEta_2D_allEvts->Add(hSPDEta_2D_triggeredEvts);
 485   hSPDEta_2D_allEvts->Divide(hTriggerCorrEta);
 486   hSPDEta_2D_NSDEvts->Add(hSPDEta_2D_triggeredEvts);
 487   hSPDEta_2D_NSDEvts->Divide(hTrigCorrEtaNSD);
 488
 489   //...event level
 490   Double_t evtTrigNoTracklets = hESDVertexvsNtracklets_trigEvts->ProjectionY("_y",1,1,"e")->Integral(0,hESDVertexvsNtracklets_trigEvts->GetNbinsY()+1);
 491   cout<<"Number of events triggered with tracklet multiplicity = 0: "<<evtTrigNoTracklets<<endl;
 492   for (Int_t ivtx=0;ivtx<20;ivtx++) {
 493     hESDVertexvsNtracklets_trigEvts->SetBinContent(1,ivtx+1,(hprod->GetBinContent(ivtx+1))*(evtTrigNoTracklets));
 494     hESDVertexvsNtracklets_trigEvts->SetBinError(1,ivtx+1,TMath::Sqrt((TMath::Power(hprod->GetBinContent(ivtx+1),2)*evtTrigNoTracklets)+(TMath::Power(evtTrigNoTracklets,2)*TMath::Power(hprod->GetBinError(ivtx+1),2))));
 495 //    cout<<"bincont "<<hESDVertexvsNtracklets_trigEvts->GetBinContent(1,ivtx+1)<<" binerr "<<hESDVertexvsNtracklets_trigEvts->GetBinError(1,ivtx+1)<<endl;
 496   }
 497
 498   hESDVertexvsNtracklets_Corr->Add(hESDVertexvsNtracklets_trigEvts);
 499   hESDVertexvsNtracklets_Corr->Divide(hTriggerCorr);
 500   hESDVertexvsNtracklets_CorrNSD->Add(hESDVertexvsNtracklets_trigEvts);
 501   hESDVertexvsNtracklets_CorrNSD->Divide(hTriggerCorrNSD);
 502
 503   //Filling normalization histograms
 504
 505   TH1D *hSPDvertexCorr_y = new TH1D("SPDvertexCorr_y","",20,-20,20);
 506   hSPDvertexCorr_y = hESDVertexvsNtracklets_Corr->ProjectionY("SPDvertexCorr_y",0,-1,"e");
 507   cout<<"Number of events (tracklets corrected)"<<hSPDvertexCorr_y->Integral(binVtxStart+1,binVtxStop)<<endl;
 508   for (Int_t ivtx=binVtxStart; ivtx<binVtxStop; ivtx++) {
 509     Double_t nEvtsivtx = hSPDvertex->GetBinContent(ivtx+1);
 510     Double_t nEvtsivtxCorr = hSPDvertexCorr_y->GetBinContent(ivtx+1);
 511     Double_t nEvtsivtxCorr_triggered = hESDVertexvsNtracklets_trigEvts->ProjectionY("_y",0,-1,"e")->GetBinContent(ivtx+1);
 512
 513     Double_t nEvtsivtxCorr_NSD = hESDVertexvsNtracklets_CorrNSD->ProjectionY("_y",0,-1,"e")->GetBinContent(ivtx+1);
 514
 515     //Double_t nEvtGen = hMCvertex->GetBinContent(ivtx+1);
 516     //cout<<"iBinVtx: "<<ivtx<<" nEvtsGen-> "<<nEvtGen<<" nEvtsCorr-> "<<nEvtsivtxCorr<<endl;
 517     for (Int_t jeta=0; jeta<60; jeta++) {
 518       if (hSPDEta_2D_bkgEffAccCorr->GetBinContent(jeta+1,ivtx+1)!=0.)
 519         hnorm_bkgEffAccCorr->Fill(hBackgroundCorr->GetXaxis()->GetBinCenter(jeta+1),nEvtsivtx);
 520       if (hSPDEta_2D_fullyCorr->GetBinContent(jeta+1,ivtx+1)!=0.)
 521         hnorm_fullyCorr->Fill(hBackgroundCorr->GetXaxis()->GetBinCenter(jeta+1),nEvtsivtx);
 522       if (hSPDEta_2D_triggeredEvts->GetBinContent(jeta+1,ivtx+1)!=0.)
 523         hnorm_triggered->Fill(hBackgroundCorr->GetXaxis()->GetBinCenter(jeta+1),nEvtsivtxCorr_triggered);
 524       if (hSPDEta_2D_allEvts->GetBinContent(jeta+1,ivtx+1)!=0.)
 525         hnorm->Fill(hBackgroundCorr->GetXaxis()->GetBinCenter(jeta+1),nEvtsivtxCorr);
 526       if (hSPDEta_2D_NSDEvts->GetBinContent(jeta+1,ivtx+1)!=0.)
 527         hnormNSD->Fill(hBackgroundCorr->GetXaxis()->GetBinCenter(jeta+1),nEvtsivtxCorr_NSD);
 528     }
 529   }
 530
 531   for (Int_t jeta=0; jeta<60; jeta++) {
 532     hnorm_bkgEffAccCorr->SetBinError(jeta+1,TMath::Sqrt(hnorm_bkgEffAccCorr->GetBinContent(jeta+1)));
 533     hnorm_fullyCorr->SetBinError(jeta+1,TMath::Sqrt(hnorm_fullyCorr->GetBinContent(jeta+1)));
 534     hnorm_triggered->SetBinError(jeta+1,TMath::Sqrt(hnorm_triggered->GetBinContent(jeta+1)));
 535     hnorm->SetBinError(jeta+1,TMath::Sqrt(hnorm->GetBinContent(jeta+1)));
 536     hnormNSD->SetBinError(jeta+1,TMath::Sqrt(hnorm->GetBinContent(jeta+1)));
 537   }
 538
 539   hdNdEta_bkgCorr->Add(hSPDEta_2D_bkgCorr->ProjectionX("SPDEta_2D_bkgCorr_x",binVtxStart+1,binVtxStop,"e"));
 540   hdNdEta_bkgCorr->Scale(1./nEvtsRec);
 541   hdNdEta_bkgEffCorr->Add(hSPDEta_2D_bkgEffCorr->ProjectionX("SPDEta_2D_bkgEffCorr_x",binVtxStart+1,binVtxStop,"e"));
 542   hdNdEta_bkgEffCorr->Scale(1./nEvtsRec);
 543   hdNdEta_bkgEffAccCorr->Divide(hSPDEta_2D_bkgEffAccCorr->ProjectionX("SPDEta_2D_bkgEffAccCorr_x",binVtxStart+1,binVtxStop,"e"),hnorm_bkgEffAccCorr);
 544   hdNdEta_fullyCorr->Divide(hSPDEta_2D_fullyCorr->ProjectionX("SPDEta_2D_fullyCorr_x",binVtxStart+1,binVtxStop,"e"),hnorm_fullyCorr);
 545   hdNdEta_triggeredEvts->Divide(hSPDEta_2D_triggeredEvts->ProjectionX("SPDEta_2D_triggeredEvts_x",binVtxStart+1,binVtxStop,"e"),hnorm_triggered);
 546   hdNdEta_allEvts->Divide(hSPDEta_2D_allEvts->ProjectionX("SPDEta_2D_allEvts_x",binVtxStart+1,binVtxStop,"e"),hnorm);
 547   hdNdEta_NSDEvts->Divide(hSPDEta_2D_NSDEvts->ProjectionX("SPDEta_2D_NSDEvts_x",binVtxStart+1,binVtxStop,"e"),hnormNSD);
 548
 549   hdNdEta_bkgCorr->Scale(10.);
 550   hdNdEta_bkgEffCorr->Scale(10.);
 551   hdNdEta_bkgEffAccCorr->Scale(10.);
 552   hdNdEta_fullyCorr->Scale(10.);
 553   hdNdEta_triggeredEvts->Scale(10.);
 554   hdNdEta_allEvts->Scale(10.);
 555   hdNdEta_NSDEvts->Scale(10.);
 556
 557   TH1F* hRationorm=new TH1F("rationorm","",60,-3,3);
 558   hRationorm->Add(hnorm);
 559   hRationorm->Scale(1./nEvtsTot10);
 560
 561   TH1D *hdNdEta_test = new TH1D("dNdEta_test","Pseudorapidity ",60,-3,3);
 562   TH2F *hMask_allEvts = new TH2F("Mask_allEvts","",60,-3,3,20,-20,20);
 563   hMask_allEvts->Sumw2();
 564
 565   TH2F *Eta_test = new TH2F("Eta_test","",60,-3.,3.,20,-20.,20.);
 566   if (kPrimariesTest) {
 567     for (Int_t kvtx=0; kvtx<20; kvtx++) {
 568       for (Int_t jeta=0; jeta<60; jeta++) {
 569         if (hSPDEta_2D_allEvts->GetBinContent(jeta+1,kvtx+1)!=0) {
 570           hMask_allEvts->SetBinContent(jeta+1,kvtx+1,1);
 571         }
 572       }
 573     }
 574     //dN/dEta gen test
 575     Eta_test->Multiply(Eta,hMask_allEvts,1.,1.);
 576     hdNdEta_test->Divide(Eta_test->ProjectionX("eta_x",binVtxStart+1,binVtxStop,"e"),hnorm,1.,1.);
 577     hdNdEta_test->Scale(10.);
 578   }
 579
 580   // Generated/corrected ratios
 581   TH2F* hRatiodNdEta2D=new TH2F("ratiodNdEta2D","",60,-3,3,20,-20,20);
 582 //  hRatiodNdEta2D->Divide(hSPDEta_2D_fullyCorr,Eta_RV,1.,1.);
 583   if (kPrimariesTest) hRatiodNdEta2D->Divide(hSPDEta_2D_allEvts,Eta_test,1.,1.);
 584   else hRatiodNdEta2D->Divide(hSPDEta_2D_allEvts,Eta,1.,1.);
 585
 586   TH1F* hRatiodNdEta=new TH1F("ratiodNdEta","",60,-3,3);
 587   hRatiodNdEta->GetXaxis()->SetTitle("Pseudorapidity #eta");
 588   hRatiodNdEta->GetYaxis()->SetTitle("Generated/Corrected");
 589
 590   TH1F* hRatiodNdEta10=new TH1F("ratiodNdEta10","",60,-3,3);
 591   hRatiodNdEta10->Divide(hdNdEta,hdNdEta10);
 592 //  new TCanvas();
 593 //  hRatiodNdEta10->Draw();
 594
 595   TH1F* hRatiodNdEtaNSD=new TH1F("ratiodNdEtaNSD","",60,-3,3);
 596   hRatiodNdEtaNSD->GetXaxis()->SetTitle("Pseudorapidity #eta");
 597   hRatiodNdEtaNSD->GetYaxis()->SetTitle("Generated/Corrected");
 598
 599   TH1F* hRatiodNdEta_vertex=new TH1F("ratiodNdEta_vertex","",60,-3,3);
 600   TH1F* hRatiodNdEta_triggered=new TH1F("ratiodNdEta_triggered","",60,-3,3);
 601
 602   for (Int_t i=0;i<60;i++) hdNdEta->SetBinError(i+1,0.);
 603   for (Int_t i=0;i<60;i++) hdNdEta_RV->SetBinError(i+1,0.);
 604
 605   if (kPrimariesTest) hRatiodNdEta->Divide(hdNdEta_test,hdNdEta_allEvts,1.,1.);
 606   else hRatiodNdEta->Divide(hdNdEta,hdNdEta_allEvts,1.,1.);
 607   hRatiodNdEtaNSD->Divide(hdNdEtaNSD,hdNdEta_NSDEvts,1.,1.);
 608   hRatiodNdEta_vertex->Divide(hdNdEta_RV,hdNdEta_fullyCorr,1.,1.);
 609   hRatiodNdEta_triggered->Divide(hdNdEta_Trigg,hdNdEta_triggeredEvts,1.,1.);
 610
 611   Float_t max = 14.;
 612
 613   //Draw a canvas with correction maps at track level
 614   TCanvas *ccorr = new TCanvas("c1","Tracklet level corrections and data");
 615   ccorr->Divide(3,2);
 616   ccorr->cd(1);
 617   hBackgroundCorr->GetXaxis()->SetTitle("Pseudorapidity #eta");
 618   hBackgroundCorr->GetYaxis()->SetTitle("z_{vertex} [cm]");
 619   hBackgroundCorr->GetXaxis()->SetRangeUser(-2,2);
 620   hBackgroundCorr->GetYaxis()->SetRangeUser(-10.,8.);
 621   hBackgroundCorr->SetMinimum(0.);
 622   hBackgroundCorr->SetMaximum(1.);
 623   hBackgroundCorr->Draw("colz");
 624   ccorr->cd(2);
 625   hTrackletRecEff->GetXaxis()->SetTitle("Pseudorapidity #eta");
 626   hTrackletRecEff->GetYaxis()->SetTitle("z_{vertex} [cm]");
 627   hTrackletRecEff->GetXaxis()->SetRangeUser(-2,2);
 628   hTrackletRecEff->GetYaxis()->SetRangeUser(-10.,8.);
 629   hTrackletRecEff->Draw("colz");
 630   ccorr->cd(3);
 631   hAccCorr->GetXaxis()->SetTitle("Pseudorapidity #eta");
 632   hAccCorr->GetYaxis()->SetTitle("z_{vertex} [cm]");
 633   hAccCorr->GetXaxis()->SetRangeUser(-2,2);
 634   hAccCorr->GetYaxis()->SetRangeUser(-10.,8.);
 635   hAccCorr->Draw("colz");
 636   ccorr->cd(4);
 637   hDecPartLay2Corr->GetXaxis()->SetTitle("Pseudorapidity #eta");
 638   hDecPartLay2Corr->GetYaxis()->SetTitle("z_{vertex} [cm]");
 639   hDecPartLay2Corr->GetXaxis()->SetRangeUser(-2,2);
 640   hDecPartLay2Corr->GetYaxis()->SetRangeUser(-10.,8.);
 641   hDecPartLay2Corr->SetMinimum(0.9);
 642   hDecPartLay2Corr->SetMaximum(1.);
 643   hDecPartLay2Corr->Draw("colz");
 644   ccorr->cd(5);
 645   hVertexTrigCorrEta->GetXaxis()->SetTitle("Pseudorapidity #eta");
 646   hVertexTrigCorrEta->GetYaxis()->SetTitle("z_{vertex} [cm]");
 647   hVertexTrigCorrEta->GetXaxis()->SetRangeUser(-2,2);
 648   hVertexTrigCorrEta->GetYaxis()->SetRangeUser(-10.,8.);
 649   hVertexTrigCorrEta->SetMinimum(0.99);
 650   hVertexTrigCorrEta->SetMaximum(1.01);
 651   hVertexTrigCorrEta->Draw("colz");
 652   ccorr->cd(6);
 653   hSPDEta_2Draw->GetXaxis()->SetTitle("Pseudorapidity #eta");
 654   hSPDEta_2Draw->GetYaxis()->SetTitle("z_{vertex} [cm]");
 655   hSPDEta_2Draw->GetXaxis()->SetRangeUser(-2,2);
 656   hSPDEta_2Draw->GetYaxis()->SetRangeUser(-10.,8.);
 657   hSPDEta_2Draw->Draw("colz");
 658
 659   //Draw a canvas with correction maps at event level
 660   TCanvas *ccorrev = new TCanvas("c2","Event level corrections and data");
 661   ccorrev->Divide(3,3);
 662   ccorrev->cd(1);
 663   hVertexCorr->GetXaxis()->SetTitle("Tracklet multiplicity");
 664   hVertexCorr->GetYaxis()->SetTitle("z_{MCvtx} [cm]");
 665   hVertexCorr->GetXaxis()->SetRangeUser(0.,20.);
 666   hVertexCorr->Draw("colz");
 667   ccorrev->cd(2);
 668   hTriggerCorr->GetXaxis()->SetTitle("Tracklet multiplicity");
 669   hTriggerCorr->GetYaxis()->SetTitle("z_{MCvtx} [cm]");
 670   hTriggerCorr->GetXaxis()->SetRangeUser(0.,20.);
 671   hTriggerCorr->Draw("colz");
 672   ccorrev->cd(3);
 673   hTriggerCorrNSD->GetXaxis()->SetTitle("Tracklet multiplicity");
 674   hTriggerCorrNSD->GetYaxis()->SetTitle("z_{MCvtx} [cm]");
 675   hTriggerCorrNSD->GetXaxis()->SetRangeUser(0.,20.);
 676   hTriggerCorrNSD->Draw("colz");
 677   ccorrev->cd(4);
 678   fHistZdistribTrigVtxEvt->GetYaxis()->SetTitle("Fraction of rec events per vertex bin");
 679   fHistZdistribTrigVtxEvt->GetXaxis()->SetTitle("z_{SPDvtx} [cm]");
 680   fHistZdistribTrigVtxEvt->Draw();
 681   ccorrev->cd(5);
 682   fHistEvtCorrFactor->GetYaxis()->SetTitle("Correction weight");
 683   fHistEvtCorrFactor->GetXaxis()->SetTitle("z_{MCvtx} [cm]");
 684   fHistEvtCorrFactor->Draw();
 685   ccorrev->cd(6);
 686   hESDVertexvsNtracklets->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 687   hESDVertexvsNtracklets->GetXaxis()->SetTitle("Tracklet mulitplicity");
 688   hESDVertexvsNtracklets->GetXaxis()->SetRangeUser(0.,20.);
 689   hESDVertexvsNtracklets->Draw("colz");
 690   ccorrev->cd(7);
 691   hESDVertexvsNtracklets_trigEvts->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 692   hESDVertexvsNtracklets_trigEvts->GetXaxis()->SetTitle("Tracklet mulitplicity");
 693   hESDVertexvsNtracklets_trigEvts->GetXaxis()->SetRangeUser(0.,20.);
 694   hESDVertexvsNtracklets_trigEvts->Draw("colz"); //bin mult=0 already distributed in z
 695
 696   // Draw dN/dEta
 697   hdNdEta_raw->SetLineColor(kGreen);
 698   hdNdEta_raw->SetLineWidth(3);
 699
 700   hdNdEta_RV->SetLineWidth(3);
 701   hdNdEta_RV->SetMinimum(0.);
 702   hdNdEta_RV->SetMaximum(max);
 703   hdNdEta_RV->SetLineColor(kRed);
 704 //  hdNdEta_RV->Draw("histo");
 705   hdNdEta->SetLineWidth(3);
 706   hdNdEta->SetMinimum(0.);
 707   hdNdEta->SetMaximum(max);
 708
 709   new TCanvas();
 710   hdNdEta->Draw("histo");
 711 //  hdNdEta_RV->Draw("histo,same");
 712
 713   hdNdEta_raw->Draw("histo,same");
 714
 715   hdNdEta_bkgCorr->SetMarkerStyle(21);
 716   hdNdEta_bkgCorr->SetMarkerColor(kBlue);
 717
 718   hdNdEta_bkgCorr->Draw("same,p,histo");
 719
 720   hdNdEta_bkgEffCorr->SetMarkerStyle(22);
 721   hdNdEta_bkgEffCorr->SetMarkerColor(kGreen);
 722
 723   hdNdEta_bkgEffAccCorr->SetMarkerStyle(23);
 724   hdNdEta_bkgEffAccCorr->SetMarkerColor(kRed);
 725
 726   hdNdEta_bkgEffCorr->Draw("p,same,e");
 727   hdNdEta_bkgEffAccCorr->Draw("p,same,e");
 728   hdNdEta_fullyCorr->SetMarkerStyle(20);
 729   hdNdEta_fullyCorr->SetMarkerColor(kBlue);
 730   hdNdEta_fullyCorr->Draw("p,same,e");
 731
 732   hdNdEta_triggeredEvts->SetMarkerStyle(21);
 733   hdNdEta_triggeredEvts->SetMarkerColor(kOrange);
 734   hdNdEta_triggeredEvts->Draw("p,same,e");
 735
 736   hdNdEta_allEvts->SetMarkerStyle(21);
 737   hdNdEta_allEvts->SetMarkerColor(kRed);
 738   hdNdEta_allEvts->Draw("p,same,e");
 739
 740   TLegend *leg1 = new TLegend(0.2,0.7,0.7,0.9,NULL,"brNDC");
 741   leg1->SetFillColor(0);
 742   leg1->SetTextFont(62);
 743   leg1->SetTextSize(0.0304569);
 744   TLegendEntry *entry1=leg1->AddEntry(hdNdEta_raw,"Reconstructed","l");
 745                 entry1=leg1->AddEntry(hdNdEta,"Generated","l");
 746                 entry1=leg1->AddEntry(hdNdEta_bkgCorr,"Bkg corrected","p");
 747                 entry1=leg1->AddEntry(hdNdEta_bkgEffCorr,"Bkg+AlgEff+SPDEff corrected","p");
 748                 entry1=leg1->AddEntry(hdNdEta_bkgEffAccCorr,"Bkg+AlgEff+SPDAccEff corrected","p");
 749                 entry1=leg1->AddEntry(hdNdEta_fullyCorr,"Corrected (triggered events with vertex)","p");
 750                 entry1=leg1->AddEntry(hdNdEta_triggeredEvts,"Corrected (triggered events)","p");
 751                 entry1=leg1->AddEntry(hdNdEta_allEvts,"Corrected (all events)","p");
 752
 753   leg1->Draw();
 754
 755   new TCanvas();
 756
 757   hdNdEta->Draw("histo");
 758   hdNdEta_NSDEvts->SetMarkerStyle(21);
 759   hdNdEtaNSD->SetLineColor(kGreen);
 760   hdNdEta_allEvts->Draw("p,same,e");
 761   hdNdEta_NSDEvts->Draw("same,p,e");
 762   hdNdEtaNSD->Draw("histo,same");
 763
 764   TLegend *leg2 = new TLegend(0.3,0.7,0.7,0.9,NULL,"brNDC");
 765   leg2->SetFillColor(0);
 766   leg2->SetTextFont(62);
 767
 768   TLegendEntry *entry2=leg2->AddEntry(hdNdEta_allEvts,"Corrected INEL","p");
 769                 entry2=leg2->AddEntry(hdNdEta_NSDEvts,"Corrected NSD","p");
 770                 entry2=leg2->AddEntry(hdNdEta,"Generated INEL","l");
 771                 entry2=leg2->AddEntry(hdNdEtaNSD,"Generated NSD","l");
 772
 773   leg2->Draw();
 774
 775   //Trigger-vertex efficiency
 776   TH1F* fHistMultAllNonDiff = (TH1F*)l_MC->At(63);
 777   TH1F* fHistMultAllSingleDiff = (TH1F*)l_MC->At(64);
 778   TH1F* fHistMultAllDoubleDiff = (TH1F*)l_MC->At(65);
 779   TH1F* fHistMultTrVtxNonDiff = (TH1F*)l_MC->At(66);
 780   TH1F* fHistMultTrVtxSingleDiff = (TH1F*)l_MC->At(67);
 781   TH1F* fHistMultTrVtxDoubleDiff = (TH1F*)l_MC->At(68);
 782
 783   TH1F* fHistTrigVtxEffvsMultEtacutNonDiff = new TH1F("fHistTrigVtxEffvsMultEtacutNonDiff","",20,-0.5,19.5);
 784   fHistTrigVtxEffvsMultEtacutNonDiff->GetXaxis()->SetTitle("Multiplicity in |#eta|<1.4");
 785   fHistTrigVtxEffvsMultEtacutNonDiff->GetYaxis()->SetTitle("Efficiency");
 786   fHistTrigVtxEffvsMultEtacutNonDiff->Sumw2();
 787   TH1F* fHistTrigVtxEffvsMultEtacutSingleDiff = new TH1F("fHistTrigVtxEffvsMultEtacutSingleDiff","",20,-0.5,19.5);
 788   fHistTrigVtxEffvsMultEtacutSingleDiff->Sumw2();
 789   TH1F* fHistTrigVtxEffvsMultEtacutDoubleDiff = new TH1F("fHistTrigVtxEffvsMultEtacutDoubleDiff","",20,-0.5,19.5);
 790   fHistTrigVtxEffvsMultEtacutDoubleDiff->Sumw2();
 791   fHistTrigVtxEffvsMultEtacutNonDiff->Divide(fHistMultTrVtxNonDiff,fHistMultAllNonDiff,1.,1.);
 792   fHistTrigVtxEffvsMultEtacutSingleDiff->Divide(fHistMultTrVtxSingleDiff,fHistMultAllSingleDiff,1.,1.);
 793   fHistTrigVtxEffvsMultEtacutDoubleDiff->Divide(fHistMultTrVtxDoubleDiff,fHistMultAllDoubleDiff,1.,1.);
 794
 795   new TCanvas();
 796
 797   fHistTrigVtxEffvsMultEtacutNonDiff->SetMarkerStyle(2);
 798   fHistTrigVtxEffvsMultEtacutSingleDiff->SetMarkerStyle(25);
 799   fHistTrigVtxEffvsMultEtacutDoubleDiff->SetMarkerStyle(24);
 800   fHistTrigVtxEffvsMultEtacutNonDiff->GetYaxis()->SetRangeUser(0.,1.1);
 801   fHistTrigVtxEffvsMultEtacutNonDiff->GetXaxis()->SetRangeUser(0.,10);
 802   fHistTrigVtxEffvsMultEtacutNonDiff->Draw("p,histo");
 803   fHistTrigVtxEffvsMultEtacutSingleDiff->Draw("same,p,histo");
 804   fHistTrigVtxEffvsMultEtacutDoubleDiff->Draw("same,p,histo");
 805
 806   TLegend *leg3 = new TLegend(0.3,0.2,0.7,0.4444,NULL,"brNDC");
 807   leg3->SetTextFont(62);
 808   leg3->SetFillColor(0);
 809   TLegendEntry *entry3=leg3->AddEntry(fHistTrigVtxEffvsMultEtacutNonDiff,"Non-diffractive","p");
 810                 entry3=leg3->AddEntry(fHistTrigVtxEffvsMultEtacutSingleDiff,"Single-diffractive","p");
 811                 entry3=leg3->AddEntry(fHistTrigVtxEffvsMultEtacutDoubleDiff,"Double-diffractive","p");
 812
 813   leg3->Draw();
 814
 815   TH1F* hProcessTypeTriggered = (TH1F*)l_MC->At(53);
 816   Float_t nInelTriggered = hProcessTypeTriggered->GetBinContent(3);
 817 //  Float_t nNSDTriggered = hProcessTypeTriggered->GetBinContent(2);
 818   Float_t nNDTriggered = hProcessTypeTriggered->GetBinContent(1);
 819   Float_t nSDTriggered = hProcessTypeTriggered->GetBinContent(4);
 820   Float_t nDDTriggered = hProcessTypeTriggered->GetBinContent(5);
 821   Float_t nND = hProcessType->GetBinContent(1);
 822   Float_t nSD = hProcessType->GetBinContent(4);
 823   Float_t nDD = hProcessType->GetBinContent(5);
 824
 825   cout<<"Number of events: Inel-->"<<hProcessType->GetBinContent(3)<<" ND-->"<<nND<<" SD-->"<<nSD<<" DD-->"<<nDD<<endl;
 826
 827   cout<<"Trigger efficiencies: Inel-->"<<nInelTriggered/nEvtsTot<<" ND-->"<<nNDTriggered/nND<<" SD-->"<<nSDTriggered/nSD<<" DD-->"<<nDDTriggered/nDD<<endl;
 828
 829   new TCanvas();
 830   hRatiodNdEtaNSD->GetYaxis()->SetRangeUser(0.9,1.1);
 831   hRatiodNdEtaNSD->GetXaxis()->SetRangeUser(-2.,1.9);
 832   hRatiodNdEtaNSD->Draw();
 833   new TCanvas();
 834   hRatiodNdEta->GetYaxis()->SetRangeUser(0.9,1.1);
 835   hRatiodNdEta->GetXaxis()->SetRangeUser(-2.,1.9);
 836   hRatiodNdEta->Draw();
 837
 838
 839   // Write histos
 840
 841   TFile *fout= new TFile("SPDdNdEta.root","RECREATE");
 842
 843   hAccCorr->Write();
 844   hAccCorrPb->Write();
 845
 846   hBackgroundCorr->Write();
 847   hTrackletRecEff->Write();
 848   hDecPartLay2Corr->Write();
 849
 850   hSPDEta_2Draw_accBin->Write();
 851   hSPDEta_2Draw->Write();
 852
 853   hSPDvertex->Write();
 854   hSPDvertexCorr_y->Write();
 855   hSPDEta_2D_bkgCorr->Write();
 856   hSPDEta_2D_bkgEffCorr->Write();
 857   hSPDEta_2D_bkgEffAccCorr->Write();
 858   hSPDEta_2D_fullyCorr->Write();
 859   hSPDEta_2D_triggeredEvts->Write();
 860   hSPDEta_2D_allEvts->Write();
 861   hSPDEta_2D_NSDEvts->Write();
 862
 863   hnorm_bkgEffAccCorr->Write();
 864   hnorm_triggered->Write();
 865   hnorm->Write();
 866   hnormNSD->Write();
 867   hnorm_gen->Write();
 868
 869   hVertexCorr->Write();
 870   hTriggerCorr->Write();
 871   hTriggerCorrNSD->Write();
 872   hTriggerCorrNSD2->Write();
 873
 874   hVertexCorrEta->Write();
 875   hTriggerCorrEta->Write();
 876   hVertexTrigCorrEta->Write();
 877   hVertexTrigCorrEtaNSD->Write();
 878   hTrigCorrEtaNSD->Write();
 879   hTrigCorrEtaNSDInv->Write();
 880
 881   fHistZdistribTrigVtxEvt->Write();
 882   fHistEvtCorrFactor->Write();
 883
 884   hESDVertexvsNtracklets->Write();
 885   hESDVertexvsNtracklets_trigEvts->Write();
 886   hESDVertexvsNtracklets_Corr->Write();
 887   hESDVertexvsNtracklets_CorrNSD->Write();
 888
 889   hdNdEta->Write();
 890   hdNdEta_RV->Write();
 891   hdNdEta_bkgCorr->Write();
 892   hdNdEta_bkgEffCorr->Write();
 893   hdNdEta_bkgEffAccCorr->Write();
 894   hdNdEta_fullyCorr->Write();
 895   hdNdEta_triggeredEvts->Write();
 896   hdNdEta_allEvts->Write();
 897   hdNdEta_NSDEvts->Write();
 898
 899   hdNdEta_raw->Write();
 900
 901   hRatiodNdEta->Write();
 902   hRatiodNdEtaNSD->Write();
 903   hRatiodNdEta_triggered->Write();
 904   hRatiodNdEta_vertex->Write();
 905   hRatiodNdEta2D->Write();
 906   hRationorm->Write();
 907
 908   Eta_test->Write();
 909   Eta->Write();
 910   Eta_RV->Write();
 911
 912   hprod->Write();
 913   hdNdEta10->Write();
 914   hRatiodNdEta10->Write();
 915
 916   fHistTrigVtxEffvsMultEtacutNonDiff->Write();
 917   fHistTrigVtxEffvsMultEtacutSingleDiff->Write();
 918   fHistTrigVtxEffvsMultEtacutDoubleDiff->Write();
 919
 920   fout->Write();
 921   fout->Close();
 922
 923 }