1 #if !defined(__CINT__) || defined(__MAKECINT__)
6 #include <TGraphErrors.h>
11 #include <TGridResult.h>
17 #include <TPaveText.h>
21 void RatioPlot(TH1F **num,TH1F **den,TString t1,TString t2,TString opt,Double_t Low[],Double_t Up[]);
22 void SetDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1);
23 void ResetOutOfRange(TH1F *histo, Int_t ipart, Double_t lowRange[3], Double_t upRange[3]);
25 void MakeCorrectedSpectraITSsaNsigma( TString period="LHC10d",
26 TString MCname="LHC10f6a",
37 Bool_t NormalizeData=1;
38 Bool_t NormalizeToINEL=0;
39 Bool_t SaveOnlyAsymm=0;
40 Bool_t CorrectOnlyForEfficiency=0;
42 Bool_t DCACorrection=1;
43 Bool_t GFCorrection=1;
46 TH1F *fHistPrimMCposBefEvSel[3];
47 TH1F *fHistPrimMCnegBefEvSel[3];
48 TH1F *fHistPrimMCpos[3];
49 TH1F *fHistPrimMCneg[3];
50 TH1F *hHistPosNSigma[3];
51 TH1F *hHistNegNSigma[3];
52 TH1F *hHistPosNSigmaMean[3];
53 TH1F *hHistNegNSigmaMean[3];
54 TH1F *hHistPosNSigmaPrim[3];
55 TH1F *hHistNegNSigmaPrim[3];
56 TH1F *hHistPosNSigmaPrimMean[3];
57 TH1F *hHistNegNSigmaPrimMean[3];
58 TH1F *hHistPosNSigmaPrimMC[3];
59 TH1F *hHistNegNSigmaPrimMC[3];
60 TH1F *hHistPosNSigmaPrimMCMean[3];
61 TH1F *hHistNegNSigmaPrimMCMean[3];
62 TH1F *hHistPosNSigmaMC[3];
63 TH1F *hHistNegNSigmaMC[3];
64 TH1F *hHistPosNSigmaMCMean[3];
65 TH1F *hHistNegNSigmaMCMean[3];
66 TH1F *hCorrFactPos[3];
67 TH1F *hCorrFactNeg[3];
68 TH1F *hCorrFactPosMean[3];
69 TH1F *hCorrFactNegMean[3];
76 TH1F *hHistPosNSigmaDATA[3];
77 TH1F *hHistNegNSigmaDATA[3];
78 TH1F *hHistPosNSigmaMeanDATA[3];
79 TH1F *hHistNegNSigmaMeanDATA[3];
82 TH1F *hITSsaPosMean[3];
83 TH1F *hITSsaNegMean[3];
84 TH1F *hITSsaRawPos[3];
85 TH1F *hITSsaRawNeg[3];
86 TH1F *hITSsaRawPosMean[3];
87 TH1F *hITSsaRawNegMean[3];
90 Double_t Low[3]={0.10,0.2,0.3};
91 Double_t Up[3]={0.4,0.5,0.7};
92 Double_t lowRange[3]={0.10,0.2,0.3};
93 Double_t upRange[3]={0.6,0.6,0.6};
97 const Int_t nbinspt=22;
98 Double_t xbins[nbinspt+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
100 TString fname_MC=Form("./%s/AnalysisResults.root",MCname.Data());
101 cout<<fname_MC.Data()<<endl;
103 TString foutMC=Form("outSpectraMC_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
104 TFile *outMC=new TFile(foutMC.Data(),"recreate");
109 if(CorrectOnlyForEfficiency){
110 foutDATA=Form("outSpectraData_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
112 foutDATA=Form("outCorrSpectraData_%s_%s_%.0fDCA.root",period.Data(),MCname.Data(),DCAcut);
115 TString fnameDATA=Form("%s/AnalysisResults.root",period.Data());
116 cout<<fnameDATA.Data()<<endl;
117 TFile *outDATA=new TFile(foutDATA.Data(),"recreate");
122 TString lnameMC=Form("clistITSsaMult%ito%i",lowmult,upmult);
123 Printf("\n------------------ READING MC %s -------------------\n",lnameMC.Data());
124 TFile *finMC = new TFile(fname_MC.Data());
125 TDirectory *dirFileMC=(TDirectory*)finMC->Get("PWG2SpectraITSsa");
126 TList * cOutputMC = (TList*)dirFileMC->Get(lnameMC.Data());
127 //cOutputMC->Print();
129 TH1F *fHistNEventsMC = (TH1F*)cOutputMC->FindObject("fHistNEvents");
130 Float_t normMC=(Float_t)1/fHistNEventsMC->GetBinContent(fHistNEventsMC->FindBin(1.));
131 printf("- NormalizationMC (NEvents): %.0f \n",fHistNEventsMC->GetBinContent(fHistNEventsMC->FindBin(1.)));
133 for(Int_t ipart=0;ipart<=2;ipart++){
134 fHistPrimMCposBefEvSel[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCposBefEvSel%i",ipart));//primaries generated (before event selection)
135 fHistPrimMCnegBefEvSel[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCnegBefEvSel%i",ipart));
136 fHistPrimMCpos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCpos%i",ipart)); //primaries generated (after event selection)
137 fHistPrimMCneg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCneg%i",ipart));
138 hHistPosNSigma[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigma%i",ipart));//NSigma recontructed, no MCtruth
139 hHistNegNSigma[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigma%i",ipart));
140 hHistPosNSigmaPrimMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMC%i",ipart));//NSigma recontructed, MCtruth
141 hHistNegNSigmaPrimMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMC%i",ipart));
142 hHistPosNSigmaMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMC%i",ipart));//NSigma recontructed, MCtruth
143 hHistNegNSigmaMC[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMC%i",ipart));
144 hHistPosNSigmaPrim[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrim%i",ipart));//NSigma recontructed, truth
145 hHistNegNSigmaPrim[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrim%i",ipart));
146 hHistPosNSigmaMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMean%i",ipart));//NSigmaMean recontructed, no MCtruth
147 hHistNegNSigmaMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMean%i",ipart));
148 hHistPosNSigmaPrimMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMCMean%i",ipart));//NSigmaMean recontructed, no MCtruth
149 hHistNegNSigmaPrimMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMCMean%i",ipart));
150 hHistPosNSigmaMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaMCMean%i",ipart));//NSigmaMean recontructed, no MCtruth
151 hHistNegNSigmaMCMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaMCMean%i",ipart));
152 hHistPosNSigmaPrimMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistPosNSigmaPrimMean%i",ipart));//NSigmaMean recontructed, no MCtruth
153 hHistNegNSigmaPrimMean[ipart]=(TH1F*)cOutputMC->FindObject(Form("hHistNegNSigmaPrimMean%i",ipart));
156 for(Int_t ipart=0;ipart<=2;ipart++){
157 hCorrFactPos[ipart]=new TH1F(Form("hCorrFactPos%i",ipart),Form("hCorrFactPos%i",ipart),nbinspt,xbins);
158 hCorrFactNeg[ipart]=new TH1F(Form("hCorrFactNeg%i",ipart),Form("hCorrFactNeg%i",ipart),nbinspt,xbins);
159 hCorrFactPosMean[ipart]=new TH1F(Form("hCorrFactPosMean%i",ipart),Form("hCorrFactPosMean%i",ipart),nbinspt,xbins);
160 hCorrFactNegMean[ipart]=new TH1F(Form("hCorrFactNegMean%i",ipart),Form("hCorrFactNegMean%i",ipart),nbinspt,xbins);
161 hEffPos[ipart]=new TH1F(Form("hEffPos%i",ipart),Form("hEffPos%i",ipart),nbinspt,xbins);
162 hEffNeg[ipart]=new TH1F(Form("hEffNeg%i",ipart),Form("hEffNeg%i",ipart),nbinspt,xbins);
163 hEffPosMean[ipart]=new TH1F(Form("hEffPosMean%i",ipart),Form("hEffPosMean%i",ipart),nbinspt,xbins);
164 hEffNegMean[ipart]=new TH1F(Form("hEffNegMean%i",ipart),Form("hEffNegMean%i",ipart),nbinspt,xbins);
167 for(Int_t ipart=0;ipart<=2;ipart++){
168 hHistPosNSigma[ipart]->Sumw2();
169 hHistNegNSigma[ipart]->Sumw2();
170 hHistPosNSigmaMean[ipart]->Sumw2();
171 hHistNegNSigmaMean[ipart]->Sumw2();
172 hHistPosNSigmaPrimMC[ipart]->Sumw2();
173 hHistNegNSigmaPrimMC[ipart]->Sumw2();
174 hHistPosNSigmaMC[ipart]->Sumw2();
175 hHistNegNSigmaMC[ipart]->Sumw2();
176 hHistPosNSigmaPrimMCMean[ipart]->Sumw2();
177 hHistNegNSigmaPrimMCMean[ipart]->Sumw2();
178 hHistPosNSigmaPrim[ipart]->Sumw2();
179 hHistNegNSigmaPrim[ipart]->Sumw2();
180 hHistPosNSigmaPrimMean[ipart]->Sumw2();
181 hHistNegNSigmaPrimMean[ipart]->Sumw2();
182 hHistPosNSigmaMCMean[ipart]->Sumw2();
183 hHistNegNSigmaMCMean[ipart]->Sumw2();
185 hCorrFactPos[ipart]->Divide(hHistPosNSigma[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
186 hCorrFactNeg[ipart]->Divide(hHistNegNSigma[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
187 hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
188 hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
190 // hCorrFactPos[ipart]->Divide(hHistPosNSigmaPrim[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
191 // hCorrFactNeg[ipart]->Divide(hHistNegNSigmaPrim[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
192 // hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaPrimMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
193 // hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaPrimMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
195 //hCorrFactPos[ipart]->Divide(hHistPosNSigmaPrimMC[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
196 //hCorrFactNeg[ipart]->Divide(hHistNegNSigmaPrimMC[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
197 //hCorrFactPosMean[ipart]->Divide(hHistPosNSigmaPrimMCMean[ipart],fHistPrimMCposBefEvSel[ipart],1,1,"B");
198 //hCorrFactNegMean[ipart]->Divide(hHistNegNSigmaPrimMCMean[ipart],fHistPrimMCnegBefEvSel[ipart],1,1,"B");
200 hCorrFactPos[ipart]->SetTitle(Form("Correction Factor NSigma Pos%i",ipart));
201 hCorrFactNeg[ipart]->SetTitle(Form("Correction Factor NSigma Neg%i",ipart));
202 hCorrFactPosMean[ipart]->SetTitle(Form("Corr Fact NSigma Pos%i",ipart));
203 hCorrFactNegMean[ipart]->SetTitle(Form("Corr Fact NSigma Neg%i",ipart));
204 hCorrFactPos[ipart]->SetName(Form("hCorrFactPos%i",ipart));
205 hCorrFactNeg[ipart]->SetName(Form("hCorrFactNeg%i",ipart));
206 hCorrFactPosMean[ipart]->SetName(Form("hCorrFactPosMean%i",ipart));
207 hCorrFactNegMean[ipart]->SetName(Form("hCorrFactNegMean%i",ipart));
208 SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hCorrFactPos[ipart]);
209 SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hCorrFactNeg[ipart]);
210 SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hCorrFactPosMean[ipart]);
211 SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hCorrFactNegMean[ipart]);
212 ResetOutOfRange(hCorrFactPos[ipart],ipart,lowRange,upRange);
213 ResetOutOfRange(hCorrFactNeg[ipart],ipart,lowRange,upRange);
214 ResetOutOfRange(hCorrFactPosMean[ipart],ipart,lowRange,upRange);
215 ResetOutOfRange(hCorrFactNegMean[ipart],ipart,lowRange,upRange);
217 hHistPosNSigmaPrimMC[ipart]->Sumw2();
218 hHistNegNSigmaPrimMC[ipart]->Sumw2();
219 hHistPosNSigmaPrimMCMean[ipart]->Sumw2();
220 hHistNegNSigmaPrimMCMean[ipart]->Sumw2();
221 hEffPos[ipart]->Divide(hHistPosNSigmaPrimMC[ipart],fHistPrimMCpos[ipart],1,1,"B");
222 hEffNeg[ipart]->Divide(hHistNegNSigmaPrimMC[ipart],fHistPrimMCneg[ipart],1,1,"B");
223 hEffPosMean[ipart]->Divide(hHistPosNSigmaPrimMCMean[ipart],fHistPrimMCpos[ipart],1,1,"B");
224 hEffNegMean[ipart]->Divide(hHistNegNSigmaPrimMCMean[ipart],fHistPrimMCneg[ipart],1,1,"B");
225 hEffPos[ipart]->SetTitle(Form("Efficiency NSigma Pos%i",ipart));
226 hEffNeg[ipart]->SetTitle(Form("Efficiency NSigma Neg%i",ipart));
227 hEffPosMean[ipart]->SetTitle(Form("Efficiency NSigma Asymm Pos%i",ipart));
228 hEffNegMean[ipart]->SetTitle(Form("Efficiency NSigma Asymm Neg%i",ipart));
229 hEffPos[ipart]->SetName(Form("hEffPos%i",ipart));
230 hEffNeg[ipart]->SetName(Form("hEffNeg%i",ipart));
231 hEffPosMean[ipart]->SetName(Form("hEffPosMean%i",ipart));
232 hEffNegMean[ipart]->SetName(Form("hEffNegMean%i",ipart));
234 SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hEffPos[ipart]);
235 SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hEffNeg[ipart]);
236 SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hEffPosMean[ipart]);
237 SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hEffNegMean[ipart]);
238 ResetOutOfRange(hEffPos[ipart],ipart,lowRange,upRange);
239 ResetOutOfRange(hEffNeg[ipart],ipart,lowRange,upRange);
240 ResetOutOfRange(hEffPosMean[ipart],ipart,lowRange,upRange);
241 ResetOutOfRange(hEffNegMean[ipart],ipart,lowRange,upRange);
244 hCorrFactPos[ipart]->SetMinimum(0.);
245 hCorrFactNeg[ipart]->SetMinimum(0.);
246 hCorrFactPosMean[ipart]->SetMinimum(0.);
247 hCorrFactNegMean[ipart]->SetMinimum(0.);
248 hEffPos[ipart]->SetMinimum(0.);
249 hEffNeg[ipart]->SetMinimum(0.);
250 hEffPosMean[ipart]->SetMinimum(0.);
251 hEffNegMean[ipart]->SetMinimum(0.);
252 hCorrFactPos[ipart]->SetMaximum(1.);
253 hCorrFactNeg[ipart]->SetMaximum(1.);
254 hCorrFactPosMean[ipart]->SetMaximum(1.);
255 hCorrFactNegMean[ipart]->SetMaximum(1.);
256 hEffPos[ipart]->SetMaximum(1.);
257 hEffNeg[ipart]->SetMaximum(1.);
258 hEffPosMean[ipart]->SetMaximum(1.);
259 hEffNegMean[ipart]->SetMaximum(1.);
263 TCanvas *cCorr=new TCanvas("cCorr","cCorr");
265 for(Int_t ipart=0;ipart<=2;ipart++){
266 if(ipart==0)hCorrFactPos[ipart]->Draw("P");
267 else hCorrFactPos[ipart]->Draw("PSAME");
268 hCorrFactNeg[ipart]->Draw("PSAME");
270 cCorr->BuildLegend();
272 TCanvas *cCorrMean=new TCanvas("cCorrMean","cCorrMean");
273 cCorrMean->SetGridy();
274 for(Int_t ipart=0;ipart<=2;ipart++){
276 TH1F * hempty=(TH1F*)hCorrFactPosMean[ipart]->Clone();
277 hempty->SetXTitle("P_{T} (GeV/c)");
278 hempty->SetYTitle("Acc x #epsilon x Contamination");
279 hempty->Reset("all");
283 hCorrFactPosMean[ipart]->Draw("PSAME");
284 hCorrFactNegMean[ipart]->Draw("PSAME");
286 cCorrMean->BuildLegend();
288 TCanvas *cCorrcfPos=new TCanvas("cCorrcfPos","cCorrcfPos");
289 cCorrcfPos->SetGridy();
290 for(Int_t ipart=0;ipart<=2;ipart++){
291 if(ipart==0)hCorrFactPos[ipart]->Draw("P");
292 else hCorrFactPos[ipart]->Draw("PSAME");
293 hCorrFactPosMean[ipart]->Draw("PSAME");
295 cCorrcfPos->BuildLegend();
297 TCanvas *cCorrcfNeg=new TCanvas("cCorrcfNeg","cCorrcfNeg");
298 cCorrcfNeg->SetGridy();
299 for(Int_t ipart=0;ipart<=2;ipart++){
300 if(ipart==0)hCorrFactNeg[ipart]->Draw("P");
301 else hCorrFactNeg[ipart]->Draw("PSAME");
302 hCorrFactNegMean[ipart]->Draw("PSAME");
304 cCorrcfNeg->BuildLegend();
306 TCanvas *cEff=new TCanvas("cEff","cEff");
308 for(Int_t ipart=0;ipart<=2;ipart++){
309 if(ipart==0)hEffPos[ipart]->Draw("P");
310 else hEffPos[ipart]->Draw("PSAME");
311 hEffNeg[ipart]->Draw("PSAME");
315 TCanvas *cEffMean=new TCanvas("cEffMean","cEffMean");
316 cEffMean->SetGridy();
317 for(Int_t ipart=0;ipart<=2;ipart++){
318 if(ipart==0)hEffPosMean[ipart]->Draw("P");
319 else hEffPosMean[ipart]->Draw("PSAME");
320 hEffNegMean[ipart]->Draw("PSAME");
322 cEffMean->BuildLegend();
324 RatioPlot(hEffNegMean,hEffPosMean,"Neg","Pos","Efficiency",Low,Up);
327 TCanvas *cEffcfPos=new TCanvas("cEffcfPos","cEffcfPos");
328 cEffcfPos->SetGridy();
329 for(Int_t ipart=0;ipart<=2;ipart++){
330 if(ipart==0)hEffPos[ipart]->Draw("P");
331 else hEffPos[ipart]->Draw("PSAME");
332 hEffPosMean[ipart]->Draw("PSAME");
334 cEffcfPos->BuildLegend();
336 TCanvas *cEffcfNeg=new TCanvas("cEffcfNeg","cEffcfNeg");
337 cEffcfNeg->SetGridy();
338 for(Int_t ipart=0;ipart<=2;ipart++){
339 if(ipart==0)hEffNeg[ipart]->Draw("P");
340 else hEffNeg[ipart]->Draw("PSAME");
341 hEffNegMean[ipart]->Draw("PSAME");
343 cEffcfNeg->BuildLegend();
349 for(Int_t ipart=0;ipart<=2;ipart++) {
350 hHistPosNSigma[ipart]->Scale(normMC,"width");
351 hHistNegNSigma[ipart]->Scale(normMC,"width");
352 hHistPosNSigmaMean[ipart]->Scale(normMC,"width");
353 hHistNegNSigmaMean[ipart]->Scale(normMC,"width");
357 //////////////////////////////Prim/All Sec/All SecStr/All
369 for(Int_t ipart=0;ipart<=2;ipart++){
370 hPrimPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCposReco%i",ipart));
371 hPrimNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistPrimMCnegReco%i",ipart));
372 hSecPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecMatMCposReco%i",ipart));
373 hSecNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecMatMCnegReco%i",ipart));
374 hSecStrPos[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecStrMCposReco%i",ipart));
375 hSecStrNeg[ipart]=(TH1F*)cOutputMC->FindObject(Form("fHistSecStrMCnegReco%i",ipart));
377 hAllSecPos[ipart]=(TH1F*)hSecPos[ipart]->Clone(Form("hAllSecPos%i",ipart));
378 hAllSecPos[ipart]->Add(hSecStrPos[ipart]);
379 hAllSecNeg[ipart]=(TH1F*)hSecNeg[ipart]->Clone(Form("hAllSecNeg%i",ipart));
380 hAllSecNeg[ipart]->Add(hSecStrNeg[ipart]);
381 hAllPos[ipart]=(TH1F*)hPrimPos[ipart]->Clone(Form("hAllPos%i",ipart));
382 hAllPos[ipart]->Add(hSecPos[ipart]);
383 hAllPos[ipart]->Add(hSecStrPos[ipart]);
384 hAllNeg[ipart]=(TH1F*)hPrimNeg[ipart]->Clone(Form("hAllNeg%i",ipart));
385 hAllNeg[ipart]->Add(hSecNeg[ipart]);
386 hAllNeg[ipart]->Add(hSecStrNeg[ipart]);
388 hPrimPos[ipart]->Divide(hAllPos[ipart]);
389 hSecPos[ipart]->Divide(hAllPos[ipart]);
390 hSecStrPos[ipart]->Divide(hAllPos[ipart]);
391 hPrimNeg[ipart]->Divide(hAllNeg[ipart]);
392 hSecNeg[ipart]->Divide(hAllNeg[ipart]);
393 hSecStrNeg[ipart]->Divide(hAllNeg[ipart]);
395 hAllSecPos[ipart]->Divide(hAllPos[ipart]);
396 hAllSecNeg[ipart]->Divide(hAllNeg[ipart]);
397 for(Int_t binnum=1;binnum<=hAllSecNeg[ipart]->GetNbinsX();binnum++){
398 hAllSecPos[ipart]->SetBinContent(binnum,1-hAllSecPos[ipart]->GetBinContent(binnum));
399 hAllSecNeg[ipart]->SetBinContent(binnum,1-hAllSecNeg[ipart]->GetBinContent(binnum));
403 /////////////////////////////////////////////
404 TFile *outMC2=new TFile(foutMC.Data(),"update");
405 TList *lMC=new TList();
407 lMC->SetName(Form("MC_Mult%ito%i",lowmult,upmult));
408 for(Int_t ipart=0;ipart<=2;ipart++){
409 lMC->Add(hCorrFactPos[ipart]);
410 lMC->Add(hCorrFactNeg[ipart]);
411 lMC->Add(hCorrFactPosMean[ipart]);
412 lMC->Add(hCorrFactNegMean[ipart]);
413 lMC->Add(hEffPos[ipart]);
414 lMC->Add(hEffNeg[ipart]);
415 lMC->Add(hEffPosMean[ipart]);
416 lMC->Add(hEffNegMean[ipart]);
417 lMC->Add(hHistPosNSigma[ipart]);
418 lMC->Add(hHistNegNSigma[ipart]);
419 lMC->Add(hHistPosNSigmaMean[ipart]);
420 lMC->Add(hHistNegNSigmaMean[ipart]);
421 lMC->Add(hPrimPos[ipart]);
422 lMC->Add(hPrimNeg[ipart]);
423 lMC->Add(hSecPos[ipart]);
424 lMC->Add(hSecNeg[ipart]);
425 lMC->Add(hSecStrPos[ipart]);
426 lMC->Add(hSecStrNeg[ipart]);
427 lMC->Add(hAllSecPos[ipart]);
428 lMC->Add(hAllSecNeg[ipart]);
430 lMC->Write(Form("MC_Mult%ito%i",lowmult,upmult),1);
435 TString lnameDATA=Form("clistITSsaMult%ito%i",lowmult,upmult);
436 Printf("\n\n------------------ READING DATA %s -------------------\n",lnameDATA.Data());
437 TFile *finDATA = new TFile(fnameDATA.Data());
438 TDirectory *dirFileDATA=(TDirectory*)finDATA->Get("PWG2SpectraITSsa");
439 TList *cOutputDATA = (TList*)dirFileDATA->Get(lnameDATA.Data());
442 for(Int_t ipart=0;ipart<=2;ipart++){
443 hHistPosNSigmaDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistPosNSigma%i",ipart));
444 hHistNegNSigmaDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistNegNSigma%i",ipart));
445 hHistPosNSigmaMeanDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistPosNSigmaMean%i",ipart));
446 hHistNegNSigmaMeanDATA[ipart]=(TH1F*)cOutputDATA->FindObject(Form("hHistNegNSigmaMean%i",ipart));
448 hITSsaPos[ipart]=(TH1F*)hHistPosNSigmaDATA[ipart]->Clone(Form("hITSsaPos%i",ipart));
449 hITSsaNeg[ipart]=(TH1F*)hHistNegNSigmaDATA[ipart]->Clone(Form("hITSsaNeg%i",ipart));
450 hITSsaPosMean[ipart]=(TH1F*)hHistPosNSigmaMeanDATA[ipart]->Clone(Form("hITSsaPosAsymm%i",ipart));
451 hITSsaNegMean[ipart]=(TH1F*)hHistNegNSigmaMeanDATA[ipart]->Clone(Form("hITSsaNegAsymm%i",ipart));
452 hITSsaRawPos[ipart]=(TH1F*)hHistPosNSigmaDATA[ipart]->Clone(Form("hITSsaRawPos%i",ipart));
453 hITSsaRawNeg[ipart]=(TH1F*)hHistNegNSigmaDATA[ipart]->Clone(Form("hITSsaRawNeg%i",ipart));
454 hITSsaRawPosMean[ipart]=(TH1F*)hHistPosNSigmaMeanDATA[ipart]->Clone(Form("hITSsaRawPosAsymm%i",ipart));
455 hITSsaRawNegMean[ipart]=(TH1F*)hHistNegNSigmaMeanDATA[ipart]->Clone(Form("hITSsaRawNegAsymm%i",ipart));
457 hITSsaPos[ipart]->Sumw2();
458 hITSsaNeg[ipart]->Sumw2();
459 hITSsaPosMean[ipart]->Sumw2();
460 hITSsaNegMean[ipart]->Sumw2();
461 hITSsaRawPos[ipart]->Sumw2();
462 hITSsaRawNeg[ipart]->Sumw2();
463 hITSsaRawPosMean[ipart]->Sumw2();
464 hITSsaRawNegMean[ipart]->Sumw2();
466 hITSsaPos[ipart]->Divide(hCorrFactPos[ipart]);
467 hITSsaNeg[ipart]->Divide(hCorrFactNeg[ipart]);
468 hITSsaPosMean[ipart]->Divide(hCorrFactPosMean[ipart]);
469 hITSsaNegMean[ipart]->Divide(hCorrFactNegMean[ipart]);
471 SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hITSsaPos[ipart]);
472 SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hITSsaNeg[ipart]);
473 SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hITSsaPosMean[ipart]);
474 SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hITSsaNegMean[ipart]);
475 SetDrawAtt(21+ipart,ipart+4,1,ipart+4,1,hITSsaRawPos[ipart]);
476 SetDrawAtt(25+ipart,ipart+4,1,ipart+4,1,hITSsaRawNeg[ipart]);
477 SetDrawAtt(21+ipart,ipart+1,1,ipart+1,1,hITSsaRawPosMean[ipart]);
478 SetDrawAtt(25+ipart,ipart+1,1,ipart+1,1,hITSsaRawNegMean[ipart]);
481 if(!CorrectOnlyForEfficiency){
485 printf("- Correction from DCA\n");
487 // Correction factor based on fit to DCA distr for P and Pibar
488 TString fnameDCAPosP=Form("./DCACorr%s_%s_%.0fDCA_PosP_TFraction.root",period.Data(),MCname.Data(),DCAcut);
489 TFile *fDCAPosP= new TFile(fnameDCAPosP.Data());
490 TH1F * hDCAPosP= (TH1F*)fDCAPosP->Get("hPrimAllDATAMC");
491 //TH1F * hDCAPosP= (TH1F*)fDCAPosP->Get("hPrimAllDATA");
492 TString fnameDCANegP=Form("./DCACorr%s_%s_%.0fDCA_NegP_TFraction.root",period.Data(),MCname.Data(),DCAcut);
493 TFile *fDCANegP= new TFile(fnameDCANegP.Data());
494 TH1F * hDCANegP= (TH1F*)fDCANegP->Get("hPrimAllDATAMC");
495 //TH1F * hDCANegP= (TH1F*)fDCANegP->Get("hPrimAllDATA");
497 for(Int_t pbin=0;pbin<=hDCAPosP->GetNbinsX();pbin++)hDCAPosP->SetBinError(pbin,0);
498 for(Int_t pbin=0;pbin<=hDCANegP->GetNbinsX();pbin++)hDCANegP->SetBinError(pbin,0);
500 hITSsaPos[2]->Multiply(hDCAPosP);
501 hITSsaNeg[2]->Multiply(hDCANegP);
502 hITSsaPosMean[2]->Multiply(hDCAPosP);
503 hITSsaNegMean[2]->Multiply(hDCANegP);
506 hDCANegP->Draw("same");
510 Printf("NO DCA correction for protons");
514 printf("- Geant/FLUKA correction\n");
515 //GeantFluka correction for Pi and K
516 TString fnameGeanFlukaPi="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSection.211.root";
517 TFile *fGeanFlukaPi= new TFile(fnameGeanFlukaPi.Data());
518 TH1F *hGeantFlukaPiPos=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionParticles");
519 TH1F *hGeantFlukaPiNeg=(TH1F*)fGeanFlukaPi->Get("gHistCorrectionForCrossSectionAntiParticles");
520 for(Int_t binPi=0;binPi<=hITSsaPos[0]->GetNbinsX();binPi++){
521 Float_t FlukaCorrPiPos=hGeantFlukaPiPos->GetBinContent(hGeantFlukaPiPos->FindBin(hITSsaPos[0]->GetBinCenter(binPi)));
522 Float_t FlukaCorrPiNeg=hGeantFlukaPiNeg->GetBinContent(hGeantFlukaPiNeg->FindBin(hITSsaNeg[0]->GetBinCenter(binPi)));
523 hITSsaPos[0]->SetBinContent(binPi,hITSsaPos[0]->GetBinContent(binPi)*FlukaCorrPiPos);
524 hITSsaPos[0]->SetBinError(binPi,hITSsaPos[0]->GetBinError(binPi)*FlukaCorrPiPos);
525 hITSsaNeg[0]->SetBinContent(binPi,hITSsaNeg[0]->GetBinContent(binPi)*FlukaCorrPiNeg);
526 hITSsaNeg[0]->SetBinError(binPi,hITSsaNeg[0]->GetBinError(binPi)*FlukaCorrPiNeg);
527 hITSsaPosMean[0]->SetBinContent(binPi,hITSsaPosMean[0]->GetBinContent(binPi)*FlukaCorrPiPos);
528 hITSsaPosMean[0]->SetBinError(binPi,hITSsaPosMean[0]->GetBinError(binPi)*FlukaCorrPiPos);
529 hITSsaNegMean[0]->SetBinContent(binPi,hITSsaNegMean[0]->GetBinContent(binPi)*FlukaCorrPiNeg);
530 hITSsaNegMean[0]->SetBinError(binPi,hITSsaNegMean[0]->GetBinError(binPi)*FlukaCorrPiNeg);
533 TString fnameGeanFlukaK="$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSection.321.root";
534 TFile *fGeanFlukaK= new TFile(fnameGeanFlukaK.Data());
535 TH1F *hGeantFlukaKPos=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionParticles");
536 TH1F *hGeantFlukaKNeg=(TH1F*)fGeanFlukaK->Get("gHistCorrectionForCrossSectionAntiParticles");
537 for(Int_t binK=0;binK<=hITSsaPos[1]->GetNbinsX();binK++){
538 Float_t FlukaCorrKPos=hGeantFlukaKPos->GetBinContent(hGeantFlukaKPos->FindBin(hITSsaPos[1]->GetBinCenter(binK)));
539 Float_t FlukaCorrKNeg=hGeantFlukaKNeg->GetBinContent(hGeantFlukaKNeg->FindBin(hITSsaNeg[1]->GetBinCenter(binK)));
540 hITSsaPos[1]->SetBinContent(binK,hITSsaPos[1]->GetBinContent(binK)*FlukaCorrKPos);
541 hITSsaNeg[1]->SetBinContent(binK,hITSsaNeg[1]->GetBinContent(binK)*FlukaCorrKNeg);
542 hITSsaPos[1]->SetBinError(binK,hITSsaPos[1]->GetBinError(binK)*FlukaCorrKPos);
543 hITSsaNeg[1]->SetBinError(binK,hITSsaNeg[1]->GetBinError(binK)*FlukaCorrKNeg);
544 hITSsaPosMean[1]->SetBinContent(binK,hITSsaPosMean[1]->GetBinContent(binK)*FlukaCorrKPos);
545 hITSsaNegMean[1]->SetBinContent(binK,hITSsaNegMean[1]->GetBinContent(binK)*FlukaCorrKNeg);
546 hITSsaPosMean[1]->SetBinError(binK,hITSsaPosMean[1]->GetBinError(binK)*FlukaCorrKPos);
547 hITSsaNegMean[1]->SetBinError(binK,hITSsaNegMean[1]->GetBinError(binK)*FlukaCorrKNeg);
551 // ITS specific file for protons/antiprotons
555 TFile* fITS = new TFile ("$ALICE_ROOT/PWGLF/SPECTRA/PiKaPr/ITSsa/RootFilesGeantFlukaCorrection/correctionForCrossSectionITS_20100719.root");
556 // TH2D * hCorrFlukaITS[kNCharge];
557 TH2D * hCorrFlukaITS[2];
558 hCorrFlukaITS[kPos] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionProtons");
559 hCorrFlukaITS[kNeg] = (TH2D*)fITS->Get("gHistCorrectionForCrossSectionAntiProtons");
560 for(Int_t icharge = 0; icharge < kNCharge; icharge++){
561 Int_t nbins = hITSsaPos[2]->GetNbinsX();
562 Int_t nbinsy=hCorrFlukaITS[icharge]->GetNbinsY();
563 //for(Int_t binP=0;binP<=hITSsaPos[2]->GetNbinsX();binP++){
564 for(Int_t ibin = 0; ibin < nbins; ibin++){
565 Float_t pt = hITSsaPos[2]->GetBinCenter(ibin);
566 Float_t minPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(1);
567 Float_t maxPtCorrection = hCorrFlukaITS[icharge]->GetYaxis()->GetBinLowEdge(nbinsy+1);
568 if (pt < minPtCorrection) pt = minPtCorrection+0.0001;
569 if (pt > maxPtCorrection) pt = maxPtCorrection;
570 Float_t correction = hCorrFlukaITS[icharge]->GetBinContent(1,hCorrFlukaITS[icharge]->GetYaxis()->FindBin(pt));
571 cout<<correction<<" charge "<<icharge<<endl;
573 if (correction != 0) {// If the bin is empty this is a 0
574 hITSsaPos[2]->SetBinContent(ibin,hITSsaPos[2]->GetBinContent(ibin)*correction);
575 hITSsaPos[2]->SetBinError(ibin,hITSsaPos[2]->GetBinError (ibin)*correction);
576 hITSsaPosMean[2]->SetBinContent(ibin,hITSsaPosMean[2]->GetBinContent(ibin)*correction);
577 hITSsaPosMean[2]->SetBinError(ibin,hITSsaPosMean[2]->GetBinError (ibin)*correction);
579 }else if (hITSsaPos[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
580 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for protons secondaries, ITS, " << endl;
581 cout << " Bin content: " << hITSsaPos[2]->GetBinContent(ibin) << endl;
585 if (correction != 0) {// If the bin is empty this is a 0
586 hITSsaNeg[2]->SetBinContent(ibin,hITSsaNeg[2]->GetBinContent(ibin)*correction);
587 hITSsaNeg[2]->SetBinError(ibin,hITSsaNeg[2]->GetBinError (ibin)*correction);
588 hITSsaNegMean[2]->SetBinContent(ibin,hITSsaNegMean[2]->GetBinContent(ibin)*correction);
589 hITSsaNegMean[2]->SetBinError(ibin,hITSsaNegMean[2]->GetBinError (ibin)*correction);
590 }else if (hITSsaNeg[2]->GetBinContent(ibin) > 0) { // If we are skipping a non-empty bin, we notify the user
591 cout << "Fluka/GEANT: Not correcting bin "<<ibin << " for Antiprotons secondaries, ITS, " << endl;
592 cout << " Bin content: " << hITSsaNeg[2]->GetBinContent(ibin) << endl;
599 Printf("SPECTRA ARE NOT CORRECTED FOR GEANT/FLUKA");
604 Printf("SPECTRA ARE NOT CORRECTED FOR DCA FIT + GEANT/FLUKA");
608 TH1F *fHistNEvents =0x0;
610 Printf("Normalization to the number of events with good vertex");
611 TFile *finMAREK = new TFile(fnameDATA.Data());
612 TDirectory *dirFileMAREK=(TDirectory*)finMAREK->Get("PWG2SpectraITSTPC");
613 if(!dirFileMAREK)Printf("!dirFileMAREK");
614 TString lnameMAREK=Form("outputlow%dup%dHI0",lowmult,upmult);
615 TList *cOutputMAREK = (TList*)dirFileMAREK->Get(lnameMAREK.Data());
616 if(!cOutputMAREK)Printf("!cOutputMAREK");
617 TH1F *StatsHistMAREK=(TH1F*)cOutputMAREK->FindObject("StatsHist");
618 norm=(Float_t)1/StatsHistMAREK->GetBinContent(3);
621 Printf("Normalization to the number of events after phys sel");
622 fHistNEvents = (TH1F*)cOutputDATA->FindObject("fHistNEvents");
623 norm=(Float_t)1/fHistNEvents->GetBinContent(fHistNEvents->FindBin(1.));
625 if(NormalizeData)printf("- Normalization (NEvents): %.0f \n",1/norm);
626 else printf("- Normalization (NEvents): %.0f !!!!!!!SPECTRA ARE NOT NORMALIZED!!!!!!\n",fHistNEvents->GetBinContent(fHistNEvents->FindBin(1.)));
627 for(Int_t ipart=0;ipart<=2;ipart++){
629 hITSsaPos[ipart]->Scale(norm,"width");
630 hITSsaNeg[ipart]->Scale(norm,"width");
631 hITSsaPosMean[ipart]->Scale(norm,"width");
632 hITSsaNegMean[ipart]->Scale(norm,"width");
633 hITSsaRawPos[ipart]->Scale(norm,"width");
634 hITSsaRawNeg[ipart]->Scale(norm,"width");
635 hITSsaRawPosMean[ipart]->Scale(norm,"width");
636 hITSsaRawNegMean[ipart]->Scale(norm,"width");
639 hITSsaPos[ipart]->Scale(1,"width");
640 hITSsaNeg[ipart]->Scale(1,"width");
641 hITSsaPosMean[ipart]->Scale(1,"width");
642 hITSsaNegMean[ipart]->Scale(1,"width");
643 hITSsaRawPos[ipart]->Scale(1,"width");
644 hITSsaRawNeg[ipart]->Scale(1,"width");
645 hITSsaRawPosMean[ipart]->Scale(1,"width");
646 hITSsaRawNegMean[ipart]->Scale(1,"width");
650 printf("Normalization to INEL : %.2f \n",INEL);
651 hITSsaPos[ipart]->Scale(INEL,"");
652 hITSsaNeg[ipart]->Scale(INEL,"");
653 hITSsaPosMean[ipart]->Scale(INEL,"");
654 hITSsaNegMean[ipart]->Scale(INEL,"");
655 hITSsaRawPos[ipart]->Scale(INEL,"");
656 hITSsaRawNeg[ipart]->Scale(INEL,"");
657 hITSsaRawPosMean[ipart]->Scale(INEL,"");
658 hITSsaRawNegMean[ipart]->Scale(INEL,"");
660 else printf("Normalization to INEL ignored \n");
662 hITSsaPos[ipart]->SetMinimum(0.00001);
663 hITSsaNeg[ipart]->SetMinimum(0.00001);
664 hITSsaPosMean[ipart]->SetMinimum(0.00001);
665 hITSsaNegMean[ipart]->SetMinimum(0.00001);
666 hITSsaRawPos[ipart]->SetMinimum(0.00001);
667 hITSsaRawNeg[ipart]->SetMinimum(0.00001);
668 hITSsaRawPosMean[ipart]->SetMinimum(0.00001);
669 hITSsaRawNegMean[ipart]->SetMinimum(0.00001);
670 hITSsaPos[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
671 hITSsaNeg[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
672 hITSsaPosMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
673 hITSsaNegMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
674 hITSsaRawPos[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
675 hITSsaRawNeg[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
676 hITSsaRawPosMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
677 hITSsaRawNegMean[ipart]->GetXaxis()->SetTitle("Pt [GeV/c]");
678 hITSsaPos[ipart]->SetTitle(Form("ITSsa NSigma Symm Pos%i",ipart));
679 hITSsaNeg[ipart]->SetTitle(Form("ITSsa NSigma Symm Neg%i",ipart));
680 hITSsaPosMean[ipart]->SetTitle(Form("ITSsa NSigma Asymm Pos%i",ipart));
681 hITSsaNegMean[ipart]->SetTitle(Form("ITSsa NSigma Asymm Neg%i",ipart));
682 hITSsaRawPos[ipart]->SetTitle(Form("ITSsaRaw NSigma Symm Pos%i",ipart));
683 hITSsaRawNeg[ipart]->SetTitle(Form("ITSsaRaw NSigma Symm Neg%i",ipart));
684 hITSsaRawPosMean[ipart]->SetTitle(Form("ITSsaRaw NSigma Asymm Pos%i",ipart));
685 hITSsaRawNegMean[ipart]->SetTitle(Form("ITSsaRaw NSigma Asymm Neg%i",ipart));
688 ResetOutOfRange(hITSsaPos[ipart],ipart,lowRange,upRange);
689 ResetOutOfRange(hITSsaNeg[ipart],ipart,lowRange,upRange);
690 ResetOutOfRange(hITSsaPosMean[ipart],ipart,lowRange,upRange);
691 ResetOutOfRange(hITSsaNegMean[ipart],ipart,lowRange,upRange);
694 TH1F *hSystematicsITSsaPos[3];
695 TH1F *hSystematicsITSsaNeg[3];
696 TH1F *hSystematicsITSsaPosMean[3];
697 TH1F *hSystematicsITSsaNegMean[3];
700 printf("- making Systematic Error\n");
701 //preliminary systematic error
702 TFile *fsys = new TFile("RootFilesSystError/ITSsa-systematics-20101014.root");
703 //fsys->Print("all");
704 //cout<<fsys->GetSize()<<endl;
707 hSystPos[0]=(TH1F*)fsys->Get("hSystTotPosPion");
708 hSystNeg[0]=(TH1F*)fsys->Get("hSystTotNegPion");
709 hSystPos[1]=(TH1F*)fsys->Get("hSystTotPosKaon");
710 hSystNeg[1]=(TH1F*)fsys->Get("hSystTotNegKaon");
711 hSystPos[2]=(TH1F*)fsys->Get("hSystTotPosProton");
712 hSystNeg[2]=(TH1F*)fsys->Get("hSystTotNegProton");
714 for(Int_t ipart=0;ipart<=2;ipart++){
715 hSystematicsITSsaPos[ipart]=(TH1F*)hITSsaPos[ipart]->Clone(Form("hSystematicsITSsaPos%i",ipart));
716 hSystematicsITSsaNeg[ipart]=(TH1F*)hITSsaNeg[ipart]->Clone(Form("hSystematicsITSsaNeg%i",ipart));
717 hSystematicsITSsaPosMean[ipart]=(TH1F*)hITSsaPosMean[ipart]->Clone(Form("hSystematicsITSsaPosAsymm%i",ipart));
718 hSystematicsITSsaNegMean[ipart]=(TH1F*)hITSsaNegMean[ipart]->Clone(Form("hSystematicsITSsaNegAsymm%i",ipart));
720 for(Int_t ibin=0;ibin<=hSystematicsITSsaPos[ipart]->GetNbinsX();ibin++){
721 Float_t syserr=hSystPos[ipart]->GetBinContent(hSystPos[ipart]->FindBin(hSystematicsITSsaPos[ipart]->GetBinCenter(ibin)));
722 hSystematicsITSsaPos[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaPos[ipart]->GetBinContent(ibin));
723 hSystematicsITSsaPosMean[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaPosMean[ipart]->GetBinContent(ibin));
725 for(Int_t ibin=0;ibin<=hSystematicsITSsaNeg[ipart]->GetNbinsX();ibin++){
726 Float_t syserr=hSystNeg[ipart]->GetBinContent(hSystNeg[ipart]->FindBin(hSystematicsITSsaNeg[ipart]->GetBinCenter(ibin)));
727 hSystematicsITSsaNeg[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaNeg[ipart]->GetBinContent(ibin));
728 hSystematicsITSsaNegMean[ipart]->SetBinError(ibin,syserr*hSystematicsITSsaNegMean[ipart]->GetBinContent(ibin));
733 TCanvas *cSpectraSymm=new TCanvas("cSpectraSymm","cSpectraSymm");
734 cSpectraSymm->SetGridy();
735 for(Int_t ipart=0;ipart<=2;ipart++){
736 if(ipart==0)hITSsaPos[ipart]->DrawCopy("P");
737 else hITSsaPos[ipart]->DrawCopy("PSAME");
738 hITSsaNeg[ipart]->DrawCopy("PSAME");
740 cSpectraSymm->BuildLegend();
742 TCanvas *cSpectraAsymm=new TCanvas("cSpectraAsymm","cSpectraAsymm");
743 cSpectraAsymm->SetGridy();
744 for(Int_t ipart=0;ipart<=2;ipart++){
745 if(ipart==0)hITSsaPosMean[ipart]->DrawCopy("P");
746 else hITSsaPosMean[ipart]->DrawCopy("PSAME");
747 hITSsaNegMean[ipart]->DrawCopy("PSAME");
749 cSpectraAsymm->BuildLegend();
753 /////////////////////////// SCALING SPECTRA FOR TRACKING EFFICIENCY !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
754 Float_t ScalingFactor=1.015;
755 Printf("SCALING FACTOR: %f",ScalingFactor);
756 for(Int_t ipart=0;ipart<=2;ipart++){
757 hITSsaPosMean[ipart]->Scale(ScalingFactor);
758 Printf("Scaling %s by %f",hITSsaPosMean[ipart]->GetName(),ScalingFactor);
759 hITSsaNegMean[ipart]->Scale(ScalingFactor);
760 Printf("Scaling %s by %f",hITSsaNegMean[ipart]->GetName(),ScalingFactor);
761 hITSsaPos[ipart]->Scale(ScalingFactor);
762 Printf("Scaling %s by %f",hITSsaPos[ipart]->GetName(),ScalingFactor);
763 hITSsaNeg[ipart]->Scale(ScalingFactor);
764 Printf("Scaling %s by %f",hITSsaNeg[ipart]->GetName(),ScalingFactor);
767 TFile *outDATA2=new TFile(foutDATA.Data(),"update");
768 TList *lDATA=new TList();
770 lDATA->SetName(Form("DATA_Mult%ito%i",lowmult,upmult));
772 for(Int_t ipart=0;ipart<=2;ipart++) lDATA->Add(hITSsaPosMean[ipart]);
773 for(Int_t ipart=0;ipart<=2;ipart++) lDATA->Add(hITSsaNegMean[ipart]);
775 for(Int_t ipart=0;ipart<=2;ipart++){
776 lDATA->Add(fHistPrimMCposBefEvSel[ipart]);
777 lDATA->Add(fHistPrimMCnegBefEvSel[ipart]);
778 lDATA->Add(hHistPosNSigmaPrimMean[ipart]);
779 lDATA->Add(hHistNegNSigmaPrimMean[ipart]);
780 lDATA->Add(hHistPosNSigmaPrimMCMean[ipart]);
781 lDATA->Add(hHistNegNSigmaPrimMCMean[ipart]);
782 lDATA->Add(hITSsaPos[ipart]);
783 lDATA->Add(hITSsaNeg[ipart]);
784 lDATA->Add(hITSsaPosMean[ipart]);
785 lDATA->Add(hITSsaNegMean[ipart]);
786 lDATA->Add(hITSsaRawPos[ipart]);
787 lDATA->Add(hITSsaRawNeg[ipart]);
788 lDATA->Add(hITSsaRawPosMean[ipart]);
789 lDATA->Add(hITSsaRawNegMean[ipart]);
792 for(Int_t ipart=0;ipart<=2;ipart++){
793 lDATA->Add(hSystematicsITSsaPos[ipart]);
794 lDATA->Add(hSystematicsITSsaNeg[ipart]);
795 lDATA->Add( hSystematicsITSsaPosMean[ipart]);
796 lDATA->Add( hSystematicsITSsaNegMean[ipart]);
800 lDATA->Write(Form("DATA_Mult%ito%i",lowmult,upmult),1);
805 RatioPlot(hITSsaNeg,hITSsaPos,"Neg","Pos","Symm NSigma",Low,Up);
806 RatioPlot(hITSsaNegMean,hITSsaPosMean,"Neg","Pos","Assymm NSigma",Low,Up);
807 RatioPlot(hITSsaPos,hITSsaPosMean,"Symm","Asymm","Positive",Low,Up);
808 RatioPlot(hITSsaNeg,hITSsaNegMean,"Symm","Asymm","Negative",Low,Up);
814 void ResetOutOfRange(TH1F *histo, Int_t ipart, Double_t lowRange[3], Double_t upRange[3]){
815 // set to -1 the bin contents out of the selcted pt range
816 for(Int_t ibin=0;ibin<=histo->GetNbinsX();ibin++){
817 if(histo->GetBinCenter(ibin)<lowRange[ipart] || histo->GetBinCenter(ibin)>upRange[ipart]) histo->SetBinContent(ibin,-1);
821 void RatioPlot(TH1F **num,TH1F **den,TString t1,TString t2,TString opt,Double_t Low[],Double_t Up[])
823 gROOT->SetStyle("Plain");
824 gStyle->SetOptTitle(0);
826 TString title=Form("%s/%s - %s",t1.Data(),t2.Data(),opt.Data());
827 TCanvas *c=new TCanvas(title.Data(),title.Data());
829 for(Int_t ipart=0;ipart<=2;ipart++){
830 num[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
831 num[ipart]->SetMarkerColor(2);
832 num[ipart]->SetLineColor(2);
833 num[ipart]->SetStats(kFALSE);
834 den[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
835 den[ipart]->SetMarkerColor(1);
836 den[ipart]->SetLineColor(1);
837 den[ipart]->SetStats(kFALSE);
841 num[ipart]->Draw("");
842 den[ipart]->Draw("same");
843 TPaveText *tpave=new TPaveText(0.1,0.2,0.5,0.29,"brNDC");
844 tpave->SetBorderSize(0);
845 tpave->SetFillStyle(0);
846 tpave->SetFillColor(0);
847 tpave->SetTextColor(2);
848 tpave->SetTextSize(.04);
849 TText *txt1=tpave->AddText(t1.Data());
850 txt1->SetTextFont(62);
851 txt1->SetTextColor(2);
852 TText *txt2=tpave->AddText(t2.Data());
853 txt2->SetTextFont(62);
854 txt2->SetTextColor(1);
856 Num[ipart]=(TH1F*)num[ipart]->Clone(title.Data());
857 Num[ipart]->SetTitle(title.Data());
858 Num[ipart]->Divide(den[ipart]);
859 Num[ipart]->GetXaxis()->SetRangeUser(Low[ipart],Up[ipart]-0.00001);
860 Num[ipart]->SetMarkerColor(4);
861 Num[ipart]->SetLineColor(4);
862 Num[ipart]->SetMaximum(1.2);
863 Num[ipart]->SetMinimum(.8);
866 Num[ipart]->SetStats(kFALSE);
868 TPaveText *tpave_1=new TPaveText(0.3,0.7,0.7,0.99,"brNDC");
869 tpave_1->SetBorderSize(0);
870 tpave_1->SetFillStyle(0);
871 tpave_1->SetFillColor(0);
872 tpave_1->SetTextColor(2);
873 tpave_1->SetTextSize(.04);
874 TText *txt1_1=tpave_1->AddText(title.Data());
875 txt1_1->SetTextFont(62);
877 // TLine line = TLine(Low[ipart],1,Up[ipart],1);
878 // line.SetLineColor(2);
879 // line.SetLineStyle(2);
880 // line.SetLineWidth(2);
881 // line.DrawClone("same");
888 void SetDrawAtt(Int_t markerstyle,Int_t markercolor,Int_t markersize,Int_t linecolor,Int_t linewidth,TH1 *h1)
890 h1->SetMarkerStyle(markerstyle);
891 h1->SetMarkerColor(markercolor);
892 h1->SetMarkerSize(markersize);
893 h1->SetLineColor(linecolor);
894 h1->SetLineWidth(linewidth);