PWG2/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 To do
  11 1. test prim
  12 2. check errors on corrections
  13 3. trigger efficiency correction
  14 */
  15
  16
  17 #include "Riostream.h"
  18 #include "TROOT.h"
  19 #include "TStyle.h"
  20 #include "TMath.h"
  21 #include "TH1.h"
  22 #include "TH2.h"
  23 #include "TFile.h"
  24 #include "TList.h"
  25 #include "TTree.h"
  26 #include "TLegend.h"
  27 #include "TLegendEntry.h"
  28 #include "TCanvas.h"
  29
  30 void correctSPDdNdEta (Char_t* fileRaw, Char_t* fileRawBkgCorr, Char_t* fileAlphaCorr, Char_t* fileMCBkgCorr, Char_t* fileMC,
  31                        Bool_t kPrimariesTest=kFALSE,
  32                        Int_t binVtxStart=5, Int_t binVtxStop=15,
  33                        Double_t cutCorr=0., Double_t cutBkg=1., Double_t cutBkgMC=1., Float_t scaleBkg=0.255, Float_t scaleBkgMC=.255)   {
  34
  35 //lim vtx: 5-15 ->[-10cm,10cm[; 0-20 ->[-20cm,20cm[
  36
  37   const Int_t nBinMultCorr = 200;
  38   Float_t nMaxMult = 20000.;
  39   Double_t binLimMultCorr[nBinMultCorr+1];
  40   binLimMultCorr[0] = 0.;
  41   binLimMultCorr[1] = 1.;
  42   for (Int_t i = 2; i<=nBinMultCorr;++i) {
  43     binLimMultCorr[i] = (i-1)*nMaxMult/nBinMultCorr;
  44   }
  45
  46   const Int_t nBinEtaCorr = 60;
  47   Float_t etaMax = 3.;
  48   Float_t etaMin = -3.;
  49   Double_t *binLimEtaCorr = new Double_t[nBinEtaCorr+1];
  50   for (Int_t i = 0; i<nBinEtaCorr+1; ++i) {
  51     binLimEtaCorr[i] = (Double_t) etaMin+i*(etaMax*2.)/nBinEtaCorr;
  52   }
  53
  54   Float_t nBinsPerPseudorapidityUnit = nBinEtaCorr/(binLimEtaCorr[nBinEtaCorr]-binLimEtaCorr[0]);
  55   cout<<"Number of bins per pseudorapidity unit"<<nBinsPerPseudorapidityUnit<<endl;
  56
  57   gStyle->SetOptLogy(kFALSE);
  58   gStyle->SetFrameLineWidth(2);
  59   gStyle->SetFrameFillColor(0);
  60   gStyle->SetCanvasColor(0);
  61   gStyle->SetPadColor(0);
  62   gStyle->SetHistLineWidth(2);
  63   gStyle->SetLabelSize(0.05, "x");
  64   gStyle->SetLabelSize(0.05, "y");
  65   gStyle->SetTitleSize(0.05, "x");
  66   gStyle->SetTitleSize(0.05, "y");
  67   gStyle->SetTitleOffset(1.1, "x");
  68   gStyle->SetTitleOffset(0.9, "y");
  69   gStyle->SetPadBottomMargin(0.14);
  70   gStyle->SetPalette(1,0);
  71   gStyle->SetTitleFillColor(0);
  72   gStyle->SetStatColor(0);
  73
  74   //Input files
  75   TFile *f_MC = new TFile(fileMC);
  76   TFile *f_acorr = new TFile(fileAlphaCorr);
  77   TFile *f_mcbkgcorr = new TFile(fileMCBkgCorr);
  78   TFile *f_rawbkg = new TFile(fileRawBkgCorr);
  79   TFile *f_raw = new TFile(fileRaw);
  80
  81   TList *l_MC = (TList*)f_MC->Get("cOutput");
  82   TList *l_acorr = (TList*)f_acorr->Get("cOutput");
  83   TList *l_mcbkgcorr = (TList*)f_mcbkgcorr->Get("cOutput");
  84   TList *l_rawbkg = (TList*)f_rawbkg->Get("cOutput");
  85   TList *l_raw = (TList*)f_raw->Get("cOutput");
  86
  87   //Histogram to be corrected at tracklet level
  88   TH2F *hSPDEta_2Draw = new TH2F(*((TH2F*)l_raw->FindObject("fHistSPDRAWEtavsZ")));
  89   hSPDEta_2Draw->SetNameTitle("SPDEta_2Draw","Reconstructed tracklets");
  90
  91   //Corrections at track level
  92   // Combinatorial background: beta correction from data
  93   TH2F *hCombBkg = (TH2F*)l_rawbkg->FindObject("fHistSPDRAWEtavsZ");
  94   hCombBkg->Scale(scaleBkg);
  95
  96   TH2F *hCombBkgCorrData = new TH2F("backgroundCorrDATA","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,20,-20.,20.);
  97   hCombBkgCorrData->GetXaxis()->SetTitle("#eta");
  98   hCombBkgCorrData->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
  99   hCombBkgCorrData->Sumw2();
 100   hCombBkgCorrData->Divide(hCombBkg,hSPDEta_2Draw,1.,1.);
 101
 102   // Combinatorial background: beta correction from MC
 103   TH2F *hCombBkgMC = (TH2F*)l_mcbkgcorr->FindObject("fHistSPDRAWEtavsZ");
 104   hCombBkgMC->Scale(scaleBkgMC);
 105   TH2F *hCombBkgMCDen = (TH2F*)l_acorr->FindObject("fHistSPDRAWEtavsZ");
 106
 107   TH2F *hCombBkgCorrMC = new TH2F("backgroundCorrMC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,20,-20.,20.);
 108   hCombBkgCorrMC->GetXaxis()->SetTitle("#eta");
 109   hCombBkgCorrMC->GetYaxis()->SetTitle("z_{SPDvtx} [cm]");
 110   hCombBkgCorrMC->Sumw2();
 111   hCombBkgCorrMC->Divide(hCombBkgMC,hCombBkgMCDen,1.,1.);
 112
 113   // 1-beta in MC
 114   TH2F *hCombBkgCorr2MC = new TH2F("backgroundCorr2MC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,20,-20.,20.);
 115   for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++)
 116     for (Int_t jz=0; jz<20; jz++)
 117       hCombBkgCorr2MC->SetBinContent(ieta+1,jz+1,1.);
 118   hCombBkgCorr2MC->Add(hCombBkgCorrMC,-1.);
 119
 120 //  new TCanvas();
 121 //  hCombBkgCorr2MC->Draw();
 122
 123   // Errors on beta
 124   //..data
 125   for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
 126     for (Int_t jz=0; jz<20; jz++) {
 127       hCombBkgCorrData->SetBinError(ieta+1,jz+1,scaleBkg*TMath::Sqrt(hCombBkg->GetBinContent(ieta+1,jz+1))/hSPDEta_2Draw->GetBinContent(ieta+1,jz+1));
 128     }
 129   }
 130
 131   //..MC--> no, computed in alpha
 132
 133
 134   // Alpha correction: MC only
 135   // alpha here is 1/alpha
 136   TH2F *hPrimaries_evtTrVtx = new TH2F(*((TH2F*) l_acorr->FindObject("fHistNonDetectableCorrNum"))); // all prim in sel events
 137
 138   TH2F *hAlphaCorr = new TH2F("AlphaCorr","",nBinEtaCorr,binLimEtaCorr,20,-20.,20.);
 139   hAlphaCorr->GetXaxis()->SetTitle("#eta");
 140   hAlphaCorr->GetYaxis()->SetTitle("z_{vtx} [cm]");
 141 //  hAlphaCorr->Divide(hPrimaries_evtTrVtx,hCombBkgMCDen,1.,1.);
 142 //  hAlphaCorr->Divide(hCombBkgCorr2MC);
 143   hAlphaCorr->Multiply(hCombBkgCorr2MC,hCombBkgMCDen);
 144   hAlphaCorr->Divide(hPrimaries_evtTrVtx);
 145
 146     // Errors on alpha
 147
 148   for (Int_t ieta=0; ieta<nBinEtaCorr; ieta++) {
 149     for (Int_t jz=0; jz<20; jz++) {
 150       hAlphaCorr->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));
 151     }
 152   }
 153
 154
 155
 156   // Number of events and normalization histogram
 157   TH1D *hSPDvertex = (TH1D*)l_raw->FindObject("fHistSPDvtxAnalysis");
 158   Double_t nEvtsRec=hSPDvertex->Integral(binVtxStart+1,binVtxStop);
 159   cout<<"Number of reconstructed events in the selected vertex range: "<<nEvtsRec <<endl;
 160
 161   //Cutting corrections at the edges of the acceptance
 162   for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 163     for (Int_t kvtx=0; kvtx<20; kvtx++) {
 164       if (hAlphaCorr->GetBinContent(jeta+1,kvtx+1)<cutCorr) {
 165         hAlphaCorr->SetBinContent(jeta+1,kvtx+1,0.);
 166         hAlphaCorr->SetBinError(jeta+1,kvtx+1,0.);
 167       }
 168       if (hCombBkgCorrData->GetBinContent(jeta+1,kvtx+1)>cutBkg) {
 169         hCombBkgCorrData->SetBinContent(jeta+1,kvtx+1,0.);
 170         hCombBkgCorrData->SetBinError(jeta+1,kvtx+1,0.);
 171       }
 172       if (hCombBkgCorrMC->GetBinContent(jeta+1,kvtx+1)>cutBkgMC) {
 173         hCombBkgCorrMC->SetBinContent(jeta+1,kvtx+1,0.);
 174         hCombBkgCorrMC->SetBinError(jeta+1,kvtx+1,0.);
 175       }
 176     }
 177   }
 178
 179   new TCanvas();
 180   hAlphaCorr->Draw();
 181   new TCanvas();
 182   hCombBkgCorrMC->Draw();
 183   new TCanvas();
 184   hCombBkgCorrData->Draw();
 185
 186
 187   // MC distribution --> ok if running on samples for centrality bins
 188
 189   TH2F *hVertexMC2D = (TH2F*)l_MC->FindObject("fHistSelEvts");  // pay att here!! change if comparing to true MC centrality
 190   TH1D *hMCvertex = new TH1D(*(hVertexMC2D->ProjectionY("MCvertex",0,-1,"e")));  //MC vertex distrib all events
 191   TH2F *Eta = (TH2F*)l_MC->FindObject("fHistNonDetectableCorrNum"); // test comparing with MC in same subset (centrality selected)
 192
 193 //  TH2F *Eta = (TH2F*)l_MC->FindObject("fHistAllPrimaries"); // to compare with MC selection!!
 194   TH1D *hdNdEta = new TH1D("dNdEta","Pseudorapidity ",nBinEtaCorr,binLimEtaCorr);
 195   hdNdEta->Sumw2();
 196   hdNdEta = Eta->ProjectionX("dNdEta",0,-1,"e");
 197   Double_t nEvtsTot=hMCvertex->Integral(0,hMCvertex->GetNbinsX()+1);
 198   Double_t nEvtsTot10= hMCvertex->Integral(binVtxStart+1,binVtxStop);
 199
 200   cout<<"Number of generated events in +-20cm: "<<nEvtsTot <<endl;
 201   cout<<"Number of generated events in the selected vertex range: "<<nEvtsTot10 <<endl;
 202
 203   hdNdEta->Scale(nBinsPerPseudorapidityUnit/nEvtsTot);
 204   hdNdEta->GetXaxis()->SetTitle("#eta");
 205   hdNdEta->GetYaxis()->SetTitle("dN_{ch}/d#eta");
 206   new TCanvas();
 207   hdNdEta->Draw();
 208
 209   //Corrected distributions
 210   TH2F *hSPDEta_2D_bkgCorr =  new TH2F("SPDEta_2D_bkgCorr", "",nBinEtaCorr,binLimEtaCorr,20,-20.,20.);
 211   hSPDEta_2D_bkgCorr->Sumw2();
 212
 213   TH2F *hSPDEta_2D_fullyCorr =  new TH2F("SPDEta_2D_fullyCorr", "",nBinEtaCorr,binLimEtaCorr,20,-20.,20.);
 214   hSPDEta_2D_fullyCorr->Sumw2();
 215
 216   TH1F *hnorm_fullyCorr =  new TH1F("norm_fullyCorr", "",nBinEtaCorr,binLimEtaCorr);
 217   hnorm_fullyCorr->Sumw2();
 218
 219
 220   //dN/deta
 221   TH1F *hdNdEta_raw = new TH1F("dNdEta_raw","",nBinEtaCorr,binLimEtaCorr);
 222   hdNdEta_raw->Sumw2();
 223   hdNdEta_raw->GetXaxis()->SetTitle("#eta");
 224   hdNdEta_raw->GetYaxis()->SetTitle("dN_{ch}/#eta");
 225   hdNdEta_raw->Add(hSPDEta_2Draw->ProjectionX("SPDEta_2Draw_x",binVtxStart+1,binVtxStop,"e"));
 226   hdNdEta_raw->Scale(nBinsPerPseudorapidityUnit/nEvtsRec);
 227   new TCanvas();
 228   hdNdEta_raw->Draw();
 229
 230   TH1F *hdNdEta_bkgCorr = new TH1F("dNdEta_bkgCorr","",nBinEtaCorr,binLimEtaCorr);
 231   hdNdEta_bkgCorr->Sumw2();
 232   TH1F *hdNdEta_fullyCorr = new TH1F("dNdEta_fullyCorr","",nBinEtaCorr,binLimEtaCorr);
 233   hdNdEta_fullyCorr->Sumw2();
 234
 235   //Apply corrections at
 236   //...track level
 237   hSPDEta_2D_bkgCorr->Add(hSPDEta_2Draw,-1);
 238   hSPDEta_2D_bkgCorr->Multiply(hCombBkgCorrData);
 239   hSPDEta_2D_bkgCorr->Add(hSPDEta_2Draw);
 240   hSPDEta_2D_fullyCorr->Add(hSPDEta_2D_bkgCorr);
 241   hSPDEta_2D_fullyCorr->Divide(hAlphaCorr);
 242   new TCanvas();
 243   hSPDEta_2D_bkgCorr->Draw();
 244   new TCanvas();
 245   hSPDEta_2D_fullyCorr->Draw();
 246
 247   //Filling normalization histogram
 248   for (Int_t ivtx=binVtxStart; ivtx<binVtxStop; ivtx++) {
 249     Double_t nEvtsivtx = hSPDvertex->GetBinContent(ivtx+1);
 250     for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 251       if (hAlphaCorr->GetBinContent(jeta+1,ivtx+1)!=0.) {
 252         hnorm_fullyCorr->Fill(hAlphaCorr->GetXaxis()->GetBinCenter(jeta+1),nEvtsivtx);
 253       }
 254     }
 255   }
 256
 257   for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 258     hnorm_fullyCorr->SetBinError(jeta+1,0);
 259   }
 260   new TCanvas();
 261   hnorm_fullyCorr->Draw();
 262
 263   hdNdEta_bkgCorr->Add(hSPDEta_2D_bkgCorr->ProjectionX("SPDEta_2D_bkgCorr_x",binVtxStart+1,binVtxStop,"e"));
 264   hdNdEta_bkgCorr->Scale(1./nEvtsRec);
 265   hdNdEta_fullyCorr->Divide(hSPDEta_2D_fullyCorr->ProjectionX("SPDEta_2D_fullyCorr_x",binVtxStart+1,binVtxStop,"e"),hnorm_fullyCorr);
 266
 267   hdNdEta_bkgCorr->Scale(nBinsPerPseudorapidityUnit);
 268   hdNdEta_fullyCorr->Scale(nBinsPerPseudorapidityUnit);
 269
 270   new TCanvas();
 271   hdNdEta_bkgCorr->Draw();
 272   new TCanvas();
 273   hdNdEta_fullyCorr->Draw();
 274   hdNdEta->Draw("same");
 275   // Create mask for MC prim in SPD acceptance
 276   TH1D *hdNdEta_test = new TH1D("dNdEta_test","Pseudorapidity ",nBinEtaCorr,binLimEtaCorr);
 277   TH2F *hMask = new TH2F("Mask","",nBinEtaCorr,binLimEtaCorr,20,-20,20);
 278   hMask->Sumw2();
 279   for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 280     for (Int_t kvtx=0; kvtx<20; kvtx++) {
 281        if (hSPDEta_2D_fullyCorr->GetBinContent(jeta+1,kvtx+1)!=0) hMask->SetBinContent(jeta+1,kvtx+1,1.);
 282
 283     }
 284   }
 285   hMask->Multiply(Eta);
 286   hdNdEta_test->Add(hMask->ProjectionX("Mask_x",binVtxStart+1,binVtxStop,"e"));
 287   hdNdEta_test->Scale(1./nEvtsTot10);
 288   hdNdEta_test->Scale(nBinsPerPseudorapidityUnit);
 289
 290   TH1F* hRatiotest=new TH1F("ratiotest","",nBinEtaCorr,binLimEtaCorr);
 291   hRatiotest->GetXaxis()->SetTitle("#eta");
 292   hRatiotest->GetYaxis()->SetTitle("Generated/Corrected");
 293
 294   for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta_test->SetBinError(i+1,0.);
 295   hRatiotest->Divide(hdNdEta_test,hdNdEta_fullyCorr,1.,1.);
 296   new TCanvas();
 297   hRatiotest->Draw("p,histo");
 298
 299
 300   // Generated/corrected ratios
 301   TH2F* hRatiodNdEta2D=new TH2F("ratiodNdEta2D","",nBinEtaCorr,binLimEtaCorr,20,-20,20);
 302   hRatiodNdEta2D->Divide(hSPDEta_2D_fullyCorr,Eta,1.,1.);
 303   new TCanvas();
 304   hRatiodNdEta2D->Draw();
 305
 306   TH1F* hRatiodNdEta=new TH1F("ratiodNdEta","",nBinEtaCorr,binLimEtaCorr);  // Add ratio for two corrected dN/dEta
 307   hRatiodNdEta->GetXaxis()->SetTitle("#eta");
 308   hRatiodNdEta->GetYaxis()->SetTitle("Generated/Corrected");
 309
 310   for (Int_t i=0;i<nBinEtaCorr;i++) hdNdEta->SetBinError(i+1,0.);
 311   hRatiodNdEta->Divide(hdNdEta,hdNdEta_fullyCorr,1.,1.);
 312   new TCanvas();
 313   hRatiodNdEta->Draw("p,histo");
 314   // Draw dN/dEta
 315   TCanvas *ccorrected = new TCanvas("c3","Corrected distributions");
 316   hdNdEta->SetLineWidth(3);
 317   hdNdEta->SetMinimum(0.);
 318   hdNdEta->SetMaximum(4000);
 319   hdNdEta->Draw("histo");
 320
 321   hdNdEta_raw->SetLineColor(kGreen);
 322   hdNdEta_raw->SetLineWidth(3);
 323   hdNdEta_raw->Draw("histo,same");
 324
 325   hdNdEta_bkgCorr->SetMarkerStyle(21);
 326   hdNdEta_bkgCorr->SetMarkerColor(kBlue);
 327
 328   hdNdEta_bkgCorr->Draw("same,p,histo");
 329
 330   hdNdEta_fullyCorr->SetMarkerStyle(20);
 331   hdNdEta_fullyCorr->SetMarkerColor(kRed);
 332   hdNdEta_fullyCorr->Draw("p,same,histo");
 333
 334   TLegend *leg1 = new TLegend(0.2,0.7,0.7,0.9,NULL,"brNDC");
 335   leg1->SetFillColor(0);
 336   leg1->SetTextFont(62);
 337   leg1->SetTextSize(0.0304569);
 338   TLegendEntry *entry1=leg1->AddEntry(hdNdEta_raw,"Reconstructed","l");
 339                 entry1=leg1->AddEntry(hdNdEta,"Generated","l");
 340                 entry1=leg1->AddEntry(hdNdEta_bkgCorr,"Bkg corrected","p");
 341                 entry1=leg1->AddEntry(hdNdEta_fullyCorr,"Corrected (triggered events with vertex)","p");
 342
 343   leg1->Draw();
 344 /*  new TCanvas();
 345   // plot the relative stat error for this correction
 346   TH2F* hStatErrTrackToPart = new TH2F("staterrperctracktopart","",nBinEtaCorr,binLimEtaCorr,20,-20.,20.);
 347   hStatErrTrackToPart->GetXaxis()->SetTitle("#eta");
 348   hStatErrTrackToPart->GetYaxis()->SetTitle("z_{vtx} [cm]");
 349   hStatErrTrackToPart->GetZaxis()->SetTitle("Statistical error (%)");
 350   for (Int_t kvtx=0; kvtx<20; kvtx++) {
 351     for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
 352       if (hTrackToParticleCorr->GetBinContent(jeta+1,kvtx+1))
 353       hStatErrTrackToPart->SetBinContent(jeta+1,kvtx+1,(hTrackToParticleCorr->GetBinError(jeta+1,kvtx+1))/(hTrackToParticleCorr->GetBinContent(jeta+1,kvtx+1))*100);
 354       else hStatErrTrackToPart->SetBinContent(jeta+1,kvtx+1,0);
 355     }
 356   }
 357   new TCanvas();
 358   hStatErrTrackToPart->Draw();
 359   new TCanvas();
 360   hRatiodNdEta_trackToParticle->Draw();
 361 new TCanvas();
 362 hTrackToParticleCorr->Draw();
 363   // Write histos
 364 */
 365   TFile *fout= new TFile("SPDdNdEta.root","RECREATE");
 366
 367
 368   // save more histos
 369   hCombBkgCorrData->Write();
 370   hCombBkgCorrMC->Write();
 371   hCombBkgCorr2MC->Write();
 372   hAlphaCorr->Write();
 373   fout->Write();
 374   fout->Close();
 375
 376 }