Changes for PbPb (Zaida)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / macros / HFPtSpectrum.C
CommitLineData
65e55bbd 1//
2// Macro to use the AliHFPtSpectrum class
3// to compute the feed-down corrections for heavy-flavor
4//
5// Z.Conesa, September 2010 (zconesa@in2p3.fr)
6//
7
8#include <Riostream.h>
9#include "TH1D.h"
10#include "TH1.h"
11#include "TH2F.h"
12#include "TFile.h"
13#include "TGraphAsymmErrors.h"
14#include "TCanvas.h"
86bdcd8c 15#include "TROOT.h"
65e55bbd 16
5541b811 17#include "AliHFSystErr.h"
65e55bbd 18#include "AliHFPtSpectrum.h"
19
20//
21// Macro execution parameters:
22// 0) filename with the theoretical predictions (direct & feed-down)
23// 1) acceptance and reconstruction efficiencies file name (direct & feed-down)
24// 2) reconstructed spectra file name
25// 3) output file name
26// 4) Set the feed-down calculation option flag: 0=none, 1=fc only, 2=Nb only
27// 5) Set the luminosity
28// 6) Set the trigger efficiency
bb427707 29// 7-14) If the efficiency histos do not have the right bin width, set the files & histo-names, they'll be computed, if the efficiencies are in file (6), don't set this parameters
65e55bbd 30//
bb427707 31void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
32 const char *efffilename="Efficiencies.root",
5f3c1b97 33 const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0",
bb427707 34 const char *outfilename="HFPtSpectrum.root",
86bdcd8c 35 Int_t option=1, Double_t lumi=1.0, Double_t effTrig=1.0,
bb427707 36 const char *directsimufilename="", const char *directsimuhistoname="CFHFccontainer0_New_3Prong_SelStep0_proj-pt",
37 const char *directrecofilename="", const char *directrecohistoname="CFHFccontainer0_New_3Prong_SelStep8_proj-pt",
38 const char *feeddownsimufilename="", const char *feeddownsimuhistoname="CFHFccontainer0allD_New_3Prong_SelStep0_proj-pt",
39 const char *feeddownrecofilename="", const char *feeddownrecohistoname="CFHFccontainer0allD_New_3Prong_SelStep8_proj-pt") {
65e55bbd 40
41 // Set if calculation considers asymmetric uncertainties or not
42 bool asym = true;
43
44 // Set the meson and decay
86bdcd8c 45 // (only D0 -> K pi, D+--> K pi pi & D* --> D0 pi implemented here)
65e55bbd 46 bool isD0Kpi = true;
47 bool isDplusKpipi = false;
86bdcd8c 48 bool isDstarD0pi = false;
49 if (isD0Kpi && isDplusKpipi && isDstarD0pi) {
50 cout << "Sorry, can not deal with more than one correction at the same time"<<endl;
65e55bbd 51 return;
52 }
53
54 if (option>2) {
55 cout<< "Bad calculation option, should be <=2"<<endl;
56 return;
57 }
58 if (option==0) asym = false;
59
60 //
61 // Get the histograms from the files
62 //
86bdcd8c 63 TH1D *hDirectMCpt=0; // Input MC c-->D spectra
64 TH1D *hFeedDownMCpt=0; // Input MC b-->D spectra
65 TH1D *hDirectMCptMax=0; // Input MC maximum c-->D spectra
66 TH1D *hDirectMCptMin=0; // Input MC minimum c-->D spectra
67 TH1D *hFeedDownMCptMax=0; // Input MC maximum b-->D spectra
68 TH1D *hFeedDownMCptMin=0; // Input MC minimum b-->D spectra
8998180c 69 TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra
86bdcd8c 70 TH1D *hDirectEffpt=0; // c-->D Acceptance and efficiency correction
71 TH1D *hFeedDownEffpt=0; // b-->D Acceptance and efficiency correction
72 TH1D *hRECpt=0; // all reconstructed D
bb427707 73 //
86bdcd8c 74 TH1D *hDirectSimulationpt=0; // Simulated c--D spectra (used to re-compute the efficiency)
75 TH1D *hDirectReconstructionpt=0; // Reconstructed c--D spectra (used to re-compute the efficiency)
76 TH1D *hFeedDownSimulationpt=0; // Simulated b--D spectra (used to re-compute the efficiency)
77 TH1D *hFeedDownReconstructionpt=0; // Reconstructed b--D spectra (used to re-compute the efficiency)
bb427707 78 //
65e55bbd 79
80 //
81 // Define/Get the input histograms
82 //
8998180c 83 Int_t decay=0;
65e55bbd 84 TFile * mcfile = new TFile(mcfilename,"read");
85 if (isD0Kpi){
8998180c 86 decay = 1;
65e55bbd 87 hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central");
88 hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central");
53f63c75 89 hDirectMCptMax = (TH1D*)mcfile->Get("hD0Kpipred_max");
90 hDirectMCptMin = (TH1D*)mcfile->Get("hD0Kpipred_min");
91 hFeedDownMCptMax = (TH1D*)mcfile->Get("hD0KpifromBpred_max");
92 hFeedDownMCptMin = (TH1D*)mcfile->Get("hD0KpifromBpred_min");
8998180c 93 gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction");
65e55bbd 94 }
95 else if (isDplusKpipi){
8998180c 96 decay = 2;
65e55bbd 97 hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central");
98 hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central");
53f63c75 99 hDirectMCptMax = (TH1D*)mcfile->Get("hDpluskpipipred_max");
100 hDirectMCptMin = (TH1D*)mcfile->Get("hDpluskpipipred_min");
101 hFeedDownMCptMax = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max");
102 hFeedDownMCptMin = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min");
8998180c 103 gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction");
65e55bbd 104 }
86bdcd8c 105 else if(isDstarD0pi){
8998180c 106 decay = 3;
86bdcd8c 107 hDirectMCpt = (TH1D*)mcfile->Get("hDstarD0pipred_central");
108 hFeedDownMCpt = (TH1D*)mcfile->Get("hDstarD0pifromBpred_central");
109 hDirectMCptMax = (TH1D*)mcfile->Get("hDstarD0pipred_max");
110 hDirectMCptMin = (TH1D*)mcfile->Get("hDstarD0pipred_min");
111 hFeedDownMCptMax = (TH1D*)mcfile->Get("hDstarD0pifromBpred_max");
112 hFeedDownMCptMin = (TH1D*)mcfile->Get("hDstarD0pifromBpred_min");
8998180c 113 gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction");
86bdcd8c 114 }
65e55bbd 115 //
116 hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra");
117 hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra");
53f63c75 118 hDirectMCptMax->SetNameTitle("hDirectMCptMax","max direct MC spectra");
119 hDirectMCptMin->SetNameTitle("hDirectMCptMin","min direct MC spectra");
120 hFeedDownMCptMax->SetNameTitle("hFeedDownMCptMax","max feed-down MC spectra");
121 hFeedDownMCptMin->SetNameTitle("hFeedDownMCptMin","min feed-down MC spectra");
65e55bbd 122 //
bb427707 123 //
124 if (strcmp(directsimufilename,"")!=0 && strcmp(directrecofilename,"")!=0 &&
125 strcmp(feeddownsimufilename,"")!=0 && strcmp(feeddownrecofilename,"")!=0 ) {
126 if (strcmp(directsimufilename,"")!=0){
127 TFile *directSimufile = new TFile(directsimufilename,"read");
128 hDirectSimulationpt = (TH1D*)directSimufile->Get(directsimuhistoname);
129 hDirectSimulationpt->SetNameTitle("hDirectSimulationpt","hDirectSimulationpt");
130 }
131 if (strcmp(directrecofilename,"")!=0){
132 TFile *directRecofile = new TFile(directrecofilename,"read");
133 hDirectReconstructionpt = (TH1D*)directRecofile->Get(directrecohistoname);
134 hDirectReconstructionpt->SetNameTitle("hDirectReconstructionpt","hDirectReconstructionpt");
135 }
136 if (strcmp(feeddownsimufilename,"")!=0){
137 TFile *feeddownSimufile = new TFile(feeddownsimufilename,"read");
138 hFeedDownSimulationpt = (TH1D*)feeddownSimufile->Get(feeddownsimuhistoname);
139 hFeedDownSimulationpt->SetNameTitle("hFeedDownSimulationpt","hFeedDownSimulationpt");
140 }
141 if (strcmp(feeddownrecofilename,"")!=0){
142 TFile *feeddownRecofile = new TFile(feeddownrecofilename,"read");
143 hFeedDownReconstructionpt = (TH1D*)feeddownRecofile->Get(feeddownrecohistoname);
144 hFeedDownReconstructionpt->SetNameTitle("hFeedDownReconstructionpt","hFeedDownReconstructionpt");
145 }
146 }
147 else {
148 TFile * efffile = new TFile(efffilename,"read");
86bdcd8c 149 hDirectEffpt = (TH1D*)efffile->Get("hEffD");
bb427707 150 hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
86bdcd8c 151 hFeedDownEffpt = (TH1D*)efffile->Get("hEffB");
bb427707 152 hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff");
153 }
154 //
65e55bbd 155 //
156 TFile * recofile = new TFile(recofilename,"read");
5f3c1b97 157 hRECpt = (TH1D*)recofile->Get(recohistoname);
65e55bbd 158 hRECpt->SetNameTitle("hRECpt","Reconstructed spectra");
159
160 //
161 // Define the output histograms
162 //
163 TFile *out = new TFile(outfilename,"recreate");
164 //
86bdcd8c 165 TH1D *histofc=0;
166 TH1D *histofcMax=0;
167 TH1D *histofcMin=0;
168 TH1D *histoYieldCorr=0;
169 TH1D *histoYieldCorrMax=0;
170 TH1D *histoYieldCorrMin=0;
171 TH1D *histoSigmaCorr=0;
172 TH1D *histoSigmaCorrMax=0;
173 TH1D *histoSigmaCorrMin=0;
65e55bbd 174 //
175 int nbins = hRECpt->GetNbinsX();
86bdcd8c 176 TGraphAsymmErrors * gYieldCorr = 0;
177 TGraphAsymmErrors * gSigmaCorr = 0;
178 TGraphAsymmErrors * gFcExtreme = 0;
179 TGraphAsymmErrors * gFcConservative = 0;
180 TGraphAsymmErrors * gYieldCorrExtreme = 0;
181 TGraphAsymmErrors * gSigmaCorrExtreme = 0;
182 TGraphAsymmErrors * gYieldCorrConservative = 0;
183 TGraphAsymmErrors * gSigmaCorrConservative = 0;
65e55bbd 184
185
186 //
187 // Main functionalities for the calculation
188 //
189
190 // Define and set the basic option flags
191 AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option);
192 spectra->SetFeedDownCalculationOption(option);
193 spectra->SetComputeAsymmetricUncertainties(asym);
194
195 // Feed the input histograms
bb427707 196 // reconstructed spectra
197 cout << " Setting the reconstructed spectrum,";
65e55bbd 198 spectra->SetReconstructedSpectrum(hRECpt);
bb427707 199 // acceptance and efficiency corrections
200 cout << " the files to compute the efficiency,";
201 if (strcmp(directsimufilename,"")!=0 && strcmp(directrecofilename,"")!=0 &&
202 strcmp(feeddownsimufilename,"")!=0 && strcmp(feeddownrecofilename,"")!=0 ) {
86bdcd8c 203 spectra->EstimateAndSetDirectEfficiencyRecoBin(hDirectSimulationpt,hDirectReconstructionpt);
204 spectra->EstimateAndSetFeedDownEfficiencyRecoBin(hFeedDownSimulationpt,hFeedDownReconstructionpt);
bb427707 205 }
206 else {
207 spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt);
208 }
209 // option specific histos (theory)
210 cout << " the theoretical spectra";
65e55bbd 211 if(option==1){
212 spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt);
53f63c75 213 if(asym) spectra->SetMCptDistributionsBounds(hDirectMCptMax,hDirectMCptMin,hFeedDownMCptMax,hFeedDownMCptMin);
65e55bbd 214 }
215 else if(option==2){
216 spectra->SetFeedDownMCptSpectra(hFeedDownMCpt);
53f63c75 217 if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCptMax,hFeedDownMCptMin);
65e55bbd 218 }
bb427707 219
220 cout << " and the normalization" <<endl;
65e55bbd 221 // Set normalization factors (uncertainties set to 0. as example)
8998180c 222 double lumiUnc = 0.10*lumi; // 10% uncertainty on the luminosity
223 spectra->SetLuminosity(lumi,lumiUnc);
53f63c75 224 spectra->SetTriggerEfficiency(effTrig,0.);
65e55bbd 225
8998180c 226 // Set the global uncertainties on the efficiencies (in percent)
227 double globalEffUnc = 0.15;
228 double globalBCEffRatioUnc = 0.15;
229 // double globalEffUnc = 0.;
230 // double globalBCEffRatioUnc = 0.;
231 spectra->SetAccEffPercentageUncertainty(globalEffUnc,globalBCEffRatioUnc);
232
65e55bbd 233 // Do the calculations
234 cout << " Doing the calculation... "<< endl;
53f63c75 235 double deltaY = 1.0;
236 double branchingRatioC = 1.0;
237 double branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction
238 spectra->ComputeHFPtSpectrum(deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
65e55bbd 239 cout << " ended the calculation, getting the histograms back " << endl;
240
8998180c 241 // Set the systematics externally
5541b811 242 AliHFSystErr *systematics = new AliHFSystErr();
243 systematics->Init(decay);
8998180c 244 bool combineFeedDown = true;
5541b811 245 spectra->ComputeSystUncertainties(systematics,combineFeedDown);
246
8998180c 247
65e55bbd 248 //
249 // Get the output histograms
250 //
251 // the corrected yield and cross-section
252 histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum();
253 histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum();
53f63c75 254 histoYieldCorrMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum();
255 histoYieldCorrMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum();
bb427707 256 histoSigmaCorrMax = (TH1D*)spectra->GetHistoUpperLimitCrossSectionFromYieldSpectrum();
257 histoSigmaCorrMin = (TH1D*)spectra->GetHistoLowerLimitCrossSectionFromYieldSpectrum();
86bdcd8c 258 histoYieldCorr->SetNameTitle("histoYieldCorr","corrected yield");
259 histoYieldCorrMax->SetNameTitle("histoYieldCorrMax","max corrected yield");
260 histoYieldCorrMin->SetNameTitle("histoYieldCorrMin","min corrected yield");
261 histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section");
262 histoSigmaCorrMax->SetNameTitle("histoSigmaCorrMax","max corrected invariant cross-section");
263 histoSigmaCorrMin->SetNameTitle("histoSigmaCorrMin","min corrected invariant cross-section");
264 // the efficiencies
265 if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection();
266 if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection();
267
268 // Get & Rename the TGraphs
65e55bbd 269 if (asym) {
270 gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
271 gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
86bdcd8c 272 gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme();
273 gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme();
274 gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative();
275 gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative();
65e55bbd 276 }
277
86bdcd8c 278 // Get & Rename the TGraphs
65e55bbd 279 if (option==0 && asym){
280 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)");
281 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)");
282 }
283 if (option==1){
284 // fc histos
53f63c75 285 histofc = (TH1D*)spectra->GetHistoFeedDownCorrectionFc();
286 histofcMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectionFc();
287 histofcMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectionFc();
86bdcd8c 288 histofc->SetNameTitle("histofc","fc correction factor");
289 histofcMax->SetNameTitle("histofcMax","max fc correction factor");
53f63c75 290 histofcMin->SetNameTitle("histofcMin","min fc correction factor");
65e55bbd 291 if (asym) {
65e55bbd 292 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)");
293 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)");
86bdcd8c 294 gFcExtreme = spectra->GetFeedDownCorrectionFcExtreme();
295 gFcExtreme->SetNameTitle("gFcExtreme","gFcExtreme");
296 gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by fc)");
297 gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by fc)");
298 gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
299 gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
300 gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by fc)");
301 gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by fc)");
65e55bbd 302 }
303 }
304 if (option==2 && asym) {
305 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)");
306 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)");
86bdcd8c 307 gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by Nb)");
308 gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by Nb)");
309 gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by Nb)");
310 gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by Nb)");
066998f0 311 gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
312 gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
65e55bbd 313 }
314
315 //
316 // Now, plot the results ! :)
317 //
318
319 cout << " Drawing the results ! " << endl;
320
5f3c1b97 321 // control plots
86bdcd8c 322 if (option==1) {
323
324 TCanvas *ceff = new TCanvas("ceff","efficiency drawing");
325 ceff->Divide(1,2);
326 ceff->cd(1);
327 hDirectEffpt->Draw();
328 ceff->cd(2);
329 hFeedDownEffpt->Draw();
330 ceff->Update();
331
332 TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
333 cTheoryRebin->Divide(1,2);
334 cTheoryRebin->cd(1);
335 hDirectMCpt->Draw("");
336 TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum();
337 hDirectMCptRebin->SetLineColor(2);
338 hDirectMCptRebin->Draw("same");
339 cTheoryRebin->cd(2);
340 hFeedDownMCpt->Draw("");
341 TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum();
342 hFeedDownRebin->SetLineColor(2);
343 hFeedDownRebin->Draw("same");
344 cTheoryRebin->Update();
345
346 TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
347 cTheoryRebinLimits->Divide(1,2);
348 cTheoryRebinLimits->cd(1);
349 hDirectMCptMax->Draw("");
350 TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum();
351 hDirectMCptMaxRebin->SetLineColor(2);
352 hDirectMCptMaxRebin->Draw("same");
353 hDirectMCptMin->Draw("same");
354 TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum();
355 hDirectMCptMinRebin->SetLineColor(2);
356 hDirectMCptMinRebin->Draw("same");
357 cTheoryRebinLimits->cd(2);
358 hFeedDownMCptMax->Draw("");
359 TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum();
360 hFeedDownMaxRebin->SetLineColor(2);
361 hFeedDownMaxRebin->Draw("same");
362 hFeedDownMCptMin->Draw("same");
363 TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum();
364 hFeedDownMinRebin->SetLineColor(2);
365 hFeedDownMinRebin->Draw("same");
366 cTheoryRebinLimits->Update();
367 }
5f3c1b97 368
65e55bbd 369 if (option==1) {
370
371 TCanvas * cfc = new TCanvas("cfc","Fc");
53f63c75 372 histofcMax->Draw("c");
65e55bbd 373 histofc->Draw("csame");
53f63c75 374 histofcMin->Draw("csame");
65e55bbd 375 cfc->Update();
376
377 if (asym) {
53f63c75 378 TH2F *histofcDraw= new TH2F("histofcDraw","histofc (for drawing)",100,0,33.25,100,0.01,1.25);
379 histofcDraw->SetStats(0);
380 histofcDraw->GetXaxis()->SetTitle("p_{T} [GeV]");
381 histofcDraw ->GetXaxis()->SetTitleSize(0.05);
382 histofcDraw->GetXaxis()->SetTitleOffset(0.95);
383 histofcDraw->GetYaxis()->SetTitle(" fc ");
384 histofcDraw->GetYaxis()->SetTitleSize(0.05);
86bdcd8c 385
386 if (gFcExtreme){
387
388// for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
389// Double_t center=0., value=0.;
390// gFcExtreme->GetPoint(item,center,value);
391// Double_t highunc = gFcExtreme->GetErrorYhigh(item) / value ;
392// Double_t lowunc = gFcExtreme->GetErrorYlow(item) / value ;
393// cout << "Fc extreme: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
394// }
395// for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
396// Double_t center=0., value=0.;
397// gFcConservative->GetPoint(item,center,value);
398// Double_t highunc = gFcConservative->GetErrorYhigh(item) / value ;
399// Double_t lowunc = gFcConservative->GetErrorYlow(item) / value ;
400// cout << "Fc conservative: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
401// }
402 TCanvas *cfcExtreme = new TCanvas("cfcExtreme","Extreme Asymmetric fc (TGraphAsymmErr)");
403 gFcExtreme->SetFillStyle(3006);
404 gFcExtreme->SetLineWidth(3);
405 gFcExtreme->SetMarkerStyle(20);
406 gFcExtreme->SetFillColor(2);
407 histofcDraw->Draw();
408 gFcExtreme->Draw("3same");
409
410 if(gFcConservative){
411 gFcConservative->SetFillStyle(3007);
412 gFcConservative->SetFillColor(4);
413 gFcConservative->Draw("3same");
414 }
415
416 cfcExtreme->Update();
417 }
65e55bbd 418 }
419
420 }
421
422 //
423 // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra)
424 //
425 TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma");
426 hDirectMCpt->SetMarkerStyle(20);
427 hDirectMCpt->SetMarkerColor(4);
428 hDirectMCpt->Draw("p");
429 histoSigmaCorr->SetMarkerStyle(21);
430 histoSigmaCorr->SetMarkerColor(2);
431 histoSigmaCorr->Draw("psame");
432 histoYieldCorr->SetMarkerStyle(22);
433 histoYieldCorr->SetMarkerColor(6);
434 histoYieldCorr->Draw("psame");
435 hRECpt->SetMarkerStyle(23);
436 hRECpt->SetMarkerColor(3);
437 hRECpt->Draw("psame");
438 cresult->SetLogy();
439 cresult->Update();
440
bb427707 441 TCanvas * cresult2 = new TCanvas("cresult2","corrected yield & sigma");
442 histoSigmaCorr->SetMarkerStyle(21);
443 histoSigmaCorr->SetMarkerColor(2);
444 histoSigmaCorr->Draw("p");
445 histoYieldCorr->SetMarkerStyle(22);
446 histoYieldCorr->SetMarkerColor(6);
447 histoYieldCorr->Draw("psame");
448 hRECpt->SetMarkerStyle(23);
449 hRECpt->SetMarkerColor(3);
450 hRECpt->Draw("psame");
451 cresult2->SetLogy();
452 cresult2->Update();
453
454
65e55bbd 455 if (asym) {
456
86bdcd8c 457 TH2F *histoDraw = new TH2F("histoDraw","histo (for drawing)",100,0,33.25,100,50.,1e7);
458 float max = 1.1*gYieldCorr->GetMaximum();
459 histoDraw->SetAxisRange(0.1,max,"Y");
460 histoDraw->SetStats(0);
461 histoDraw->GetXaxis()->SetTitle("p_{T} [GeV]");
462 histoDraw->GetXaxis()->SetTitleSize(0.05);
463 histoDraw->GetXaxis()->SetTitleOffset(0.95);
464 histoDraw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]");
465 histoDraw->GetYaxis()->SetTitleSize(0.05);
466 TCanvas * cyieldAsym = new TCanvas("cyieldAsym","Asymmetric corrected yield (TGraphAsymmErr)");
467 gYieldCorr->SetFillStyle(3001);
468 gYieldCorr->SetLineWidth(3);
469 gYieldCorr->SetMarkerStyle(20);
470 gYieldCorr->SetFillColor(3);
471 histoDraw->Draw();
472 gYieldCorr->Draw("3LPsame");
473 gYieldCorr->Draw("Xsame");
474 cyieldAsym->SetLogy();
475 cyieldAsym->Update();
476
477 TCanvas * cyieldExtreme = new TCanvas("cyieldExtreme","Extreme Asymmetric corrected yield (TGraphAsymmErr)");
478 histoYieldCorr->Draw();
479 gYieldCorrExtreme->SetFillStyle(3002);
480 gYieldCorrExtreme->SetLineWidth(3);
481 gYieldCorrExtreme->SetMarkerStyle(20);
482 gYieldCorrExtreme->SetFillColor(2);
483 histoYieldCorr->Draw();
484 gYieldCorr->Draw("3same");
485 gYieldCorrExtreme->Draw("3same");
486 cyieldExtreme->SetLogy();
487 cyieldExtreme->Update();
488
489 TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9);
490 max = 1.1*gSigmaCorr->GetMaximum();
491 histo2Draw->SetAxisRange(0.1,max,"Y");
492 histo2Draw->SetStats(0);
493 histo2Draw->GetXaxis()->SetTitle("p_{T} [GeV]");
494 histo2Draw->GetXaxis()->SetTitleSize(0.05);
495 histo2Draw->GetXaxis()->SetTitleOffset(0.95);
496 histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}");
497 histo2Draw->GetYaxis()->SetTitleSize(0.05);
498 TCanvas * csigmaAsym = new TCanvas("csigmaAsym","Asymmetric corrected sigma (TGraphAsymmErr)");
499 gSigmaCorr->SetFillStyle(3001);
500 gSigmaCorr->SetLineWidth(3);
501 gSigmaCorr->SetMarkerStyle(21);
502 gSigmaCorr->SetFillColor(3);
503 histo2Draw->Draw();
504 gSigmaCorr->Draw("3LPsame");
505 gSigmaCorr->Draw("Xsame");
506 csigmaAsym->SetLogy();
507 csigmaAsym->Update();
508
509// cout << endl <<" Sytematics (stat approach) " <<endl;
510// for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
511// Double_t center=0., value=0.;
512// gSigmaCorr->GetPoint(item,center,value);
513// Double_t highunc = gSigmaCorr->GetErrorYhigh(item) / value ;
514// Double_t lowunc = gSigmaCorr->GetErrorYlow(item) / value ;
515// cout << "Sigma syst (stat), i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
516// }
517
518 TCanvas * csigmaExtreme = new TCanvas("csigmaExtreme","Asymmetric extreme corrected sigma (TGraphAsymmErr)");
519 histoSigmaCorr->Draw();
520 gSigmaCorr->Draw("3Psame");
521 gSigmaCorrExtreme->SetFillStyle(3002);
522 gSigmaCorrExtreme->SetLineWidth(3);
523 gSigmaCorrExtreme->SetMarkerStyle(21);
524 gSigmaCorrExtreme->SetFillColor(2);
525 gSigmaCorrExtreme->Draw("3Psame");
526 csigmaExtreme->SetLogy();
527 csigmaExtreme->Update();
528
529// cout << endl << " Sytematics (Extreme approach)" <<endl;
530// for(Int_t item=0; item<gSigmaCorrExtreme->GetN(); item++){
531// Double_t center=0., value=0.;
532// gSigmaCorrExtreme->GetPoint(item,center,value);
533// Double_t highunc = gSigmaCorrExtreme->GetErrorYhigh(item) / value ;
534// Double_t lowunc = gSigmaCorrExtreme->GetErrorYlow(item) / value ;
535// cout << "Sigma syst (extreme) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
536// }
537
538// cout << endl << " Sytematics (Conservative approach)" <<endl;
539// for(Int_t item=0; item<gSigmaCorrConservative->GetN(); item++){
540// Double_t center=0., value=0.;
541// gSigmaCorrConservative->GetPoint(item,center,value);
542// Double_t highunc = gSigmaCorrConservative->GetErrorYhigh(item) / value ;
543// Double_t lowunc = gSigmaCorrConservative->GetErrorYlow(item) / value ;
544// cout << "Sigma syst (conservative) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
545// }
546
547 }
65e55bbd 548
549
550
551 //
552 // Write the histograms to the output file
553 //
5f3c1b97 554 cout << " Saving the results ! " << endl<< endl;
555
65e55bbd 556 out->cd();
557 //
86bdcd8c 558 hDirectMCpt->Write(); hFeedDownMCpt->Write();
53f63c75 559 hDirectMCptMax->Write(); hDirectMCptMin->Write();
560 hFeedDownMCptMax->Write(); hFeedDownMCptMin->Write();
86bdcd8c 561 if(hDirectEffpt) hDirectEffpt->Write(); if(hFeedDownEffpt) hFeedDownEffpt->Write();
65e55bbd 562 hRECpt->Write();
563 //
564 histoYieldCorr->Write();
53f63c75 565 histoYieldCorrMax->Write(); histoYieldCorrMin->Write();
65e55bbd 566 histoSigmaCorr->Write();
bb427707 567 histoSigmaCorrMax->Write(); histoSigmaCorrMin->Write();
65e55bbd 568
569 if(asym){
570 gYieldCorr->Write();
571 gSigmaCorr->Write();
86bdcd8c 572 if(gYieldCorrExtreme) gYieldCorrExtreme->Write();
573 if(gSigmaCorrExtreme) gSigmaCorrExtreme->Write();
574 if(gYieldCorrConservative) gYieldCorrConservative->Write();
575 if(gSigmaCorrConservative) gSigmaCorrConservative->Write();
066998f0 576 if(asym && gFcConservative) gFcConservative->Write();
65e55bbd 577 }
578
579 if(option==1){
580 histofc->Write();
86bdcd8c 581 histofcMax->Write(); histofcMin->Write();
582 if(asym && gFcExtreme) gFcExtreme->Write();
65e55bbd 583 }
584
8998180c 585
586 // Draw the cross-section
066998f0 587 // spectra->DrawSpectrum(gPrediction);
8998180c 588
5f3c1b97 589 // out->Close();
65e55bbd 590
591}