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