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.006245;
 145     } else if ( cc == kpPb2040 ) {
 146       tab = 0.134; tabUnc = 0.004899;
 147     } else if ( cc == kpPb4060 ) {
 148       tab = 0.092; tabUnc = 0.004796;
 149     } else if ( cc == kpPb60100 ) {
 150       tab = 0.041; tabUnc = 0.008832;
 151     }
 152   }
 153   else if( ccestimator == kZNA ){
 154     if ( cc == kpPb020 ) {
 155       tab = 0.164; tabUnc = 0.010724;
 156     } else if ( cc == kpPb2040 ) {
 157       tab = 0.137; tabUnc = 0.005099;
 158     } else if ( cc == kpPb4060 ) {
 159       tab = 0.1011; tabUnc = 0.006;
 160     } else if ( cc == kpPb60100 ) {
 161       tab = 0.0459; tabUnc = 0.003162;
 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     if(ccestimator==kV0A) {
 421       if(cc == kpPb020) systematics->SetCentrality("020V0A");
 422       else if(cc == kpPb2040) systematics->SetCentrality("2040V0A");
 423       else if(cc == kpPb4060) systematics->SetCentrality("4060V0A");
 424       else if(cc == kpPb60100) systematics->SetCentrality("60100V0A");
 425     } else if (ccestimator==kZNA) {
 426       if(cc == kpPb020) systematics->SetCentrality("020ZNA");
 427       else if(cc == kpPb2040) systematics->SetCentrality("2040ZNA");
 428       else if(cc == kpPb4060) systematics->SetCentrality("4060ZNA");
 429       else if(cc == kpPb60100) systematics->SetCentrality("60100ZNA");
 430     } else {
 431       if(!(cc == kpPb0100)) {
 432         cout <<" Error on the pPb options"<<endl;
 433         return;
 434       }
 435     }
 436   }
 437   //
 438   else if( cc!=kpp7 )  {
 439     systematics->SetCollisionType(1);
 440     if ( cc == k07half ) systematics->SetCentrality("07half");
 441     else if ( cc == k010 )  systematics->SetCentrality("010");
 442     else if ( cc == k1020 )  systematics->SetCentrality("1020");
 443     else if ( cc == k020 )  systematics->SetCentrality("020");
 444     else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
 445       systematics->SetCentrality("2040");
 446       systematics->SetIsPbPb2010EnergyScan(true);
 447     }
 448     else if ( cc == k3050 ) {
 449       if (isRaavsEP == kPhiIntegrated) systematics->SetCentrality("4080");
 450       else if (isRaavsEP == kInPlane) systematics->SetCentrality("3050InPlane");
 451       else if (isRaavsEP == kOutOfPlane) systematics->SetCentrality("3050OutOfPlane");
 452     }
 453     else if ( cc == k4060 || cc == k4050 || cc == k5060 )  systematics->SetCentrality("4060");
 454     else if ( cc == k6080 )  systematics->SetCentrality("6080");
 455     else if ( cc == k4080 ) systematics->SetCentrality("4080");
 456     else {
 457       cout << " Systematics not yet implemented " << endl;
 458       return;
 459     }
 460   } else { systematics->SetCollisionType(0); }
 461   //
 462   systematics->Init(decay);
 463   spectra->ComputeSystUncertainties(systematics,combineFeedDown);
 464
 465   //
 466   // Get the output histograms
 467   //
 468   // the corrected yield and cross-section
 469   histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum();
 470   histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum();
 471   histoYieldCorrMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum();
 472   histoYieldCorrMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum();
 473   histoSigmaCorrMax = (TH1D*)spectra->GetHistoUpperLimitCrossSectionFromYieldSpectrum();
 474   histoSigmaCorrMin = (TH1D*)spectra->GetHistoLowerLimitCrossSectionFromYieldSpectrum();
 475   histoYieldCorr->SetNameTitle("histoYieldCorr","corrected yield");
 476   histoYieldCorrMax->SetNameTitle("histoYieldCorrMax","max corrected yield");
 477   histoYieldCorrMin->SetNameTitle("histoYieldCorrMin","min corrected yield");
 478   histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section");
 479   histoSigmaCorrMax->SetNameTitle("histoSigmaCorrMax","max corrected invariant cross-section");
 480   histoSigmaCorrMin->SetNameTitle("histoSigmaCorrMin","min corrected invariant cross-section");
 481   // the efficiencies
 482   if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection();
 483   if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection();
 484   // Get the PbPb Eloss hypothesis histograms
 485   if(PbPbEloss){
 486     histofcRcb = spectra->GetHistoFeedDownCorrectionFcVsEloss();
 487     histoYieldCorrRcb = spectra->GetHistoFeedDownCorrectedSpectrumVsEloss();
 488     histoSigmaCorrRcb = spectra->GetHistoCrossSectionFromYieldSpectrumVsEloss();
 489     histofcRcb->SetName("histofcRcb");
 490     histoYieldCorrRcb->SetName("histoYieldCorrRcb");
 491     histoSigmaCorrRcb->SetName("histoSigmaCorrRcb");
 492   }
 493
 494   // Get & Rename the TGraphs
 495   gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
 496   gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
 497   if (asym) {
 498     gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme();
 499     gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme();
 500     gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative();
 501     gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative();
 502   }
 503
 504   // Get & Rename the TGraphs
 505   if (option==0){
 506     gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)");
 507     gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)");
 508   }
 509   if (option==1){
 510     // fc histos
 511     histofc = (TH1D*)spectra->GetHistoFeedDownCorrectionFc();
 512     histofcMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectionFc();
 513     histofcMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectionFc();
 514     histofc->SetNameTitle("histofc","fc correction factor");
 515     histofcMax->SetNameTitle("histofcMax","max fc correction factor");
 516     histofcMin->SetNameTitle("histofcMin","min fc correction factor");
 517     if (asym) {
 518       gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)");
 519       gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)");
 520       gFcExtreme = spectra->GetFeedDownCorrectionFcExtreme();
 521       gFcExtreme->SetNameTitle("gFcExtreme","gFcExtreme");
 522       gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by fc)");
 523       gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by fc)");
 524       gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
 525       gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
 526       gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by fc)");
 527       gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by fc)");
 528     }
 529   }
 530   if (option==2 && asym) {
 531     gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)");
 532     gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)");
 533     gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by Nb)");
 534     gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by Nb)");
 535     gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by Nb)");
 536     gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by Nb)");
 537     gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
 538     gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
 539   }
 540
 541   if(PbPbEloss){
 542     nSigma = spectra->GetNtupleCrossSectionVsEloss();
 543   }
 544
 545   //
 546   // Now, plot the results ! :)
 547   //
 548
 549   gROOT->SetStyle("Plain");
 550
 551   cout << " Drawing the results ! " << endl;
 552
 553   // control plots
 554   if (option==1) {
 555
 556     TCanvas *ceff = new TCanvas("ceff","efficiency drawing");
 557     ceff->Divide(1,2);
 558     ceff->cd(1);
 559     hDirectEffpt->Draw();
 560     ceff->cd(2);
 561     hFeedDownEffpt->Draw();
 562     ceff->Update();
 563
 564     TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
 565     cTheoryRebin->Divide(1,2);
 566     cTheoryRebin->cd(1);
 567     hDirectMCpt->Draw("");
 568     TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum();
 569     hDirectMCptRebin->SetLineColor(2);
 570     hDirectMCptRebin->Draw("same");
 571     cTheoryRebin->cd(2);
 572     hFeedDownMCpt->Draw("");
 573     TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum();
 574     hFeedDownRebin->SetLineColor(2);
 575     hFeedDownRebin->Draw("same");
 576     cTheoryRebin->Update();
 577
 578     TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
 579     cTheoryRebinLimits->Divide(1,2);
 580     cTheoryRebinLimits->cd(1);
 581     hDirectMCptMax->Draw("");
 582     TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum();
 583     hDirectMCptMaxRebin->SetLineColor(2);
 584     hDirectMCptMaxRebin->Draw("same");
 585     hDirectMCptMin->Draw("same");
 586     TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum();
 587     hDirectMCptMinRebin->SetLineColor(2);
 588     hDirectMCptMinRebin->Draw("same");
 589     cTheoryRebinLimits->cd(2);
 590     hFeedDownMCptMax->Draw("");
 591     TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum();
 592     hFeedDownMaxRebin->SetLineColor(2);
 593     hFeedDownMaxRebin->Draw("same");
 594     hFeedDownMCptMin->Draw("same");
 595     TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum();
 596     hFeedDownMinRebin->SetLineColor(2);
 597     hFeedDownMinRebin->Draw("same");
 598     cTheoryRebinLimits->Update();
 599   }
 600
 601   if (option==1) {
 602
 603     TCanvas * cfc = new TCanvas("cfc","Fc");
 604     histofcMax->Draw("c");
 605     histofc->Draw("csame");
 606     histofcMin->Draw("csame");
 607     cfc->Update();
 608
 609     if (asym) {
 610       TH2F *histofcDraw= new TH2F("histofcDraw","histofc (for drawing)",100,0,33.25,100,0.01,1.25);
 611       histofcDraw->SetStats(0);
 612       histofcDraw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 613       histofcDraw ->GetXaxis()->SetTitleSize(0.05);
 614       histofcDraw->GetXaxis()->SetTitleOffset(0.95);
 615       histofcDraw->GetYaxis()->SetTitle(" fc ");
 616       histofcDraw->GetYaxis()->SetTitleSize(0.05);
 617
 618       if (gFcExtreme){
 619
 620 //      for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 621 //        Double_t center=0., value=0.;
 622 //        gFcExtreme->GetPoint(item,center,value);
 623 //        Double_t highunc = gFcExtreme->GetErrorYhigh(item) / value ;
 624 //        Double_t lowunc = gFcExtreme->GetErrorYlow(item) / value ;
 625 //        cout << "Fc extreme: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 626 //      }
 627 //      for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 628 //        Double_t center=0., value=0.;
 629 //        gFcConservative->GetPoint(item,center,value);
 630 //        Double_t highunc = gFcConservative->GetErrorYhigh(item) / value ;
 631 //        Double_t lowunc = gFcConservative->GetErrorYlow(item) / value ;
 632 //        cout << "Fc conservative: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 633 //      }
 634         TCanvas *cfcExtreme = new TCanvas("cfcExtreme","Extreme Asymmetric fc (TGraphAsymmErr)");
 635         gFcExtreme->SetFillStyle(3006);
 636         gFcExtreme->SetLineWidth(3);
 637         gFcExtreme->SetMarkerStyle(20);
 638         gFcExtreme->SetFillColor(2);
 639         histofcDraw->Draw();
 640         gFcExtreme->Draw("3same");
 641
 642         if(gFcConservative){
 643           gFcConservative->SetFillStyle(3007);
 644           gFcConservative->SetFillColor(4);
 645           gFcConservative->Draw("3same");
 646         }
 647
 648         cfcExtreme->Update();
 649       }
 650     }
 651
 652   }
 653
 654   //
 655   // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra)
 656   //
 657   TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma");
 658   hDirectMCpt->SetMarkerStyle(20);
 659   hDirectMCpt->SetMarkerColor(4);
 660   hDirectMCpt->Draw("p");
 661   histoSigmaCorr->SetMarkerStyle(21);
 662   histoSigmaCorr->SetMarkerColor(2);
 663   histoSigmaCorr->Draw("psame");
 664   histoYieldCorr->SetMarkerStyle(22);
 665   histoYieldCorr->SetMarkerColor(6);
 666   histoYieldCorr->Draw("psame");
 667   hRECpt->SetMarkerStyle(23);
 668   hRECpt->SetMarkerColor(3);
 669   hRECpt->Draw("psame");
 670   cresult->SetLogy();
 671   cresult->Update();
 672
 673   TCanvas * cresult2 = new TCanvas("cresult2","corrected yield & sigma");
 674   histoSigmaCorr->SetMarkerStyle(21);
 675   histoSigmaCorr->SetMarkerColor(2);
 676   histoSigmaCorr->Draw("p");
 677   histoYieldCorr->SetMarkerStyle(22);
 678   histoYieldCorr->SetMarkerColor(6);
 679   histoYieldCorr->Draw("psame");
 680   hRECpt->SetMarkerStyle(23);
 681   hRECpt->SetMarkerColor(3);
 682   hRECpt->Draw("psame");
 683   cresult2->SetLogy();
 684   cresult2->Update();
 685
 686
 687   if (asym) {
 688
 689     TH2F *histoDraw = new TH2F("histoDraw","histo (for drawing)",100,0,33.25,100,50.,1e7);
 690     float max = 1.1*gYieldCorr->GetMaximum();
 691     histoDraw->SetAxisRange(0.1,max,"Y");
 692     histoDraw->SetStats(0);
 693     histoDraw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 694     histoDraw->GetXaxis()->SetTitleSize(0.05);
 695     histoDraw->GetXaxis()->SetTitleOffset(0.95);
 696     histoDraw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]");
 697     histoDraw->GetYaxis()->SetTitleSize(0.05);
 698     TCanvas * cyieldAsym = new TCanvas("cyieldAsym","Asymmetric corrected yield (TGraphAsymmErr)");
 699     gYieldCorr->SetFillStyle(3001);
 700     gYieldCorr->SetLineWidth(3);
 701     gYieldCorr->SetMarkerStyle(20);
 702     gYieldCorr->SetFillColor(3);
 703     histoDraw->Draw();
 704     gYieldCorr->Draw("3LPsame");
 705     gYieldCorr->Draw("Xsame");
 706     cyieldAsym->SetLogy();
 707     cyieldAsym->Update();
 708
 709     TCanvas * cyieldExtreme = new TCanvas("cyieldExtreme","Extreme Asymmetric corrected yield (TGraphAsymmErr)");
 710     histoYieldCorr->Draw();
 711     gYieldCorrExtreme->SetFillStyle(3002);
 712     gYieldCorrExtreme->SetLineWidth(3);
 713     gYieldCorrExtreme->SetMarkerStyle(20);
 714     gYieldCorrExtreme->SetFillColor(2);
 715     histoYieldCorr->Draw();
 716     gYieldCorr->Draw("3same");
 717     gYieldCorrExtreme->Draw("3same");
 718     cyieldExtreme->SetLogy();
 719     cyieldExtreme->Update();
 720
 721     TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9);
 722     max = 1.1*gSigmaCorr->GetMaximum();
 723     histo2Draw->SetAxisRange(0.1,max,"Y");
 724     histo2Draw->SetStats(0);
 725     histo2Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 726     histo2Draw->GetXaxis()->SetTitleSize(0.05);
 727     histo2Draw->GetXaxis()->SetTitleOffset(0.95);
 728     histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}");
 729     histo2Draw->GetYaxis()->SetTitleSize(0.05);
 730     TCanvas * csigmaAsym = new TCanvas("csigmaAsym","Asymmetric corrected sigma (TGraphAsymmErr)");
 731     gSigmaCorr->SetFillStyle(3001);
 732     gSigmaCorr->SetLineWidth(3);
 733     gSigmaCorr->SetMarkerStyle(21);
 734     gSigmaCorr->SetFillColor(3);
 735     histo2Draw->Draw();
 736     gSigmaCorr->Draw("3LPsame");
 737     gSigmaCorr->Draw("Xsame");
 738     csigmaAsym->SetLogy();
 739     csigmaAsym->Update();
 740
 741 //     cout << endl <<" Sytematics (stat approach) " <<endl;
 742 //     for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 743 //       Double_t center=0., value=0.;
 744 //       gSigmaCorr->GetPoint(item,center,value);
 745 //       Double_t highunc = gSigmaCorr->GetErrorYhigh(item) / value ;
 746 //       Double_t lowunc = gSigmaCorr->GetErrorYlow(item) / value ;
 747 //       cout << "Sigma syst (stat), i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 748 //     }
 749
 750     TCanvas * csigmaExtreme = new TCanvas("csigmaExtreme","Asymmetric extreme corrected sigma (TGraphAsymmErr)");
 751     histoSigmaCorr->Draw();
 752     gSigmaCorr->Draw("3Psame");
 753     gSigmaCorrExtreme->SetFillStyle(3002);
 754     gSigmaCorrExtreme->SetLineWidth(3);
 755     gSigmaCorrExtreme->SetMarkerStyle(21);
 756     gSigmaCorrExtreme->SetFillColor(2);
 757     gSigmaCorrExtreme->Draw("3Psame");
 758     csigmaExtreme->SetLogy();
 759     csigmaExtreme->Update();
 760
 761 //     cout << endl << " Sytematics (Extreme approach)" <<endl;
 762 //     for(Int_t item=0; item<gSigmaCorrExtreme->GetN(); item++){
 763 //       Double_t center=0., value=0.;
 764 //       gSigmaCorrExtreme->GetPoint(item,center,value);
 765 //       Double_t highunc = gSigmaCorrExtreme->GetErrorYhigh(item) / value ;
 766 //       Double_t lowunc = gSigmaCorrExtreme->GetErrorYlow(item) / value ;
 767 //       cout << "Sigma syst (extreme) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 768 //     }
 769
 770 //     cout << endl << " Sytematics (Conservative approach)" <<endl;
 771 //     for(Int_t item=0; item<gSigmaCorrConservative->GetN(); item++){
 772 //       Double_t center=0., value=0.;
 773 //       gSigmaCorrConservative->GetPoint(item,center,value);
 774 //       Double_t highunc = gSigmaCorrConservative->GetErrorYhigh(item) / value ;
 775 //       Double_t lowunc = gSigmaCorrConservative->GetErrorYlow(item) / value ;
 776 //       cout << "Sigma syst (conservative) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 777 //     }
 778
 779  }
 780
 781   // Draw the PbPb Eloss hypothesis histograms
 782   if(PbPbEloss){
 783     AliHFPtSpectrum *CalcBins=NULL;
 784     gStyle->SetPalette(1);
 785     TCanvas *canvasfcRcb = new TCanvas("canvasfcRcb","fc vs pt vs Rcb");
 786     //    histofcRcb->Draw("cont4z");
 787     histofcRcb->Draw("colz");
 788     canvasfcRcb->Update();
 789     canvasfcRcb->cd(2);
 790     TCanvas *canvasfcRcb1 = new TCanvas("canvasfcRcb1","fc vs pt vs Rcb=1");
 791     histofcRcb_px = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px",40,40);
 792     histofcRcb_px->SetLineColor(2);
 793     if (option==1) {
 794       histofc->Draw();
 795       histofcRcb_px->Draw("same");
 796     } else histofcRcb_px->Draw("");
 797     canvasfcRcb1->Update();
 798     TCanvas *canvasfcRcb2 = new TCanvas("canvasfcRcb2","fc vs pt vs Rcb fixed Rcb");
 799     Int_t bin0 = CalcBins->FindTH2YBin(histofcRcb,0.25);
 800     Int_t bin1 = CalcBins->FindTH2YBin(histofcRcb,0.5);
 801     Int_t bin2 = CalcBins->FindTH2YBin(histofcRcb,1.0);
 802     Int_t bin3 = CalcBins->FindTH2YBin(histofcRcb,1.5);
 803     Int_t bin4 = CalcBins->FindTH2YBin(histofcRcb,2.0);
 804     Int_t bin5 = CalcBins->FindTH2YBin(histofcRcb,3.0);
 805     Int_t bin6 = CalcBins->FindTH2YBin(histofcRcb,4.0);
 806     TH1D * histofcRcb_px0a = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0a",bin0,bin0);
 807     TH1D * histofcRcb_px0 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0",bin1,bin1);
 808     TH1D * histofcRcb_px1 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px1",bin2,bin2);
 809     TH1D * histofcRcb_px2 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px2",bin3,bin3);
 810     TH1D * histofcRcb_px3 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px3",bin4,bin4);
 811     TH1D * histofcRcb_px4 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px4",bin5,bin5);
 812     TH1D * histofcRcb_px5 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px5",bin6,bin6);
 813     if (option==1) {
 814       histofc->Draw();
 815       //      histofcRcb_px->Draw("same");
 816     } else {
 817       //      histofcRcb_px->Draw("");
 818       histofcRcb_px0a->SetLineColor(2);
 819       histofcRcb_px0a->Draw("");
 820     }
 821     histofcRcb_px0a->SetLineColor(2);
 822     histofcRcb_px0a->Draw("same");
 823     histofcRcb_px0->SetLineColor(4);
 824     histofcRcb_px0->Draw("same");
 825     histofcRcb_px1->SetLineColor(3);
 826     histofcRcb_px1->Draw("same");
 827     histofcRcb_px2->SetLineColor(kCyan);
 828     histofcRcb_px2->Draw("same");
 829     histofcRcb_px3->SetLineColor(kMagenta+1);
 830     histofcRcb_px3->Draw("same");
 831     histofcRcb_px4->SetLineColor(kOrange+7);
 832     histofcRcb_px4->Draw("same");
 833     histofcRcb_px5->SetLineColor(kGreen+3);
 834     histofcRcb_px5->Draw("same");
 835     TLegend *legrcc = new TLegend(0.8,0.8,0.95,0.9);
 836     legrcc->SetFillColor(0);
 837     if (option==1) {
 838       legrcc->AddEntry(histofcRcb_px0a,"Rc/b=0.25","l");
 839       legrcc->AddEntry(histofcRcb_px0,"Rc/b=0.5","l");
 840       legrcc->AddEntry(histofcRcb_px1,"Rc/b=1.0","l");
 841       legrcc->AddEntry(histofcRcb_px2,"Rc/b=1.5","l");
 842       legrcc->AddEntry(histofcRcb_px3,"Rc/b=2.0","l");
 843       legrcc->AddEntry(histofcRcb_px4,"Rc/b=3.0","l");
 844       legrcc->AddEntry(histofcRcb_px5,"Rc/b=4.0","l");
 845     }else{
 846       legrcc->AddEntry(histofcRcb_px0a,"Rb=0.25","l");
 847       legrcc->AddEntry(histofcRcb_px0,"Rb=0.5","l");
 848       legrcc->AddEntry(histofcRcb_px1,"Rb=1.0","l");
 849       legrcc->AddEntry(histofcRcb_px2,"Rb=1.5","l");
 850       legrcc->AddEntry(histofcRcb_px3,"Rb=2.0","l");
 851       legrcc->AddEntry(histofcRcb_px4,"Rb=3.0","l");
 852       legrcc->AddEntry(histofcRcb_px5,"Rb=4.0","l");
 853     }
 854     legrcc->Draw();
 855     canvasfcRcb2->Update();
 856     TCanvas *canvasYRcb = new TCanvas("canvasYRcb","corrected yield vs pt vs Rcb");
 857     histoYieldCorrRcb->Draw("cont4z");
 858     canvasYRcb->Update();
 859     TCanvas *canvasSRcb = new TCanvas("canvasSRcb","sigma vs pt vs Rcb");
 860     histoSigmaCorrRcb->Draw("cont4z");
 861     canvasSRcb->Update();
 862     TCanvas *canvasSRcb1 = new TCanvas("canvasSRcb1","sigma vs pt vs Rcb fixed Rcb");
 863     TH1D * histoSigmaCorrRcb_px0a = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0a",bin0,bin0);
 864     TH1D * histoSigmaCorrRcb_px0 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0",bin1,bin1);
 865     TH1D * histoSigmaCorrRcb_px1 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px1",bin2,bin2);
 866     TH1D * histoSigmaCorrRcb_px2 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px2",bin3,bin3);
 867     TH1D * histoSigmaCorrRcb_px3 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px3",bin4,bin4);
 868     TH1D * histoSigmaCorrRcb_px4 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px4",bin5,bin5);
 869     TH1D * histoSigmaCorrRcb_px5 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px5",bin6,bin6);
 870     histoSigmaCorr->Draw();
 871     histoSigmaCorrRcb_px0a->SetLineColor(2);
 872     histoSigmaCorrRcb_px0a->Draw("hsame");
 873     histoSigmaCorrRcb_px0->SetLineColor(4);
 874     histoSigmaCorrRcb_px0->Draw("hsame");
 875     histoSigmaCorrRcb_px1->SetLineColor(3);
 876     histoSigmaCorrRcb_px1->Draw("hsame");
 877     histoSigmaCorrRcb_px2->SetLineColor(kCyan);
 878     histoSigmaCorrRcb_px2->Draw("hsame");
 879     histoSigmaCorrRcb_px3->SetLineColor(kMagenta+1);
 880     histoSigmaCorrRcb_px3->Draw("hsame");
 881     histoSigmaCorrRcb_px4->SetLineColor(kOrange+7);
 882     histoSigmaCorrRcb_px4->Draw("same");
 883     histoSigmaCorrRcb_px5->SetLineColor(kGreen+3);
 884     histoSigmaCorrRcb_px5->Draw("same");
 885     TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
 886     legrcb->SetFillColor(0);
 887     if (option==1) {
 888       legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rc/b=0.25","l");
 889       legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rc/b=0.5","l");
 890       legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rc/b=1.0","l");
 891       legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rc/b=1.5","l");
 892       legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rc/b=2.0","l");
 893       legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rc/b=3.0","l");
 894       legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rc/b=4.0","l");
 895     }else{
 896       legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rb=0.25","l");
 897       legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rb=0.5","l");
 898       legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rb=1.0","l");
 899       legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rb=1.5","l");
 900       legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rb=2.0","l");
 901       legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rb=3.0","l");
 902       legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rb=4.0","l");
 903     }
 904     legrcb->Draw();
 905     canvasSRcb1->Update();
 906   }
 907
 908
 909   //
 910   // Write the histograms to the output file
 911   //
 912   cout << " Saving the results ! " << endl<< endl;
 913
 914   out->cd();
 915   //
 916   hDirectMCpt->Write();        hFeedDownMCpt->Write();
 917   hDirectMCptMax->Write();     hDirectMCptMin->Write();
 918   hFeedDownMCptMax->Write();   hFeedDownMCptMin->Write();
 919   if(hDirectEffpt) hDirectEffpt->Write();        if(hFeedDownEffpt) hFeedDownEffpt->Write();
 920   hRECpt->Write();
 921   //
 922   histoYieldCorr->Write();
 923   histoYieldCorrMax->Write();     histoYieldCorrMin->Write();
 924   histoSigmaCorr->Write();
 925   histoSigmaCorrMax->Write();     histoSigmaCorrMin->Write();
 926
 927   if(PbPbEloss){
 928     histofcRcb->Write();    histofcRcb_px->Write();
 929     histoYieldCorrRcb->Write();
 930     histoSigmaCorrRcb->Write();
 931     nSigma->Write();
 932   }
 933
 934   gYieldCorr->Write();
 935   gSigmaCorr->Write();
 936   if(asym){
 937     if(gYieldCorrExtreme) gYieldCorrExtreme->Write();
 938     if(gSigmaCorrExtreme) gSigmaCorrExtreme->Write();
 939     if(gYieldCorrConservative) gYieldCorrConservative->Write();
 940     if(gSigmaCorrConservative) gSigmaCorrConservative->Write();
 941     if(asym && gFcConservative) gFcConservative->Write();
 942   }
 943
 944   if(option==1){
 945     histofc->Write();
 946     histofcMax->Write();     histofcMin->Write();
 947     if(asym && gFcExtreme) gFcExtreme->Write();
 948   }
 949
 950
 951   TH1D * hStatUncEffcSigma = spectra->GetDirectStatEffUncOnSigma();
 952   TH1D * hStatUncEffbSigma = spectra->GetFeedDownStatEffUncOnSigma();
 953   hStatUncEffcSigma->Write();
 954   hStatUncEffbSigma->Write();
 955   if(option!=0){
 956     TH1D * hStatUncEffcFD = spectra->GetDirectStatEffUncOnFc();
 957     TH1D * hStatUncEffbFD = spectra->GetFeedDownStatEffUncOnFc();
 958     hStatUncEffcFD->Write();
 959     hStatUncEffbFD->Write();
 960   }
 961
 962   // Draw the cross-section
 963   //  spectra->DrawSpectrum(gPrediction);
 964
 965   //  out->Close();
 966
 967 }