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