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