PWGHF/vertexingHF/macros/HFPtSpectrum.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <Riostream.h>
   3 #include "TH1D.h"
   4 #include "TH1.h"
   5 #include "TH2F.h"
   6 #include "TNtuple.h"
   7 #include "TFile.h"
   8 #include "TGraphAsymmErrors.h"
   9 #include "TCanvas.h"
  10 #include "TROOT.h"
  11 #include "TStyle.h"
  12 #include "TLegend.h"
  13 #include "AliHFSystErr.h"
  14 #include "AliHFPtSpectrum.h"
  15 #endif
  16
  17 /* $Id$ */
  18
  19 //
  20 // Macro to use the AliHFPtSpectrum class
  21 //  to compute the feed-down corrections for heavy-flavor
  22 //
  23 //  Z.Conesa, September 2010 (zconesa@in2p3.fr)
  24 //
  25
  26
  27
  28 //
  29 // Macro execution parameters:
  30 //  0) filename with the theoretical predictions  (direct & feed-down)
  31 //  1) acceptance and reconstruction efficiencies file name (direct & feed-down)
  32 //  2) reconstructed spectra file name
  33 //  3) output file name
  34 //  4) Set the feed-down calculation option flag: knone=none, kfc=fc only, kNb=Nb only
  35 //  5-6) Set the luminosity: the number of events analyzed, and the cross-section of the sample [nb]
  36 //  7) Set whether the yield is for particle + anti-particles or only one of the 'charges'
  37 //  8) Set the centrality class
  38 //  9) Flag to decide if there is need to evaluate the dependence on the energy loss
  39 //
  40
  41 enum centrality{ kpp7, kpp276, k010, k1020, k020, k2040, k4060, k6080, k4080, k80100 };
  42 enum BFDSubtrMethod { knone, kfc, kNb };
  43
  44 void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
  45                     const char *efffilename="Efficiencies.root",
  46                     const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0",
  47                     const char *outfilename="HFPtSpectrum.root",
  48                     Int_t fdMethod=kNb, Double_t nevents=1.0, Double_t sigma=1.0, // sigma[nb]
  49                     Bool_t isParticlePlusAntiParticleYield=true, Int_t cc=kpp7, Bool_t PbPbEloss=false ) {
  50
  51
  52   gROOT->Macro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
  53
  54   //  Set if calculation considers asymmetric uncertainties or not
  55   Bool_t asym = true;
  56
  57   // Set the meson and decay
  58   // (only D0 -> K pi, D+--> K pi pi & D* --> D0 pi & D+s -->KKpi implemented here)
  59   Bool_t isD0Kpi = true;
  60   Bool_t isDplusKpipi = false;
  61   Bool_t isDstarD0pi = false;
  62   Bool_t isDsKKpi = false;
  63   if (isD0Kpi && isDplusKpipi && isDstarD0pi && isDsKKpi) {
  64     cout << "Sorry, can not deal with more than one correction at the same time"<<endl;
  65     return;
  66   }
  67
  68   Int_t option=3;
  69   if (fdMethod==kfc) option=1;
  70   else if (fdMethod==kNb) option=2;
  71   else if (fdMethod==knone) option=0;
  72   else option=3;
  73
  74   if (option>2) {
  75     cout<< "Bad calculation option, should be <=2"<<endl;
  76     return;
  77   }
  78   if (option==0) asym = false;
  79
  80
  81   //
  82   // Defining the Tab values for the given centrality class
  83   // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
  84   //
  85   Double_t tab = 1., tabUnc = 0.;
  86   if ( cc == k010 ) {
  87     tab = 23.48; tabUnc = 0.97;
  88   } else if ( cc == k1020 ) {
  89     tab = 14.4318; tabUnc = 0.5733;
  90   } else if ( cc == k020 ) {
  91     tab = 18.93; tabUnc = 0.74;
  92   } else if ( cc == k2040 ) {
  93     tab = 6.86; tabUnc = 0.28;
  94   } else if ( cc == k4060 ) {
  95     tab = 2.00;  tabUnc= 0.11;
  96   } else if ( cc == k4080 ) {
  97     tab = 1.20451; tabUnc = 0.071843;
  98   } else if ( cc == k6080 ) {
  99     tab = 0.419; tabUnc = 0.033;
 100   } else if ( cc == k80100 ){
 101     tab = 0.0690; tabUnc = 0.0062;
 102   }
 103   tab *= 1e-9; // to pass from mb^{-1} to pb^{-1}
 104   tabUnc *= 1e-9;
 105
 106
 107
 108   //
 109   // Get the histograms from the files
 110   //
 111   TH1D *hDirectMCpt=0;           // Input MC c-->D spectra
 112   TH1D *hFeedDownMCpt=0;         // Input MC b-->D spectra
 113   TH1D *hDirectMCptMax=0;        // Input MC maximum c-->D spectra
 114   TH1D *hDirectMCptMin=0;        // Input MC minimum c-->D spectra
 115   TH1D *hFeedDownMCptMax=0;      // Input MC maximum b-->D spectra
 116   TH1D *hFeedDownMCptMin=0;      // Input MC minimum b-->D spectra
 117   //  TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra
 118   TH1D *hDirectEffpt=0;          // c-->D Acceptance and efficiency correction
 119   TH1D *hFeedDownEffpt=0;        // b-->D Acceptance and efficiency correction
 120   TH1D *hRECpt=0;                // all reconstructed D
 121
 122   //
 123   // Define/Get the input histograms
 124   //
 125   Int_t decay=0;
 126   TFile * mcfile = new TFile(mcfilename,"read");
 127   if (isD0Kpi){
 128     decay = 1;
 129     hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central");
 130     hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central_corr");
 131     hDirectMCptMax = (TH1D*)mcfile->Get("hD0Kpipred_max");
 132     hDirectMCptMin = (TH1D*)mcfile->Get("hD0Kpipred_min");
 133     hFeedDownMCptMax = (TH1D*)mcfile->Get("hD0KpifromBpred_max_corr");
 134     hFeedDownMCptMin = (TH1D*)mcfile->Get("hD0KpifromBpred_min_corr");
 135     //    gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction");
 136   }
 137   else if (isDplusKpipi){
 138     decay = 2;
 139     hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central");
 140     hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central_corr");
 141     hDirectMCptMax = (TH1D*)mcfile->Get("hDpluskpipipred_max");
 142     hDirectMCptMin = (TH1D*)mcfile->Get("hDpluskpipipred_min");
 143     hFeedDownMCptMax = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max_corr");
 144     hFeedDownMCptMin = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min_corr");
 145     //    gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction");
 146   }
 147   else if(isDstarD0pi){
 148     decay = 3;
 149     hDirectMCpt = (TH1D*)mcfile->Get("hDstarD0pipred_central");
 150     hFeedDownMCpt = (TH1D*)mcfile->Get("hDstarD0pifromBpred_central_corr");
 151     hDirectMCptMax = (TH1D*)mcfile->Get("hDstarD0pipred_max");
 152     hDirectMCptMin = (TH1D*)mcfile->Get("hDstarD0pipred_min");
 153     hFeedDownMCptMax = (TH1D*)mcfile->Get("hDstarD0pifromBpred_max_corr");
 154     hFeedDownMCptMin = (TH1D*)mcfile->Get("hDstarD0pifromBpred_min_corr");
 155     //    gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction");
 156   }
 157   else if (isDsKKpi){
 158     decay = 4;
 159     hDirectMCpt = (TH1D*)mcfile->Get("hDsKkpipred_central");
 160     hFeedDownMCpt = (TH1D*)mcfile->Get("hDsKkpifromBpred_central");
 161     hDirectMCptMax = (TH1D*)mcfile->Get("hDsKkpipred_max");
 162     hDirectMCptMin = (TH1D*)mcfile->Get("hDsKkpipred_min");
 163     hFeedDownMCptMax = (TH1D*)mcfile->Get("hDsKkpifromBpred_max");
 164     hFeedDownMCptMin = (TH1D*)mcfile->Get("hDsKkpifromBpred_min");
 165   }
 166   //
 167   hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra");
 168   hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra");
 169   hDirectMCptMax->SetNameTitle("hDirectMCptMax","max direct MC spectra");
 170   hDirectMCptMin->SetNameTitle("hDirectMCptMin","min direct MC spectra");
 171   hFeedDownMCptMax->SetNameTitle("hFeedDownMCptMax","max feed-down MC spectra");
 172   hFeedDownMCptMin->SetNameTitle("hFeedDownMCptMin","min feed-down MC spectra");
 173   //
 174   //
 175   TFile * efffile = new TFile(efffilename,"read");
 176   hDirectEffpt = (TH1D*)efffile->Get("hDirectEffpt");
 177   hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
 178   hFeedDownEffpt = (TH1D*)efffile->Get("hFeedDownEffpt");
 179   hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff");
 180   //
 181   //
 182   TFile * recofile = new TFile(recofilename,"read");
 183   hRECpt = (TH1D*)recofile->Get(recohistoname);
 184   hRECpt->SetNameTitle("hRECpt","Reconstructed spectra");
 185
 186   //
 187   // Define the output histograms
 188   //
 189   TFile *out = new TFile(outfilename,"recreate");
 190   //
 191   TH1D *histofc=0;
 192   TH1D *histofcMax=0;
 193   TH1D *histofcMin=0;
 194   TH1D *histoYieldCorr=0;
 195   TH1D *histoYieldCorrMax=0;
 196   TH1D *histoYieldCorrMin=0;
 197   TH1D *histoSigmaCorr=0;
 198   TH1D *histoSigmaCorrMax=0;
 199   TH1D *histoSigmaCorrMin=0;
 200   //
 201   TH2D *histofcRcb=0;
 202   TH1D *histofcRcb_px=0;
 203   TH2D *histoYieldCorrRcb=0;
 204   TH2D *histoSigmaCorrRcb=0;
 205   //
 206   TGraphAsymmErrors * gYieldCorr = 0;
 207   TGraphAsymmErrors * gSigmaCorr = 0;
 208   TGraphAsymmErrors * gFcExtreme = 0;
 209   TGraphAsymmErrors * gFcConservative = 0;
 210   TGraphAsymmErrors * gYieldCorrExtreme = 0;
 211   TGraphAsymmErrors * gSigmaCorrExtreme = 0;
 212   TGraphAsymmErrors * gYieldCorrConservative = 0;
 213   TGraphAsymmErrors * gSigmaCorrConservative = 0;
 214   //
 215   TNtuple * nSigma = 0;
 216
 217
 218   //
 219   // Main functionalities for the calculation
 220   //
 221
 222   // Define and set the basic option flags
 223   AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option);
 224   spectra->SetFeedDownCalculationOption(option);
 225   spectra->SetComputeAsymmetricUncertainties(asym);
 226   // Set flag on whether to additional PbPb Eloss hypothesis have to be computed
 227   spectra->SetComputeElossHypothesis(PbPbEloss);
 228
 229   // Feed the input histograms
 230   //  reconstructed spectra
 231   cout << " Setting the reconstructed spectrum,";
 232   spectra->SetReconstructedSpectrum(hRECpt);
 233   // acceptance and efficiency corrections
 234   cout << " the efficiency,";
 235   spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt);
 236   //    spectra->SetfIsStatUncEff(false);
 237   // option specific histos (theory)
 238   cout << " the theoretical spectra";
 239   if(option==1){
 240     spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt);
 241     if(asym) spectra->SetMCptDistributionsBounds(hDirectMCptMax,hDirectMCptMin,hFeedDownMCptMax,hFeedDownMCptMin);
 242   }
 243   else if(option==2){
 244     spectra->SetFeedDownMCptSpectra(hFeedDownMCpt);
 245     if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCptMax,hFeedDownMCptMin);
 246   }
 247
 248   cout << " and the normalization" <<endl;
 249   // Set normalization factors (uncertainties set to 0. as example)
 250   spectra->SetNormalization(nevents,sigma);
 251   Double_t lumi = nevents / sigma ;
 252   Double_t lumiUnc = 0.04*lumi; // 10% uncertainty on the luminosity
 253   spectra->SetLuminosity(lumi,lumiUnc);
 254   Double_t effTrig = 1.0;
 255   spectra->SetTriggerEfficiency(effTrig,0.);
 256
 257   // Set the global uncertainties on the efficiencies (in percent)
 258   Double_t globalEffUnc = 0.15;
 259   Double_t globalBCEffRatioUnc = 0.15;
 260   spectra->SetAccEffPercentageUncertainty(globalEffUnc,globalBCEffRatioUnc);
 261
 262   // Set if the yield is for particle+anti-particle or only one type
 263   spectra->SetIsParticlePlusAntiParticleYield(isParticlePlusAntiParticleYield);
 264
 265   // Set the Tab parameter and uncertainties
 266   if ( (cc != kpp7) && (cc != kpp276) ) {
 267     spectra->SetTabParameter(tab,tabUnc);
 268   }
 269
 270   // Do the calculations
 271   cout << " Doing the calculation... "<< endl;
 272   Double_t deltaY = 1.0;
 273   Double_t branchingRatioC = 1.0;
 274   Double_t branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction
 275   spectra->ComputeHFPtSpectrum(deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
 276   cout << "   ended the calculation, getting the histograms back " << endl;
 277
 278   // Set the systematics externally
 279   Bool_t combineFeedDown = true;
 280   AliHFSystErr *systematics = new AliHFSystErr();
 281   if( cc==kpp276 ) {
 282     systematics->SetIsLowEnergy(true);
 283   } else if( cc!=kpp7 )  {
 284     systematics->SetCollisionType(1);
 285     if ( cc == k010 )  systematics->SetCentrality("010");
 286     else if ( cc == k1020 )  systematics->SetCentrality("1020");
 287     else if ( cc == k020 )  systematics->SetCentrality("020");
 288     else if ( cc == k2040 ) {
 289       systematics->SetCentrality("2040");
 290       systematics->SetIsPbPb2010EnergyScan(true);
 291     }
 292     else if ( cc == k4060 )  systematics->SetCentrality("4060");
 293     else if ( cc == k6080 )  systematics->SetCentrality("6080");
 294     else if ( cc == k4080 ) systematics->SetCentrality("4080");
 295     else {
 296       cout << " Systematics not yet implemented " << endl;
 297       return;
 298     }
 299   } else { systematics->SetCollisionType(0); }
 300   systematics->Init(decay);
 301   spectra->ComputeSystUncertainties(systematics,combineFeedDown);
 302
 303   //
 304   // Get the output histograms
 305   //
 306   // the corrected yield and cross-section
 307   histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum();
 308   histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum();
 309   histoYieldCorrMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum();
 310   histoYieldCorrMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum();
 311   histoSigmaCorrMax = (TH1D*)spectra->GetHistoUpperLimitCrossSectionFromYieldSpectrum();
 312   histoSigmaCorrMin = (TH1D*)spectra->GetHistoLowerLimitCrossSectionFromYieldSpectrum();
 313   histoYieldCorr->SetNameTitle("histoYieldCorr","corrected yield");
 314   histoYieldCorrMax->SetNameTitle("histoYieldCorrMax","max corrected yield");
 315   histoYieldCorrMin->SetNameTitle("histoYieldCorrMin","min corrected yield");
 316   histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section");
 317   histoSigmaCorrMax->SetNameTitle("histoSigmaCorrMax","max corrected invariant cross-section");
 318   histoSigmaCorrMin->SetNameTitle("histoSigmaCorrMin","min corrected invariant cross-section");
 319   // the efficiencies
 320   if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection();
 321   if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection();
 322   // Get the PbPb Eloss hypothesis histograms
 323   if(PbPbEloss){
 324     histofcRcb = spectra->GetHistoFeedDownCorrectionFcVsEloss();
 325     histoYieldCorrRcb = spectra->GetHistoFeedDownCorrectedSpectrumVsEloss();
 326     histoSigmaCorrRcb = spectra->GetHistoCrossSectionFromYieldSpectrumVsEloss();
 327     histofcRcb->SetName("histofcRcb");
 328     histoYieldCorrRcb->SetName("histoYieldCorrRcb");
 329     histoSigmaCorrRcb->SetName("histoSigmaCorrRcb");
 330   }
 331
 332   // Get & Rename the TGraphs
 333   if (asym) {
 334     gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
 335     gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
 336     gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme();
 337     gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme();
 338     gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative();
 339     gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative();
 340   }
 341
 342   // Get & Rename the TGraphs
 343   if (option==0 && asym){
 344     gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)");
 345     gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)");
 346   }
 347   if (option==1){
 348     // fc histos
 349     histofc = (TH1D*)spectra->GetHistoFeedDownCorrectionFc();
 350     histofcMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectionFc();
 351     histofcMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectionFc();
 352     histofc->SetNameTitle("histofc","fc correction factor");
 353     histofcMax->SetNameTitle("histofcMax","max fc correction factor");
 354     histofcMin->SetNameTitle("histofcMin","min fc correction factor");
 355     if (asym) {
 356       gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)");
 357       gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)");
 358       gFcExtreme = spectra->GetFeedDownCorrectionFcExtreme();
 359       gFcExtreme->SetNameTitle("gFcExtreme","gFcExtreme");
 360       gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by fc)");
 361       gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by fc)");
 362       gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
 363       gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
 364       gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by fc)");
 365       gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by fc)");
 366     }
 367   }
 368   if (option==2 && asym) {
 369     gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)");
 370     gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)");
 371     gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by Nb)");
 372     gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by Nb)");
 373     gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by Nb)");
 374     gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by Nb)");
 375     gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
 376     gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
 377   }
 378
 379   if(PbPbEloss){
 380     nSigma = spectra->GetNtupleCrossSectionVsEloss();
 381   }
 382
 383   //
 384   // Now, plot the results ! :)
 385   //
 386
 387   gROOT->SetStyle("Plain");
 388
 389   cout << " Drawing the results ! " << endl;
 390
 391   // control plots
 392   if (option==1) {
 393
 394     TCanvas *ceff = new TCanvas("ceff","efficiency drawing");
 395     ceff->Divide(1,2);
 396     ceff->cd(1);
 397     hDirectEffpt->Draw();
 398     ceff->cd(2);
 399     hFeedDownEffpt->Draw();
 400     ceff->Update();
 401
 402     TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
 403     cTheoryRebin->Divide(1,2);
 404     cTheoryRebin->cd(1);
 405     hDirectMCpt->Draw("");
 406     TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum();
 407     hDirectMCptRebin->SetLineColor(2);
 408     hDirectMCptRebin->Draw("same");
 409     cTheoryRebin->cd(2);
 410     hFeedDownMCpt->Draw("");
 411     TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum();
 412     hFeedDownRebin->SetLineColor(2);
 413     hFeedDownRebin->Draw("same");
 414     cTheoryRebin->Update();
 415
 416     TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
 417     cTheoryRebinLimits->Divide(1,2);
 418     cTheoryRebinLimits->cd(1);
 419     hDirectMCptMax->Draw("");
 420     TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum();
 421     hDirectMCptMaxRebin->SetLineColor(2);
 422     hDirectMCptMaxRebin->Draw("same");
 423     hDirectMCptMin->Draw("same");
 424     TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum();
 425     hDirectMCptMinRebin->SetLineColor(2);
 426     hDirectMCptMinRebin->Draw("same");
 427     cTheoryRebinLimits->cd(2);
 428     hFeedDownMCptMax->Draw("");
 429     TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum();
 430     hFeedDownMaxRebin->SetLineColor(2);
 431     hFeedDownMaxRebin->Draw("same");
 432     hFeedDownMCptMin->Draw("same");
 433     TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum();
 434     hFeedDownMinRebin->SetLineColor(2);
 435     hFeedDownMinRebin->Draw("same");
 436     cTheoryRebinLimits->Update();
 437   }
 438
 439   if (option==1) {
 440
 441     TCanvas * cfc = new TCanvas("cfc","Fc");
 442     histofcMax->Draw("c");
 443     histofc->Draw("csame");
 444     histofcMin->Draw("csame");
 445     cfc->Update();
 446
 447     if (asym) {
 448       TH2F *histofcDraw= new TH2F("histofcDraw","histofc (for drawing)",100,0,33.25,100,0.01,1.25);
 449       histofcDraw->SetStats(0);
 450       histofcDraw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 451       histofcDraw ->GetXaxis()->SetTitleSize(0.05);
 452       histofcDraw->GetXaxis()->SetTitleOffset(0.95);
 453       histofcDraw->GetYaxis()->SetTitle(" fc ");
 454       histofcDraw->GetYaxis()->SetTitleSize(0.05);
 455
 456       if (gFcExtreme){
 457
 458 //      for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 459 //        Double_t center=0., value=0.;
 460 //        gFcExtreme->GetPoint(item,center,value);
 461 //        Double_t highunc = gFcExtreme->GetErrorYhigh(item) / value ;
 462 //        Double_t lowunc = gFcExtreme->GetErrorYlow(item) / value ;
 463 //        cout << "Fc extreme: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 464 //      }
 465 //      for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 466 //        Double_t center=0., value=0.;
 467 //        gFcConservative->GetPoint(item,center,value);
 468 //        Double_t highunc = gFcConservative->GetErrorYhigh(item) / value ;
 469 //        Double_t lowunc = gFcConservative->GetErrorYlow(item) / value ;
 470 //        cout << "Fc conservative: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 471 //      }
 472         TCanvas *cfcExtreme = new TCanvas("cfcExtreme","Extreme Asymmetric fc (TGraphAsymmErr)");
 473         gFcExtreme->SetFillStyle(3006);
 474         gFcExtreme->SetLineWidth(3);
 475         gFcExtreme->SetMarkerStyle(20);
 476         gFcExtreme->SetFillColor(2);
 477         histofcDraw->Draw();
 478         gFcExtreme->Draw("3same");
 479
 480         if(gFcConservative){
 481           gFcConservative->SetFillStyle(3007);
 482           gFcConservative->SetFillColor(4);
 483           gFcConservative->Draw("3same");
 484         }
 485
 486         cfcExtreme->Update();
 487       }
 488     }
 489
 490   }
 491
 492   //
 493   // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra)
 494   //
 495   TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma");
 496   hDirectMCpt->SetMarkerStyle(20);
 497   hDirectMCpt->SetMarkerColor(4);
 498   hDirectMCpt->Draw("p");
 499   histoSigmaCorr->SetMarkerStyle(21);
 500   histoSigmaCorr->SetMarkerColor(2);
 501   histoSigmaCorr->Draw("psame");
 502   histoYieldCorr->SetMarkerStyle(22);
 503   histoYieldCorr->SetMarkerColor(6);
 504   histoYieldCorr->Draw("psame");
 505   hRECpt->SetMarkerStyle(23);
 506   hRECpt->SetMarkerColor(3);
 507   hRECpt->Draw("psame");
 508   cresult->SetLogy();
 509   cresult->Update();
 510
 511   TCanvas * cresult2 = new TCanvas("cresult2","corrected yield & sigma");
 512   histoSigmaCorr->SetMarkerStyle(21);
 513   histoSigmaCorr->SetMarkerColor(2);
 514   histoSigmaCorr->Draw("p");
 515   histoYieldCorr->SetMarkerStyle(22);
 516   histoYieldCorr->SetMarkerColor(6);
 517   histoYieldCorr->Draw("psame");
 518   hRECpt->SetMarkerStyle(23);
 519   hRECpt->SetMarkerColor(3);
 520   hRECpt->Draw("psame");
 521   cresult2->SetLogy();
 522   cresult2->Update();
 523
 524
 525   if (asym) {
 526
 527     TH2F *histoDraw = new TH2F("histoDraw","histo (for drawing)",100,0,33.25,100,50.,1e7);
 528     float max = 1.1*gYieldCorr->GetMaximum();
 529     histoDraw->SetAxisRange(0.1,max,"Y");
 530     histoDraw->SetStats(0);
 531     histoDraw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 532     histoDraw->GetXaxis()->SetTitleSize(0.05);
 533     histoDraw->GetXaxis()->SetTitleOffset(0.95);
 534     histoDraw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]");
 535     histoDraw->GetYaxis()->SetTitleSize(0.05);
 536     TCanvas * cyieldAsym = new TCanvas("cyieldAsym","Asymmetric corrected yield (TGraphAsymmErr)");
 537     gYieldCorr->SetFillStyle(3001);
 538     gYieldCorr->SetLineWidth(3);
 539     gYieldCorr->SetMarkerStyle(20);
 540     gYieldCorr->SetFillColor(3);
 541     histoDraw->Draw();
 542     gYieldCorr->Draw("3LPsame");
 543     gYieldCorr->Draw("Xsame");
 544     cyieldAsym->SetLogy();
 545     cyieldAsym->Update();
 546
 547     TCanvas * cyieldExtreme = new TCanvas("cyieldExtreme","Extreme Asymmetric corrected yield (TGraphAsymmErr)");
 548     histoYieldCorr->Draw();
 549     gYieldCorrExtreme->SetFillStyle(3002);
 550     gYieldCorrExtreme->SetLineWidth(3);
 551     gYieldCorrExtreme->SetMarkerStyle(20);
 552     gYieldCorrExtreme->SetFillColor(2);
 553     histoYieldCorr->Draw();
 554     gYieldCorr->Draw("3same");
 555     gYieldCorrExtreme->Draw("3same");
 556     cyieldExtreme->SetLogy();
 557     cyieldExtreme->Update();
 558
 559     TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9);
 560     max = 1.1*gSigmaCorr->GetMaximum();
 561     histo2Draw->SetAxisRange(0.1,max,"Y");
 562     histo2Draw->SetStats(0);
 563     histo2Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 564     histo2Draw->GetXaxis()->SetTitleSize(0.05);
 565     histo2Draw->GetXaxis()->SetTitleOffset(0.95);
 566     histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}");
 567     histo2Draw->GetYaxis()->SetTitleSize(0.05);
 568     TCanvas * csigmaAsym = new TCanvas("csigmaAsym","Asymmetric corrected sigma (TGraphAsymmErr)");
 569     gSigmaCorr->SetFillStyle(3001);
 570     gSigmaCorr->SetLineWidth(3);
 571     gSigmaCorr->SetMarkerStyle(21);
 572     gSigmaCorr->SetFillColor(3);
 573     histo2Draw->Draw();
 574     gSigmaCorr->Draw("3LPsame");
 575     gSigmaCorr->Draw("Xsame");
 576     csigmaAsym->SetLogy();
 577     csigmaAsym->Update();
 578
 579 //     cout << endl <<" Sytematics (stat approach) " <<endl;
 580 //     for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 581 //       Double_t center=0., value=0.;
 582 //       gSigmaCorr->GetPoint(item,center,value);
 583 //       Double_t highunc = gSigmaCorr->GetErrorYhigh(item) / value ;
 584 //       Double_t lowunc = gSigmaCorr->GetErrorYlow(item) / value ;
 585 //       cout << "Sigma syst (stat), i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 586 //     }
 587
 588     TCanvas * csigmaExtreme = new TCanvas("csigmaExtreme","Asymmetric extreme corrected sigma (TGraphAsymmErr)");
 589     histoSigmaCorr->Draw();
 590     gSigmaCorr->Draw("3Psame");
 591     gSigmaCorrExtreme->SetFillStyle(3002);
 592     gSigmaCorrExtreme->SetLineWidth(3);
 593     gSigmaCorrExtreme->SetMarkerStyle(21);
 594     gSigmaCorrExtreme->SetFillColor(2);
 595     gSigmaCorrExtreme->Draw("3Psame");
 596     csigmaExtreme->SetLogy();
 597     csigmaExtreme->Update();
 598
 599 //     cout << endl << " Sytematics (Extreme approach)" <<endl;
 600 //     for(Int_t item=0; item<gSigmaCorrExtreme->GetN(); item++){
 601 //       Double_t center=0., value=0.;
 602 //       gSigmaCorrExtreme->GetPoint(item,center,value);
 603 //       Double_t highunc = gSigmaCorrExtreme->GetErrorYhigh(item) / value ;
 604 //       Double_t lowunc = gSigmaCorrExtreme->GetErrorYlow(item) / value ;
 605 //       cout << "Sigma syst (extreme) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 606 //     }
 607
 608 //     cout << endl << " Sytematics (Conservative approach)" <<endl;
 609 //     for(Int_t item=0; item<gSigmaCorrConservative->GetN(); item++){
 610 //       Double_t center=0., value=0.;
 611 //       gSigmaCorrConservative->GetPoint(item,center,value);
 612 //       Double_t highunc = gSigmaCorrConservative->GetErrorYhigh(item) / value ;
 613 //       Double_t lowunc = gSigmaCorrConservative->GetErrorYlow(item) / value ;
 614 //       cout << "Sigma syst (conservative) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 615 //     }
 616
 617  }
 618
 619   // Draw the PbPb Eloss hypothesis histograms
 620   if(PbPbEloss){
 621     AliHFPtSpectrum *CalcBins;
 622     gStyle->SetPalette(1);
 623     TCanvas *canvasfcRcb = new TCanvas("canvasfcRcb","fc vs pt vs Rcb");
 624     //    histofcRcb->Draw("cont4z");
 625     histofcRcb->Draw("colz");
 626     canvasfcRcb->Update();
 627     canvasfcRcb->cd(2);
 628     TCanvas *canvasfcRcb1 = new TCanvas("canvasfcRcb1","fc vs pt vs Rcb=1");
 629     histofcRcb_px = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px",40,40);
 630     histofcRcb_px->SetLineColor(2);
 631     if (option==1) {
 632       histofc->Draw();
 633       histofcRcb_px->Draw("same");
 634     } else histofcRcb_px->Draw("");
 635     canvasfcRcb1->Update();
 636     TCanvas *canvasfcRcb2 = new TCanvas("canvasfcRcb2","fc vs pt vs Rcb fixed Rcb");
 637     Int_t bin0 = CalcBins->FindTH2YBin(histofcRcb,0.25);
 638     Int_t bin1 = CalcBins->FindTH2YBin(histofcRcb,0.5);
 639     Int_t bin2 = CalcBins->FindTH2YBin(histofcRcb,1.0);
 640     Int_t bin3 = CalcBins->FindTH2YBin(histofcRcb,1.5);
 641     Int_t bin4 = CalcBins->FindTH2YBin(histofcRcb,2.0);
 642     Int_t bin5 = CalcBins->FindTH2YBin(histofcRcb,3.0);
 643     Int_t bin6 = CalcBins->FindTH2YBin(histofcRcb,4.0);
 644     TH1D * histofcRcb_px0a = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0a",bin0,bin0);
 645     TH1D * histofcRcb_px0 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0",bin1,bin1);
 646     TH1D * histofcRcb_px1 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px1",bin2,bin2);
 647     TH1D * histofcRcb_px2 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px2",bin3,bin3);
 648     TH1D * histofcRcb_px3 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px3",bin4,bin4);
 649     TH1D * histofcRcb_px4 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px4",bin5,bin5);
 650     TH1D * histofcRcb_px5 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px5",bin6,bin6);
 651     if (option==1) {
 652       histofc->Draw();
 653       //      histofcRcb_px->Draw("same");
 654     } else {
 655       //      histofcRcb_px->Draw("");
 656       histofcRcb_px0a->SetLineColor(2);
 657       histofcRcb_px0a->Draw("");
 658     }
 659     histofcRcb_px0a->SetLineColor(2);
 660     histofcRcb_px0a->Draw("same");
 661     histofcRcb_px0->SetLineColor(4);
 662     histofcRcb_px0->Draw("same");
 663     histofcRcb_px1->SetLineColor(3);
 664     histofcRcb_px1->Draw("same");
 665     histofcRcb_px2->SetLineColor(kCyan);
 666     histofcRcb_px2->Draw("same");
 667     histofcRcb_px3->SetLineColor(kMagenta+1);
 668     histofcRcb_px3->Draw("same");
 669     histofcRcb_px4->SetLineColor(kOrange+7);
 670     histofcRcb_px4->Draw("same");
 671     histofcRcb_px5->SetLineColor(kGreen+3);
 672     histofcRcb_px5->Draw("same");
 673     TLegend *legrcc = new TLegend(0.8,0.8,0.95,0.9);
 674     legrcc->SetFillColor(0);
 675     if (option==1) {
 676       legrcc->AddEntry(histofcRcb_px0a,"Rc/b=0.25","l");
 677       legrcc->AddEntry(histofcRcb_px0,"Rc/b=0.5","l");
 678       legrcc->AddEntry(histofcRcb_px1,"Rc/b=1.0","l");
 679       legrcc->AddEntry(histofcRcb_px2,"Rc/b=1.5","l");
 680       legrcc->AddEntry(histofcRcb_px3,"Rc/b=2.0","l");
 681       legrcc->AddEntry(histofcRcb_px4,"Rc/b=3.0","l");
 682       legrcc->AddEntry(histofcRcb_px5,"Rc/b=4.0","l");
 683     }else{
 684       legrcc->AddEntry(histofcRcb_px0a,"Rb=0.25","l");
 685       legrcc->AddEntry(histofcRcb_px0,"Rb=0.5","l");
 686       legrcc->AddEntry(histofcRcb_px1,"Rb=1.0","l");
 687       legrcc->AddEntry(histofcRcb_px2,"Rb=1.5","l");
 688       legrcc->AddEntry(histofcRcb_px3,"Rb=2.0","l");
 689       legrcc->AddEntry(histofcRcb_px4,"Rb=3.0","l");
 690       legrcc->AddEntry(histofcRcb_px5,"Rb=4.0","l");
 691     }
 692     legrcc->Draw();
 693     canvasfcRcb2->Update();
 694     TCanvas *canvasYRcb = new TCanvas("canvasYRcb","corrected yield vs pt vs Rcb");
 695     histoYieldCorrRcb->Draw("cont4z");
 696     canvasYRcb->Update();
 697     TCanvas *canvasSRcb = new TCanvas("canvasSRcb","sigma vs pt vs Rcb");
 698     histoSigmaCorrRcb->Draw("cont4z");
 699     canvasSRcb->Update();
 700     TCanvas *canvasSRcb1 = new TCanvas("canvasSRcb1","sigma vs pt vs Rcb fixed Rcb");
 701     TH1D * histoSigmaCorrRcb_px0a = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0a",bin0,bin0);
 702     TH1D * histoSigmaCorrRcb_px0 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0",bin1,bin1);
 703     TH1D * histoSigmaCorrRcb_px1 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px1",bin2,bin2);
 704     TH1D * histoSigmaCorrRcb_px2 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px2",bin3,bin3);
 705     TH1D * histoSigmaCorrRcb_px3 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px3",bin4,bin4);
 706     TH1D * histoSigmaCorrRcb_px4 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px4",bin5,bin5);
 707     TH1D * histoSigmaCorrRcb_px5 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px5",bin6,bin6);
 708     histoSigmaCorr->Draw();
 709     histoSigmaCorrRcb_px0a->SetLineColor(2);
 710     histoSigmaCorrRcb_px0a->Draw("hsame");
 711     histoSigmaCorrRcb_px0->SetLineColor(4);
 712     histoSigmaCorrRcb_px0->Draw("hsame");
 713     histoSigmaCorrRcb_px1->SetLineColor(3);
 714     histoSigmaCorrRcb_px1->Draw("hsame");
 715     histoSigmaCorrRcb_px2->SetLineColor(kCyan);
 716     histoSigmaCorrRcb_px2->Draw("hsame");
 717     histoSigmaCorrRcb_px3->SetLineColor(kMagenta+1);
 718     histoSigmaCorrRcb_px3->Draw("hsame");
 719     histoSigmaCorrRcb_px4->SetLineColor(kOrange+7);
 720     histoSigmaCorrRcb_px4->Draw("same");
 721     histoSigmaCorrRcb_px5->SetLineColor(kGreen+3);
 722     histoSigmaCorrRcb_px5->Draw("same");
 723     TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
 724     legrcb->SetFillColor(0);
 725     if (option==1) {
 726       legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rc/b=0.25","l");
 727       legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rc/b=0.5","l");
 728       legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rc/b=1.0","l");
 729       legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rc/b=1.5","l");
 730       legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rc/b=2.0","l");
 731       legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rc/b=3.0","l");
 732       legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rc/b=4.0","l");
 733     }else{
 734       legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rb=0.25","l");
 735       legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rb=0.5","l");
 736       legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rb=1.0","l");
 737       legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rb=1.5","l");
 738       legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rb=2.0","l");
 739       legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rb=3.0","l");
 740       legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rb=4.0","l");
 741     }
 742     legrcb->Draw();
 743     canvasSRcb1->Update();
 744   }
 745
 746
 747   //
 748   // Write the histograms to the output file
 749   //
 750   cout << " Saving the results ! " << endl<< endl;
 751
 752   out->cd();
 753   //
 754   hDirectMCpt->Write();        hFeedDownMCpt->Write();
 755   hDirectMCptMax->Write();     hDirectMCptMin->Write();
 756   hFeedDownMCptMax->Write();   hFeedDownMCptMin->Write();
 757   if(hDirectEffpt) hDirectEffpt->Write();        if(hFeedDownEffpt) hFeedDownEffpt->Write();
 758   hRECpt->Write();
 759   //
 760   histoYieldCorr->Write();
 761   histoYieldCorrMax->Write();     histoYieldCorrMin->Write();
 762   histoSigmaCorr->Write();
 763   histoSigmaCorrMax->Write();     histoSigmaCorrMin->Write();
 764
 765   if(PbPbEloss){
 766     histofcRcb->Write();    histofcRcb_px->Write();
 767     histoYieldCorrRcb->Write();
 768     histoSigmaCorrRcb->Write();
 769     nSigma->Write();
 770   }
 771
 772   if(asym){
 773     gYieldCorr->Write();
 774     gSigmaCorr->Write();
 775     if(gYieldCorrExtreme) gYieldCorrExtreme->Write();
 776     if(gSigmaCorrExtreme) gSigmaCorrExtreme->Write();
 777     if(gYieldCorrConservative) gYieldCorrConservative->Write();
 778     if(gSigmaCorrConservative) gSigmaCorrConservative->Write();
 779     if(asym && gFcConservative) gFcConservative->Write();
 780   }
 781
 782   if(option==1){
 783     histofc->Write();
 784     histofcMax->Write();     histofcMin->Write();
 785     if(asym && gFcExtreme) gFcExtreme->Write();
 786   }
 787
 788
 789   TH1D * hStatUncEffcSigma = spectra->GetDirectStatEffUncOnSigma();
 790   TH1D * hStatUncEffbSigma = spectra->GetFeedDownStatEffUncOnSigma();
 791   TH1D * hStatUncEffcFD = spectra->GetDirectStatEffUncOnFc();
 792   TH1D * hStatUncEffbFD = spectra->GetFeedDownStatEffUncOnFc();
 793   hStatUncEffcSigma->Write();
 794   hStatUncEffbSigma->Write();
 795   hStatUncEffcFD->Write();
 796   hStatUncEffbFD->Write();
 797
 798   // Draw the cross-section
 799   //  spectra->DrawSpectrum(gPrediction);
 800
 801   //  out->Close();
 802
 803 }