PWGPP/EVCHAR/correctSPDdNdEta.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 Float_t CalculateErrordNdEtaPerPartPair(Float_t a, Float_t b, Float_t da, Float_t db);
  24 void combscalingfactors (TList *lall, TList *lcomb, TList *lmc, Double_t *scaleNormRange_rot, Double_t *scaleNormRange_labels);
  25
  26 void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr,
  27                        Char_t* fileAlphaCorr, Char_t* fileMCBkgCorr, Char_t* fileMC,
  28                        Int_t binVtxStart=5, Int_t binVtxStop=15,
  29                        Double_t cutCorr=10, Double_t cutBkg=1., Double_t cutBkgMC=1.,
  30                        Bool_t bkgFromMCLabels = kFALSE,
  31                        Bool_t changeStrangeness = kFALSE)   {
  32
  33 // vtx lim 13-27 for +-7cm
  34 // vtx lim 15-25 for +-5cm
  35
  36
  37   Int_t nBinVtx = 40;
  38   Double_t MaxVtx = 20.;
  39
  40   const Int_t nBinMultCorr = 200;
  41   Float_t nMaxMult = 20000.;
  42   Double_t binLimMultCorr[nBinMultCorr+1];
  43   binLimMultCorr[0] = 0.;
  44   binLimMultCorr[1] = 1.;
  45   for (Int_t i = 2; i<=nBinMultCorr;++i) {
  46     binLimMultCorr[i] = (i-1)*nMaxMult/nBinMultCorr;
  47   }
  48
  49   const Int_t nBinEtaCorr = 60;
  50   Float_t etaMax = 3.;
  51   Float_t etaMin = -3.;
  52   Double_t *binLimEtaCorr = new Double_t[nBinEtaCorr+1];
  53   for (Int_t i = 0; i<nBinEtaCorr+1; ++i) {
  54     binLimEtaCorr[i] = (Double_t) etaMin+i*(etaMax*2.)/nBinEtaCorr;
  55   }
  56
  57   Float_t nBinsPerPseudorapidityUnit = nBinEtaCorr/(binLimEtaCorr[nBinEtaCorr]-binLimEtaCorr[0]);
  58   cout<<"Number of bins per pseudorapidity unit"<<nBinsPerPseudorapidityUnit<<endl;
  59
  60   gStyle->SetOptLogy(kFALSE);
  61   gStyle->SetFrameLineWidth(2);
  62   gStyle->SetFrameFillColor(0);
  63   gStyle->SetCanvasColor(0);
  64   gStyle->SetPadColor(0);
  65   gStyle->SetHistLineWidth(2);
  66   gStyle->SetLabelSize(0.05, "x");
  67   gStyle->SetLabelSize(0.05, "y");
  68   gStyle->SetTitleSize(0.05, "x");
  69   gStyle->SetTitleSize(0.05, "y");
  70   gStyle->SetTitleOffset(1.1, "x");
  71   gStyle->SetTitleOffset(0.9, "y");
  72   gStyle->SetPadBottomMargin(0.14);
  73   gStyle->SetPalette(1,0);
  74   gStyle->SetTitleFillColor(0);
  75   gStyle->SetStatColor(0);
  76
  77   //Input files
  78   TFile *f_MC = new TFile(fileMC);
  79   TFile *f_acorr = new TFile(fileAlphaCorr);
  80   TFile *f_mcbkgcorr = new TFile(fileMCBkgCorr);
  81   TFile *f_rawbkg = new TFile(fileRawBkgCorr);
  82   TFile *f_raw = new TFile(fileRaw);
  83
  84   TList *l_MC = (TList*)f_MC->Get("cOutput");
  85   TList *l_acorr = (TList*)f_acorr->Get("cOutput");
  86   TList *l_mcbkgcorr = (TList*)f_mcbkgcorr->Get("cOutput");
  87   TList *l_rawbkg = (TList*)f_rawbkg->Get("cOutput");
  88   TList *l_raw = (TList*)f_raw->Get("cOutput");
  89
  90   Double_t scaleBkg = 1.;
  91   Double_t scaleBkgMC = 1.;
  92   Double_t scaleBkgLabels = 1.;
  93
  94   cout<<" Background scaling factors for MC: "<<endl;
  95   combscalingfactors (l_acorr,l_mcbkgcorr,l_acorr,&scaleBkgMC,&scaleBkgLabels);
  96   cout<<" Background scaling factors for data: "<<endl;
  97   combscalingfactors (l_raw,l_rawbkg,l_acorr,&scaleBkg,&scaleBkgLabels);
  98
  99   //Histogram to be corrected at tracklet level
 100   TH2F *hSPDEta_2Draw = new TH2F(*((TH2F*)l_raw->FindObject("fHistSPDRAWEtavsZ")));
 101   hSPDEta_2Draw->SetNameTitle("SPDEta_2Draw","Reconstructed tracklets");
 102
 103   //Corrections at track level
 104   // Combinatorial background: beta correction from data
 105   TH2F *hCombBkg = (TH2F*)l_rawbkg->FindObject("fHistSPDRAWEtavsZ");
 106   hCombBkg->Scale(scaleBkg);
 107
 108   TH2F *hCombBkgCorrData = new TH2F("backgroundCorrDATA","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 109   hCombBkgCorrData->GetXaxis()->SetTitle("#eta");
 110   hCombBkgCorrData->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 111   hCombBkgCorrData->Sumw2();
 112
 113   // Combinatorial background: beta correction from MC
 114   TH2F *hCombBkgMC = (TH2F*)l_mcbkgcorr->FindObject("fHistSPDRAWEtavsZ");
 115   hCombBkgMC->Scale(scaleBkgMC);
 116   TH2F *hCombBkgMCDen = (TH2F*)l_acorr->FindObject("fHistSPDRAWEtavsZ");
 117
 118   TH2F *hCombBkgCorrMC = new TH2F("backgroundCorrMC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 119   hCombBkgCorrMC->GetXaxis()->SetTitle("#eta");
 120   hCombBkgCorrMC->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 121   hCombBkgCorrMC->Sumw2();
 122
 123   // Background correction from MC labels
 124   TH2F *hBkgCombLabels = (TH2F*)l_acorr->FindObject("fHistBkgCombLabels");
 125   if (bkgFromMCLabels) {
 126     hCombBkgCorrData->Divide(hBkgCombLabels,hSPDEta_2Draw,1.,1.);
 127     hCombBkgCorrData->Scale(scaleBkgLabels);
 128   } else {
 129     hCombBkgCorrData->Divide(hCombBkg,hSPDEta_2Draw,1.,1.);
 130   }
 131
 132   // 1-beta in DATA
 133   TH2F *hCombBkgCorr2Data = new TH2F("backgroundCorr2Data","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 134
 135   for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
 136     for (Int_t jz=0; jz<nBinVtx; jz++) {
 137       hCombBkgCorr2Data->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrData->GetBinContent(ieta+1,jz+1));
 138       hCombBkgCorr2Data->SetBinError(ieta+1,jz+1,hCombBkgCorrData->GetBinError(ieta+1,jz+1));
 139     }
 140   }
 141
 142   // Errors on beta
 143   //..data
 144   for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
 145     for (Int_t jz=0; jz<nBinVtx; jz++) {
 146       hCombBkgCorrData->SetBinError(ieta+1,jz+1,scaleBkg*TMath::Sqrt(hCombBkg->GetBinContent(ieta+1,jz+1))/hSPDEta_2Draw->GetBinContent(ieta+1,jz+1));
 147 //      hCombBkgCorrData->SetBinError(ieta+1,jz+1,hCombBkgCorrData->GetBincontent(ieta+1,jz+1)*TMath::Sqrt(hSPDEta_2Draw->GetBinContent(ieta+1,jz+1))/hSPDEta_2Draw->GetBinContent(ieta+1,jz+1));
 148     }
 149   }
 150
 151   // Alpha correction: MC only
 152   TH2F *hPrimaries_evtTrVtx = new TH2F(*((TH2F*) l_acorr->FindObject("fHistNonDetectableCorrNum"))); // all prim in sel events
 153   new TCanvas();
 154   hPrimaries_evtTrVtx->Draw();
 155
 156   TH2F *hInvAlphaCorr = new TH2F("InvAlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 157   hInvAlphaCorr->Sumw2();
 158
 159   TH2F *hGenProtons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistProtons")));
 160   TH2F *hGenKaons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistKaons")));
 161   TH2F *hGenPions = new TH2F(*((TH2F*) l_acorr->FindObject("fHistPions")));
 162   TH2F *hGenOthers = new TH2F(*((TH2F*) l_acorr->FindObject("fHistOthers")));
 163
 164   TH2F *hPrimReweighted = new TH2F("primReweighted","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 165   hPrimReweighted->Add(hGenKaons);
 166   hPrimReweighted->Scale(2.25); //double kaon fraction
 167   hPrimReweighted->Add(hGenProtons);
 168   hPrimReweighted->Add(hGenPions);
 169   hPrimReweighted->Add(hGenOthers); //use this in alpha if flag for syst true
 170
 171   TH2F *hRecProtons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedProtons")));
 172   TH2F *hRecKaons = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedKaons")));
 173   TH2F *hRecPions = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedPions")));
 174   TH2F *hRecOthers = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedOthers")));
 175   TH2F *hRecSec = new TH2F(*((TH2F*) l_acorr->FindObject("fHistReconstructedSec")));
 176
 177   TH2F *hRecReweighted = new TH2F("recReweighted","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 178   hRecReweighted->Add(hRecKaons);
 179   hRecReweighted->Scale(2.25);
 180   hRecReweighted->Add(hRecProtons);
 181   hRecReweighted->Add(hRecPions);
 182   hRecReweighted->Add(hRecOthers); //use this in alpha if flag for syst true
 183   hRecReweighted->Add(hRecSec);
 184   hRecReweighted->Add(hBkgCombLabels);
 185
 186   // 1-beta in MC
 187   TH2F *hCombBkgCorr2MC = new TH2F("backgroundCorr2MC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 188
 189   if (changeStrangeness) {
 190     if (bkgFromMCLabels) hCombBkgCorrMC->Divide(hBkgCombLabels,hRecReweighted,1.,1.);
 191     else hCombBkgCorrMC->Divide(hCombBkgMC,hRecReweighted,1.,1.);
 192
 193     for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++)
 194       for (Int_t jz=0; jz<nBinVtx; jz++)
 195         hCombBkgCorr2MC->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrMC->GetBinContent(ieta+1,jz+1));
 196
 197     hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hRecReweighted);
 198     hInvAlphaCorr->Divide(hPrimReweighted);
 199   } else {
 200
 201     if (bkgFromMCLabels) hCombBkgCorrMC->Divide(hBkgCombLabels,hCombBkgMCDen,1.,1.);
 202     else hCombBkgCorrMC->Divide(hCombBkgMC,hCombBkgMCDen,1.,1.);
 203
 204     for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++)
 205       for (Int_t jz=0; jz<nBinVtx; jz++)
 206         hCombBkgCorr2MC->SetBinContent(ieta+1,jz+1,1-hCombBkgCorrMC->GetBinContent(ieta+1,jz+1));
 207
 208     hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hCombBkgMCDen);
 209     hInvAlphaCorr->Divide(hPrimaries_evtTrVtx);
 210
 211   }
 212
 213 //  new TCanvas();
 214 //  hCombBkgCorr2MC->Draw();
 215
 216
 217   // Efficiency for particle species
 218   TH2F *hEffProtons = new TH2F("effProtons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 219   TH2F *hEffKaons   = new TH2F("effKaons","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 220   TH2F *hEffPions   = new TH2F("effPions","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 221   TH2F *hEffOthers  = new TH2F("effOthers","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 222
 223   hEffProtons->Divide(hRecProtons,hGenProtons,1.,1.);
 224   hEffKaons->Divide(hRecKaons,hGenKaons,1.,1.);
 225   hEffPions->Divide(hRecPions,hGenPions,1.,1.);
 226   hEffOthers->Divide(hRecOthers,hGenOthers,1.,1.);
 227
 228 //  hEffKaons->Divide(hEffPions);
 229
 230
 231   new TCanvas();
 232   hInvAlphaCorr->Draw();
 233
 234   // Errors on alpha
 235
 236   for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
 237     for (Int_t jz=0; jz<nBinVtx; jz++) {
 238       hInvAlphaCorr->SetBinError(ieta+1,jz+1,TMath::Sqrt((hCombBkgCorr2MC->GetBinContent(ieta+1,jz+1))*(hCombBkgMCDen->GetBinContent(ieta+1,jz+1)))/hPrimaries_evtTrVtx->GetBinContent(ieta+1,jz+1));
 239     }
 240   }
 241
 242   TH2F *hAlphaCorr = new TH2F("AlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 243   hAlphaCorr->GetXaxis()->SetTitle("#eta");
 244   hAlphaCorr->GetYaxis()->SetTitle("z_{vtx} [cm]");
 245
 246   for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
 247     for (Int_t jz=0; jz<nBinVtx; jz++) {
 248        hAlphaCorr->SetBinContent(ieta+1,jz+1,1./hInvAlphaCorr->GetBinContent(ieta+1,jz+1));
 249     }
 250   }
 251
 252   new TCanvas();
 253   hAlphaCorr->Draw();
 254
 255
 256   //////////////////////// Factorize alpha ///////////////////////////
 257
 258   ////////////////////////////////////////////////////////////////////
 259
 260   // Number of events and normalization histograms
 261   TH1D *hSPDvertex = (TH1D*)l_raw->FindObject("fHistSPDvtxAnalysis");
 262   Double_t nEvtsRec=hSPDvertex->Integral(binVtxStart+1,binVtxStop);
 263   cout<<"Number of reconstructed events in the selected vertex range: "<<nEvtsRec <<endl;
 264
 265   //Cutting corrections at the edges of the acceptance
 266   for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 267     for (Int_t kvtx=0; kvtx<nBinVtx; kvtx++) {
 268       if (hAlphaCorr->GetBinContent(jeta+1,kvtx+1)>cutCorr||hAlphaCorr->GetBinContent(jeta+1,kvtx+1)<0) {
 269         hAlphaCorr->SetBinContent(jeta+1,kvtx+1,0.);
 270         hAlphaCorr->SetBinError(jeta+1,kvtx+1,0.);
 271         hInvAlphaCorr->SetBinContent(jeta+1,kvtx+1,0.);
 272         hInvAlphaCorr->SetBinError(jeta+1,kvtx+1,0.);
 273       }
 274       if (hCombBkgCorrData->GetBinContent(jeta+1,kvtx+1)>cutBkg) {
 275         hCombBkgCorrData->SetBinContent(jeta+1,kvtx+1,0.);
 276         hCombBkgCorrData->SetBinError(jeta+1,kvtx+1,0.);
 277       }
 278       if (hCombBkgCorrMC->GetBinContent(jeta+1,kvtx+1)>cutBkgMC) {
 279         hCombBkgCorrMC->SetBinContent(jeta+1,kvtx+1,0.);
 280         hCombBkgCorrMC->SetBinError(jeta+1,kvtx+1,0.);
 281       }
 282     }
 283   }
 284
 285   new TCanvas();
 286   hCombBkgCorrMC->Draw();
 287   new TCanvas();
 288   hCombBkgCorrData->Draw();
 289
 290   // MC distribution
 291
 292   //... to compare corrected distribution to MC in selected events
 293   TH2F *hVertexMC2D = (TH2F*)l_MC->FindObject("fHistSelEvts");
 294   TH1D *hMCvertex = new TH1D(*(hVertexMC2D->ProjectionY("MCvertex",0,-1,"e")));  //MC vertex distrib
 295   new TCanvas();
 296   hMCvertex->Draw();
 297
 298   TH2F *Eta = (TH2F*)l_MC->FindObject("fHistNonDetectableCorrNum");
 299
 300   TH1D *hdNdEta = new TH1D("dNdEta","Pseudorapidity ",nBinEtaCorr,binLimEtaCorr);
 301   hdNdEta = Eta->ProjectionX("dNdEta",0,-1,"e"); //here already all events
 302   Double_t nEvtsTot=hMCvertex->Integral(0,hMCvertex->GetNbinsX()+1);
 303
 304   Double_t nEvtsTotSel= hMCvertex->Integral(binVtxStart+1,binVtxStop);
 305
 306   cout<<"Number of generated events: "<<nEvtsTot<<endl;
 307   cout<<"Number of generated events in the selected vertex range: "<<nEvtsTotSel <<endl;
 308
 309   Float_t nprimcentraleta = Eta->ProjectionX("dNdEta",0,-1,"e")->Integral(26,35); // project all or only in the sel vertex range?
 310   hdNdEta->Scale(nBinsPerPseudorapidityUnit/nEvtsTot);  // if so change number of events
 311   hdNdEta->GetXaxis()->SetTitle("#eta");
 312   hdNdEta->GetYaxis()->SetTitle("dN_{ch}/d#eta");
 313   for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta->SetBinError(i+1,0.);
 314   // dNdEta in +- 0.5
 315   cout<<"Monte Carlo dN/dEta in |eta|<0.5  (selected events for analysis and all events!): "<<nprimcentraleta/nEvtsTot<<endl;
 316   cout<<"Monte Carlo dN/dEta in |eta|<0.5  (selected events for analysis and events in vtx range!): "<<(Eta->ProjectionX("_x",binVtxStart+1,binVtxStop,"e")->Integral(26,35))/nEvtsTotSel<<endl;
 317   new TCanvas();
 318   hdNdEta->Draw();
 319
 320
 321   //Corrected distributions
 322   TH2F *hSPDEta_2D_bkgCorr =  new TH2F("SPDEta_2D_bkgCorr", "",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 323   hSPDEta_2D_bkgCorr->Sumw2();
 324
 325   TH2F *hSPDEta_2D_fullyCorr =  new TH2F("SPDEta_2D_fullyCorr", "",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 326   hSPDEta_2D_fullyCorr->Sumw2();
 327
 328   TH1F *hnorm_fullyCorr =  new TH1F("norm_fullyCorr", "",nBinEtaCorr,binLimEtaCorr);
 329
 330   //dN/deta
 331   TH1F *hdNdEta_raw = new TH1F("dNdEta_raw","",nBinEtaCorr,binLimEtaCorr);
 332   hdNdEta_raw->Sumw2();
 333   hdNdEta_raw->GetXaxis()->SetTitle("#eta");
 334   hdNdEta_raw->GetYaxis()->SetTitle("dN_{ch}/d#eta");
 335   hdNdEta_raw->Add(hSPDEta_2Draw->ProjectionX("SPDEta_2Draw_x",binVtxStart+1,binVtxStop,"e"));
 336   hdNdEta_raw->Scale(nBinsPerPseudorapidityUnit/nEvtsRec);
 337   new TCanvas();
 338   hdNdEta_raw->Draw();
 339
 340   TH1F *hdNdEta_bkgCorr = new TH1F("dNdEta_bkgCorr","",nBinEtaCorr,binLimEtaCorr);
 341   hdNdEta_bkgCorr->Sumw2();
 342   TH1F *hdNdEta_fullyCorr = new TH1F("dNdEta_fullyCorr","",nBinEtaCorr,binLimEtaCorr);
 343   hdNdEta_fullyCorr->Sumw2();
 344
 345   //Apply corrections at
 346   //...track level
 347   hSPDEta_2D_bkgCorr->Add(hSPDEta_2Draw);
 348   hSPDEta_2D_bkgCorr->Multiply(hCombBkgCorr2Data);
 349   hSPDEta_2D_fullyCorr->Add(hSPDEta_2D_bkgCorr);
 350   hSPDEta_2D_fullyCorr->Divide(hInvAlphaCorr);
 351   new TCanvas();
 352   hSPDEta_2D_bkgCorr->Draw();
 353   new TCanvas();
 354   hSPDEta_2D_fullyCorr->Draw();
 355
 356   //Filling normalization histogram
 357   for (Int_t ivtx=binVtxStart; ivtx<binVtxStop; ivtx++) {
 358     Double_t nEvtsivtx = hSPDvertex->GetBinContent(ivtx+1);
 359     for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 360       if (hAlphaCorr->GetBinContent(jeta+1,ivtx+1)!=0.) {
 361         hnorm_fullyCorr->Fill(hAlphaCorr->GetXaxis()->GetBinCenter(jeta+1),nEvtsivtx);
 362       }
 363     }
 364   }
 365
 366   for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 367     hnorm_fullyCorr->SetBinError(jeta+1,0);
 368   }
 369   new TCanvas();
 370   hnorm_fullyCorr->Draw();
 371
 372   hdNdEta_bkgCorr->Add(hSPDEta_2D_bkgCorr->ProjectionX("SPDEta_2D_bkgCorr_x",binVtxStart+1,binVtxStop,"e"));
 373   hdNdEta_bkgCorr->Scale(1./nEvtsRec);
 374   hdNdEta_fullyCorr->Divide(hSPDEta_2D_fullyCorr->ProjectionX("SPDEta_2D_fullyCorr_x",binVtxStart+1,binVtxStop,"e"),hnorm_fullyCorr);
 375
 376   TH1F *hCorrecteddNdEtaCentr = new TH1F("correcteddNdEtaCentr","",10,-0.5,0.5);
 377   hCorrecteddNdEtaCentr->Sumw2();
 378   for (Int_t jeta=0; jeta<10; jeta++) {
 379     hCorrecteddNdEtaCentr->SetBinContent(jeta+1,hdNdEta_fullyCorr->GetBinContent(26+jeta));
 380     hCorrecteddNdEtaCentr->SetBinError(jeta+1,hdNdEta_fullyCorr->GetBinError(26+jeta));
 381   }
 382   hCorrecteddNdEtaCentr->Rebin(10);
 383
 384   new TCanvas();
 385   hCorrecteddNdEtaCentr->Draw();
 386
 387   Float_t ncorrtrackletscentr = hSPDEta_2D_fullyCorr->ProjectionX("SPDEta_2D_fullyCorr_x",binVtxStart+1,binVtxStop,"e")->Integral(26,35);
 388   Float_t nevtscentr = hnorm_fullyCorr->Integral(26,35)/10.;
 389   cout<<""<<endl;
 390   cout<<"Corrected dN/dEta in |eta|<0.5: "<<ncorrtrackletscentr/nevtscentr<<endl;
 391   cout<<"Error on corrected tracklets in |eta|<0.5: +-"<<hCorrecteddNdEtaCentr->GetBinError(1)<<endl;
 392   // dN/dEta per participant pairs
 393
 394   Float_t NpartV0Glauber = 381.188;
 395   Float_t NpartCl2Glauber = 378.5;
 396   Float_t NpartCl2Hij = 365.3;
 397
 398   Float_t errNpartV0Glauber = 18.4951;
 399   Float_t errNpartCl2Glauber = 7.57; // 2% not used
 400   Float_t errNpartCl2Hij = 7.306;
 401
 402   Float_t dNdEtaCorr = ncorrtrackletscentr/nevtscentr;
 403   Float_t errdNdEtaCorr = hCorrecteddNdEtaCentr->GetBinError(1);
 404
 405   Float_t dNdEtaPartPairsV0G = dNdEtaCorr/(.5*NpartV0Glauber);
 406   Float_t dNdEtaPartPairsCl2G = dNdEtaCorr/(.5*NpartCl2Glauber);
 407   Float_t dNdEtaPartPairsCl2H = dNdEtaCorr/(.5*NpartCl2Hij);
 408
 409   Float_t errdNdEtaPartPairsV0G  = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartV0Glauber,errdNdEtaCorr,errNpartV0Glauber);
 410   Float_t errdNdEtaPartPairsCl2G = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartCl2Glauber,errdNdEtaCorr,errNpartCl2Glauber);
 411   Float_t errdNdEtaPartPairsCl2H = CalculateErrordNdEtaPerPartPair(dNdEtaCorr,NpartCl2Hij,errdNdEtaCorr,errNpartCl2Hij);
 412
 413   cout<<"V0 Glauber dN/dEta/.5<Npart> "<<dNdEtaPartPairsV0G<<" +- "<<errdNdEtaPartPairsV0G<<endl;
 414   cout<<"Cl2 Glauber dN/dEta/.5<Npart> "<<dNdEtaPartPairsCl2G<<" +- "<<errdNdEtaPartPairsCl2G<<endl;
 415   cout<<"Cl2 HIJING dN/dEta/.5<Npart> "<<dNdEtaPartPairsCl2H<<" +- "<<errdNdEtaPartPairsCl2H<<endl;
 416
 417   hdNdEta_bkgCorr->Scale(nBinsPerPseudorapidityUnit);
 418   hdNdEta_fullyCorr->Scale(nBinsPerPseudorapidityUnit);
 419
 420   // Create mask for MC prim in SPD acceptance
 421   TH1D *hdNdEta_test = new TH1D("dNdEta_test","Pseudorapidity ",nBinEtaCorr,binLimEtaCorr);
 422   TH2F *hMask = new TH2F("Mask","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 423   for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 424     for (Int_t kvtx=0; kvtx<nBinVtx; kvtx++) {
 425        if (hSPDEta_2D_fullyCorr->GetBinContent(jeta+1,kvtx+1)!=0) hMask->SetBinContent(jeta+1,kvtx+1,1.);
 426     }
 427   }
 428   hMask->Multiply(Eta);
 429   hdNdEta_test->Add(hMask->ProjectionX("Mask_x",binVtxStart+1,binVtxStop,"e"));
 430   hdNdEta_test->Divide(hnorm_fullyCorr);
 431 //  hdNdEta_test->Scale(1./nEvtsTotSel); // -> number of MC events
 432   hdNdEta_test->Scale(nBinsPerPseudorapidityUnit);
 433   for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta_test->SetBinError(i+1,0.);
 434
 435   TH1F* hRatiotest=new TH1F("ratiotest","",nBinEtaCorr,binLimEtaCorr);
 436   hRatiotest->GetXaxis()->SetTitle("#eta");
 437   hRatiotest->GetYaxis()->SetTitle("Generated/Corrected");
 438   hRatiotest->Divide(hdNdEta_test,hdNdEta_fullyCorr,1.,1.);
 439   new TCanvas();
 440   hRatiotest->Draw("p");
 441
 442   // Generated/corrected ratios for consistency checks
 443   TH2F* hRatiodNdEta2D=new TH2F("ratiodNdEta2D","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 444   hRatiodNdEta2D->Divide(hSPDEta_2D_fullyCorr,Eta,1.,1.);
 445   new TCanvas();
 446   hRatiodNdEta2D->Draw();
 447
 448   TH1F* hRatiodNdEta=new TH1F("ratiodNdEta","",nBinEtaCorr,binLimEtaCorr);
 449   hRatiodNdEta->GetXaxis()->SetTitle("#eta");
 450   hRatiodNdEta->GetYaxis()->SetTitle("Generated/Corrected");
 451   hRatiodNdEta->Divide(hdNdEta,hdNdEta_fullyCorr,1.,1.);
 452   new TCanvas();
 453   hRatiodNdEta->Draw("p");
 454
 455   // Draw dN/dEta
 456   new TCanvas();
 457 //  hdNdEta->SetLineWidth(3);
 458 //  hdNdEta->Draw("histo");
 459
 460   hdNdEta_raw->SetMinimum(0.);
 461   hdNdEta_raw->SetMaximum(3000);
 462   hdNdEta_raw->SetLineColor(kGreen);
 463   hdNdEta_raw->SetLineWidth(3);
 464   hdNdEta_raw->Draw("histo");
 465
 466   hdNdEta_bkgCorr->SetMarkerStyle(21);
 467   hdNdEta_bkgCorr->SetMarkerColor(kBlue);
 468
 469   hdNdEta_bkgCorr->Draw("same,p");
 470
 471   hdNdEta_fullyCorr->SetMarkerStyle(20);
 472   hdNdEta_fullyCorr->SetMarkerColor(kRed);
 473   hdNdEta_fullyCorr->Draw("p,same");
 474
 475   TLegend *leg1 = new TLegend(0.2,0.7,0.7,0.9,NULL,"brNDC");
 476   leg1->SetFillColor(0);
 477   leg1->SetTextFont(62);
 478   leg1->SetTextSize(0.0304569);
 479   TLegendEntry *entry1=leg1->AddEntry(hdNdEta_raw,"Raw","l");
 480 //                entry1=leg1->AddEntry(hdNdEta,"Generated","l");
 481                 entry1=leg1->AddEntry(hdNdEta_bkgCorr,"Comb bkg corrected","p");
 482                 entry1=leg1->AddEntry(hdNdEta_fullyCorr,"Eff/acc corrected","p");
 483
 484   leg1->Draw();
 485
 486 /*  new TCanvas();
 487   // plot the relative stat error for this correction
 488   TH2F* hStatErrTrackToPart = new TH2F("staterrperctracktopart","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
 489   hStatErrTrackToPart->GetXaxis()->SetTitle("#eta");
 490   hStatErrTrackToPart->GetYaxis()->SetTitle("z_{vtx} [cm]");
 491   hStatErrTrackToPart->GetZaxis()->SetTitle("Statistical error (%)");
 492   for (Int_t kvtx=0; kvtx<nBinVtx; kvtx++) {
 493     for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 494       if (hTrackToParticleCorr->GetBinContent(jeta+1,kvtx+1))
 495       hStatErrTrackToPart->SetBinContent(jeta+1,kvtx+1,(hTrackToParticleCorr->GetBinError(jeta+1,kvtx+1))/(hTrackToParticleCorr->GetBinContent(jeta+1,kvtx+1))*100);
 496       else hStatErrTrackToPart->SetBinContent(jeta+1,kvtx+1,0);
 497     }
 498   }
 499   new TCanvas();
 500   hStatErrTrackToPart->Draw();
 501   new TCanvas();
 502   hRatiodNdEta_trackToParticle->Draw();
 503   new TCanvas();
 504   hTrackToParticleCorr->Draw();*/
 505
 506   // Write histos
 507   TFile *fout= new TFile("SPDdNdEta.root","RECREATE");
 508
 509   hCombBkgCorrData->Write();
 510   hCombBkgCorrMC->Write();
 511   hCombBkgCorr2MC->Write();
 512   hAlphaCorr->Write();
 513
 514   hdNdEta_raw->Write();
 515   hdNdEta_bkgCorr->Write();
 516   hdNdEta_fullyCorr->Write();
 517
 518   hCorrecteddNdEtaCentr->Write();
 519
 520   hEffProtons->Write();
 521   hEffKaons->Write();
 522   hEffPions->Write();
 523   hEffOthers->Write();
 524
 525   fout->Write();
 526   fout->Close();
 527
 528 }
 529
 530 Float_t CalculateErrordNdEtaPerPartPair(Float_t a, Float_t b, Float_t da, Float_t db) {
 531
 532   Float_t errdndetaperpartpairs = (a/b)*TMath::Sqrt(TMath::Power(da,2)/TMath::Power(a,2)+TMath::Power(db,2)/TMath::Power(b,2));
 533   return errdndetaperpartpairs;
 534
 535 }
 536
 537 void combscalingfactors (TList *lall, TList *lcomb, TList *lmc, Double_t *scaleNormRange_rot, Double_t *scaleNormRange_labels) {
 538
 539
 540  TH2F *hall = (TH2F*) lall->FindObject("fHistSPDdePhideTheta");
 541  TH2F *hcomb = (TH2F*) lcomb->FindObject("fHistSPDdePhideTheta");
 542
 543  TH2F *hlabel = (TH2F*) lmc->FindObject("fHistdPhidThetaComb");
 544
 545  TH1D* hproall = new TH1D("_xall","",2000,-1.,1.);
 546  hproall->Sumw2();
 547  hproall = hall->ProjectionX("_xa",0,-1,"e");
 548
 549  TH1D* hprocomb = new TH1D("_xcomb","",2000,-1.,1.);
 550  hprocomb->Sumw2();
 551  hprocomb = hcomb->ProjectionX("_xc",0,-1,"e");
 552
 553  TH1D *hprolabel = new TH1D("_xcomblabel","",2000,-1.,1.);
 554  hprolabel = hlabel->ProjectionX("_xcombl",0,-1,"e");
 555
 556  // Normalize to integral in [0.08,0.20]
 557  Double_t denNormRange = hproall->Integral(801,920) + hproall->Integral(1081,1200);
 558  Double_t numNormRange_rot   = hprocomb->Integral(801,920) + hprocomb->Integral(1081,1200);
 559  *scaleNormRange_rot = denNormRange / numNormRange_rot;
 560
 561  new TCanvas();
 562  hproall->Draw("histo");
 563
 564  hprocomb->Scale(*scaleNormRange_rot);
 565  cout<<"Scaling factor normalizing to close tails: "<<*scaleNormRange_rot<<endl;
 566  cout<<"Percentage of background"<< hprocomb->Integral(1,2000)/hproall->Integral(1,2000)<<endl;
 567  cout<<"Percentage of background in |#D#phi|<0.08 "<< hprocomb->Integral(921,1080)/hproall->Integral(921,1080)<<endl;
 568  cout<<"Percentage of background in |#D#phi|<0.06 "<< hprocomb->Integral(941,1060)/hproall->Integral(941,1060)<<endl;
 569  hprocomb->SetLineColor(kBlue);
 570  hprocomb->Draw("histo,same");
 571
 572  // Fit the label bkg comb distribution to the data distribution
 573  // Normalize to integral in [0.08,0.20]
 574  Double_t numNormRange_labels   = hprolabel->Integral(801,920) + hprolabel->Integral(1081,1200);
 575  *scaleNormRange_labels = denNormRange / numNormRange_labels;
 576
 577  hprolabel->Scale(*scaleNormRange_labels);
 578  cout<<"Scaling factor labels: "<<*scaleNormRange_labels<<endl;
 579  cout<<"Percentage of background with labels "<< (hprolabel->Integral(1,2000)/hproall->Integral(1,2000))<<endl;
 580  cout<<"Percentage of background with labels in |#D#phi|<0.08 "<< (hprolabel->Integral(921,1080)/hproall->Integral(921,1080))<<endl;
 581
 582  hprolabel->SetLineColor(kGreen);
 583  hproall->Draw("histo,same");
 584 /*
 585  TFile *fout= new TFile("scaling.root","RECREATE");
 586  hproall->Write();
 587  hprocomb->Write();
 588  hprolabel->Write();
 589  fout->Write();
 590  fout->Close();
 591 */
 592 }
 593
 594