1 /*************************************************************************
2 * Macro correctSPDdNdEta *
3 * To read data and correction histograms produced with *
4 * AliAnalysisTaskSPDdNdEta and apply corrections to data *
6 * Author: M. Nicassio (INFN Bari) *
7 * Contact: Maria.Nicassio@ba.infn.it, Domenico.Elia@ba.infn.it *
8 **************************************************************************/
10 #include "Riostream.h"
20 #include "TLegendEntry.h"
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);
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) {
33 // vtx lim 13-27 for +-7cm
34 // vtx lim 15-25 for +-5cm
38 Double_t MaxVtx = 20.;
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;
49 const Int_t nBinEtaCorr = 60;
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;
57 Float_t nBinsPerPseudorapidityUnit = nBinEtaCorr/(binLimEtaCorr[nBinEtaCorr]-binLimEtaCorr[0]);
58 cout<<"Number of bins per pseudorapidity unit"<<nBinsPerPseudorapidityUnit<<endl;
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);
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);
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");
90 Double_t scaleBkg = 1.;
91 Double_t scaleBkgMC = 1.;
92 Double_t scaleBkgLabels = 1.;
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);
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");
103 //Corrections at track level
104 // Combinatorial background: beta correction from data
105 TH2F *hCombBkg = (TH2F*)l_rawbkg->FindObject("fHistSPDRAWEtavsZ");
106 hCombBkg->Scale(scaleBkg);
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();
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");
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();
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);
129 hCombBkgCorrData->Divide(hCombBkg,hSPDEta_2Draw,1.,1.);
133 TH2F *hCombBkgCorr2Data = new TH2F("backgroundCorr2Data","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
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));
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));
151 // Alpha correction: MC only
152 TH2F *hPrimaries_evtTrVtx = new TH2F(*((TH2F*) l_acorr->FindObject("fHistNonDetectableCorrNum"))); // all prim in sel events
154 hPrimaries_evtTrVtx->Draw();
156 TH2F *hInvAlphaCorr = new TH2F("InvAlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
157 hInvAlphaCorr->Sumw2();
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")));
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
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")));
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);
187 TH2F *hCombBkgCorr2MC = new TH2F("backgroundCorr2MC","Combinatorial background correction",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
189 if (changeStrangeness) {
190 if (bkgFromMCLabels) hCombBkgCorrMC->Divide(hBkgCombLabels,hRecReweighted,1.,1.);
191 else hCombBkgCorrMC->Divide(hCombBkgMC,hRecReweighted,1.,1.);
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));
197 hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hRecReweighted);
198 hInvAlphaCorr->Divide(hPrimReweighted);
201 if (bkgFromMCLabels) hCombBkgCorrMC->Divide(hBkgCombLabels,hCombBkgMCDen,1.,1.);
202 else hCombBkgCorrMC->Divide(hCombBkgMC,hCombBkgMCDen,1.,1.);
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));
208 hInvAlphaCorr->Multiply(hCombBkgCorr2MC,hCombBkgMCDen);
209 hInvAlphaCorr->Divide(hPrimaries_evtTrVtx);
214 // hCombBkgCorr2MC->Draw();
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);
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.);
228 // hEffKaons->Divide(hEffPions);
232 hInvAlphaCorr->Draw();
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));
242 TH2F *hAlphaCorr = new TH2F("AlphaCorr","",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
243 hAlphaCorr->GetXaxis()->SetTitle("#eta");
244 hAlphaCorr->GetYaxis()->SetTitle("z_{vtx} [cm]");
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));
256 //////////////////////// Factorize alpha ///////////////////////////
258 ////////////////////////////////////////////////////////////////////
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;
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.);
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.);
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.);
286 hCombBkgCorrMC->Draw();
288 hCombBkgCorrData->Draw();
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
298 TH2F *Eta = (TH2F*)l_MC->FindObject("fHistNonDetectableCorrNum");
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);
304 Double_t nEvtsTotSel= hMCvertex->Integral(binVtxStart+1,binVtxStop);
306 cout<<"Number of generated events: "<<nEvtsTot<<endl;
307 cout<<"Number of generated events in the selected vertex range: "<<nEvtsTotSel <<endl;
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.);
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;
321 //Corrected distributions
322 TH2F *hSPDEta_2D_bkgCorr = new TH2F("SPDEta_2D_bkgCorr", "",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
323 hSPDEta_2D_bkgCorr->Sumw2();
325 TH2F *hSPDEta_2D_fullyCorr = new TH2F("SPDEta_2D_fullyCorr", "",nBinEtaCorr,binLimEtaCorr,nBinVtx,-MaxVtx,MaxVtx);
326 hSPDEta_2D_fullyCorr->Sumw2();
328 TH1F *hnorm_fullyCorr = new TH1F("norm_fullyCorr", "",nBinEtaCorr,binLimEtaCorr);
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);
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();
345 //Apply corrections at
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);
352 hSPDEta_2D_bkgCorr->Draw();
354 hSPDEta_2D_fullyCorr->Draw();
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);
366 for (Int_t jeta=0; jeta<nBinEtaCorr; jeta++) {
367 hnorm_fullyCorr->SetBinError(jeta+1,0);
370 hnorm_fullyCorr->Draw();
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);
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));
382 hCorrecteddNdEtaCentr->Rebin(10);
385 hCorrecteddNdEtaCentr->Draw();
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.;
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
394 Float_t NpartV0Glauber = 381.188;
395 Float_t NpartCl2Glauber = 378.5;
396 Float_t NpartCl2Hij = 365.3;
398 Float_t errNpartV0Glauber = 18.4951;
399 Float_t errNpartCl2Glauber = 7.57; // 2% not used
400 Float_t errNpartCl2Hij = 7.306;
402 Float_t dNdEtaCorr = ncorrtrackletscentr/nevtscentr;
403 Float_t errdNdEtaCorr = hCorrecteddNdEtaCentr->GetBinError(1);
405 Float_t dNdEtaPartPairsV0G = dNdEtaCorr/(.5*NpartV0Glauber);
406 Float_t dNdEtaPartPairsCl2G = dNdEtaCorr/(.5*NpartCl2Glauber);
407 Float_t dNdEtaPartPairsCl2H = dNdEtaCorr/(.5*NpartCl2Hij);
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);
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;
417 hdNdEta_bkgCorr->Scale(nBinsPerPseudorapidityUnit);
418 hdNdEta_fullyCorr->Scale(nBinsPerPseudorapidityUnit);
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.);
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.);
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.);
440 hRatiotest->Draw("p");
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.);
446 hRatiodNdEta2D->Draw();
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.);
453 hRatiodNdEta->Draw("p");
457 // hdNdEta->SetLineWidth(3);
458 // hdNdEta->Draw("histo");
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");
466 hdNdEta_bkgCorr->SetMarkerStyle(21);
467 hdNdEta_bkgCorr->SetMarkerColor(kBlue);
469 hdNdEta_bkgCorr->Draw("same,p");
471 hdNdEta_fullyCorr->SetMarkerStyle(20);
472 hdNdEta_fullyCorr->SetMarkerColor(kRed);
473 hdNdEta_fullyCorr->Draw("p,same");
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");
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);
500 hStatErrTrackToPart->Draw();
502 hRatiodNdEta_trackToParticle->Draw();
504 hTrackToParticleCorr->Draw();*/
507 TFile *fout= new TFile("SPDdNdEta.root","RECREATE");
509 hCombBkgCorrData->Write();
510 hCombBkgCorrMC->Write();
511 hCombBkgCorr2MC->Write();
514 hdNdEta_raw->Write();
515 hdNdEta_bkgCorr->Write();
516 hdNdEta_fullyCorr->Write();
518 hCorrecteddNdEtaCentr->Write();
520 hEffProtons->Write();
530 Float_t CalculateErrordNdEtaPerPartPair(Float_t a, Float_t b, Float_t da, Float_t db) {
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;
537 void combscalingfactors (TList *lall, TList *lcomb, TList *lmc, Double_t *scaleNormRange_rot, Double_t *scaleNormRange_labels) {
540 TH2F *hall = (TH2F*) lall->FindObject("fHistSPDdePhideTheta");
541 TH2F *hcomb = (TH2F*) lcomb->FindObject("fHistSPDdePhideTheta");
543 TH2F *hlabel = (TH2F*) lmc->FindObject("fHistdPhidThetaComb");
545 TH1D* hproall = new TH1D("_xall","",2000,-1.,1.);
547 hproall = hall->ProjectionX("_xa",0,-1,"e");
549 TH1D* hprocomb = new TH1D("_xcomb","",2000,-1.,1.);
551 hprocomb = hcomb->ProjectionX("_xc",0,-1,"e");
553 TH1D *hprolabel = new TH1D("_xcomblabel","",2000,-1.,1.);
554 hprolabel = hlabel->ProjectionX("_xcombl",0,-1,"e");
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;
562 hproall->Draw("histo");
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");
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;
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;
582 hprolabel->SetLineColor(kGreen);
583 hproall->Draw("histo,same");
585 TFile *fout= new TFile("scaling.root","RECREATE");