PWGHF/vertexingHF/macros/HFPtSpectrum.C

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