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