Cosmetics
[u/mrichter/AliRoot.git] / 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 #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