PWG3/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 "TFile.h"
  13 #include "TGraphAsymmErrors.h"
  14 #include "TCanvas.h"
  15 #include "TROOT.h"
  16
  17 #include "AliHFSystErr.h"
  18 #include "AliHFPtSpectrum.h"
  19
  20 //
  21 // Macro execution parameters:
  22 //  0) filename with the theoretical predictions  (direct & feed-down)
  23 //  1) acceptance and reconstruction efficiencies file name (direct & feed-down)
  24 //  2) reconstructed spectra file name
  25 //  3) output file name
  26 //  4) Set the feed-down calculation option flag: 0=none, 1=fc only, 2=Nb only
  27 //  5) Set the luminosity
  28 //  6) Set the trigger efficiency
  29 //  7-14) 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 this parameters
  30 //
  31 void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
  32                     const char *efffilename="Efficiencies.root",
  33                     const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0",
  34                     const char *outfilename="HFPtSpectrum.root",
  35                     Int_t option=1, Double_t lumi=1.0, Double_t effTrig=1.0,
  36                     const char *directsimufilename="", const char *directsimuhistoname="CFHFccontainer0_New_3Prong_SelStep0_proj-pt",
  37                     const char *directrecofilename="", const char *directrecohistoname="CFHFccontainer0_New_3Prong_SelStep8_proj-pt",
  38                     const char *feeddownsimufilename="", const char *feeddownsimuhistoname="CFHFccontainer0allD_New_3Prong_SelStep0_proj-pt",
  39                     const char *feeddownrecofilename="", const char *feeddownrecohistoname="CFHFccontainer0allD_New_3Prong_SelStep8_proj-pt") {
  40
  41   //  Set if calculation considers asymmetric uncertainties or not
  42   bool asym = true;
  43
  44   // Set the meson and decay
  45   // (only D0 -> K pi, D+--> K pi pi & D* --> D0 pi implemented here)
  46   bool isD0Kpi = true;
  47   bool isDplusKpipi = false;
  48   bool isDstarD0pi = false;
  49   if (isD0Kpi && isDplusKpipi && isDstarD0pi) {
  50     cout << "Sorry, can not deal with more than one correction at the same time"<<endl;
  51     return;
  52   }
  53
  54   if (option>2) {
  55     cout<< "Bad calculation option, should be <=2"<<endl;
  56     return;
  57   }
  58   if (option==0) asym = false;
  59
  60   //
  61   // Get the histograms from the files
  62   //
  63   TH1D *hDirectMCpt=0;           // Input MC c-->D spectra
  64   TH1D *hFeedDownMCpt=0;         // Input MC b-->D spectra
  65   TH1D *hDirectMCptMax=0;        // Input MC maximum c-->D spectra
  66   TH1D *hDirectMCptMin=0;        // Input MC minimum c-->D spectra
  67   TH1D *hFeedDownMCptMax=0;      // Input MC maximum b-->D spectra
  68   TH1D *hFeedDownMCptMin=0;      // Input MC minimum b-->D spectra
  69   TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra
  70   TH1D *hDirectEffpt=0;          // c-->D Acceptance and efficiency correction
  71   TH1D *hFeedDownEffpt=0;        // b-->D Acceptance and efficiency correction
  72   TH1D *hRECpt=0;                // all reconstructed D
  73   //
  74   TH1D *hDirectSimulationpt=0;       // Simulated c--D spectra (used to re-compute the efficiency)
  75   TH1D *hDirectReconstructionpt=0;   // Reconstructed c--D spectra (used to re-compute the efficiency)
  76   TH1D *hFeedDownSimulationpt=0;     // Simulated b--D spectra (used to re-compute the efficiency)
  77   TH1D *hFeedDownReconstructionpt=0; // Reconstructed b--D spectra (used to re-compute the efficiency)
  78   //
  79
  80   //
  81   // Define/Get the input histograms
  82   //
  83   Int_t decay=0;
  84   TFile * mcfile = new TFile(mcfilename,"read");
  85   if (isD0Kpi){
  86     decay = 1;
  87     hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central");
  88     hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central");
  89     hDirectMCptMax = (TH1D*)mcfile->Get("hD0Kpipred_max");
  90     hDirectMCptMin = (TH1D*)mcfile->Get("hD0Kpipred_min");
  91     hFeedDownMCptMax = (TH1D*)mcfile->Get("hD0KpifromBpred_max");
  92     hFeedDownMCptMin = (TH1D*)mcfile->Get("hD0KpifromBpred_min");
  93     gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction");
  94   }
  95   else if (isDplusKpipi){
  96     decay = 2;
  97     hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central");
  98     hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central");
  99     hDirectMCptMax = (TH1D*)mcfile->Get("hDpluskpipipred_max");
 100     hDirectMCptMin = (TH1D*)mcfile->Get("hDpluskpipipred_min");
 101     hFeedDownMCptMax = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max");
 102     hFeedDownMCptMin = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min");
 103     gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction");
 104   }
 105   else if(isDstarD0pi){
 106     decay = 3;
 107     hDirectMCpt = (TH1D*)mcfile->Get("hDstarD0pipred_central");
 108     hFeedDownMCpt = (TH1D*)mcfile->Get("hDstarD0pifromBpred_central");
 109     hDirectMCptMax = (TH1D*)mcfile->Get("hDstarD0pipred_max");
 110     hDirectMCptMin = (TH1D*)mcfile->Get("hDstarD0pipred_min");
 111     hFeedDownMCptMax = (TH1D*)mcfile->Get("hDstarD0pifromBpred_max");
 112     hFeedDownMCptMin = (TH1D*)mcfile->Get("hDstarD0pifromBpred_min");
 113     gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction");
 114   }
 115   //
 116   hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra");
 117   hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra");
 118   hDirectMCptMax->SetNameTitle("hDirectMCptMax","max direct MC spectra");
 119   hDirectMCptMin->SetNameTitle("hDirectMCptMin","min direct MC spectra");
 120   hFeedDownMCptMax->SetNameTitle("hFeedDownMCptMax","max feed-down MC spectra");
 121   hFeedDownMCptMin->SetNameTitle("hFeedDownMCptMin","min feed-down MC spectra");
 122   //
 123   //
 124   if (strcmp(directsimufilename,"")!=0 && strcmp(directrecofilename,"")!=0 &&
 125       strcmp(feeddownsimufilename,"")!=0 && strcmp(feeddownrecofilename,"")!=0 ) {
 126     if (strcmp(directsimufilename,"")!=0){
 127       TFile *directSimufile = new TFile(directsimufilename,"read");
 128       hDirectSimulationpt = (TH1D*)directSimufile->Get(directsimuhistoname);
 129       hDirectSimulationpt->SetNameTitle("hDirectSimulationpt","hDirectSimulationpt");
 130     }
 131     if (strcmp(directrecofilename,"")!=0){
 132       TFile *directRecofile = new TFile(directrecofilename,"read");
 133       hDirectReconstructionpt = (TH1D*)directRecofile->Get(directrecohistoname);
 134       hDirectReconstructionpt->SetNameTitle("hDirectReconstructionpt","hDirectReconstructionpt");
 135     }
 136     if (strcmp(feeddownsimufilename,"")!=0){
 137       TFile *feeddownSimufile = new TFile(feeddownsimufilename,"read");
 138       hFeedDownSimulationpt = (TH1D*)feeddownSimufile->Get(feeddownsimuhistoname);
 139       hFeedDownSimulationpt->SetNameTitle("hFeedDownSimulationpt","hFeedDownSimulationpt");
 140     }
 141     if (strcmp(feeddownrecofilename,"")!=0){
 142       TFile *feeddownRecofile = new TFile(feeddownrecofilename,"read");
 143       hFeedDownReconstructionpt = (TH1D*)feeddownRecofile->Get(feeddownrecohistoname);
 144       hFeedDownReconstructionpt->SetNameTitle("hFeedDownReconstructionpt","hFeedDownReconstructionpt");
 145     }
 146   }
 147   else {
 148     TFile * efffile = new TFile(efffilename,"read");
 149     hDirectEffpt = (TH1D*)efffile->Get("hEffD");
 150     hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
 151     hFeedDownEffpt = (TH1D*)efffile->Get("hEffB");
 152     hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff");
 153   }
 154   //
 155   //
 156   TFile * recofile = new TFile(recofilename,"read");
 157   hRECpt = (TH1D*)recofile->Get(recohistoname);
 158   hRECpt->SetNameTitle("hRECpt","Reconstructed spectra");
 159
 160   //
 161   // Define the output histograms
 162   //
 163   TFile *out = new TFile(outfilename,"recreate");
 164   //
 165   TH1D *histofc=0;
 166   TH1D *histofcMax=0;
 167   TH1D *histofcMin=0;
 168   TH1D *histoYieldCorr=0;
 169   TH1D *histoYieldCorrMax=0;
 170   TH1D *histoYieldCorrMin=0;
 171   TH1D *histoSigmaCorr=0;
 172   TH1D *histoSigmaCorrMax=0;
 173   TH1D *histoSigmaCorrMin=0;
 174   //
 175   int nbins = hRECpt->GetNbinsX();
 176   TGraphAsymmErrors * gYieldCorr = 0;
 177   TGraphAsymmErrors * gSigmaCorr = 0;
 178   TGraphAsymmErrors * gFcExtreme = 0;
 179   TGraphAsymmErrors * gFcConservative = 0;
 180   TGraphAsymmErrors * gYieldCorrExtreme = 0;
 181   TGraphAsymmErrors * gSigmaCorrExtreme = 0;
 182   TGraphAsymmErrors * gYieldCorrConservative = 0;
 183   TGraphAsymmErrors * gSigmaCorrConservative = 0;
 184
 185
 186   //
 187   // Main functionalities for the calculation
 188   //
 189
 190   // Define and set the basic option flags
 191   AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option);
 192   spectra->SetFeedDownCalculationOption(option);
 193   spectra->SetComputeAsymmetricUncertainties(asym);
 194
 195   // Feed the input histograms
 196   //  reconstructed spectra
 197   cout << " Setting the reconstructed spectrum,";
 198   spectra->SetReconstructedSpectrum(hRECpt);
 199   // acceptance and efficiency corrections
 200   cout << " the files to compute the efficiency,";
 201   if (strcmp(directsimufilename,"")!=0 && strcmp(directrecofilename,"")!=0 &&
 202       strcmp(feeddownsimufilename,"")!=0 && strcmp(feeddownrecofilename,"")!=0 ) {
 203     spectra->EstimateAndSetDirectEfficiencyRecoBin(hDirectSimulationpt,hDirectReconstructionpt);
 204     spectra->EstimateAndSetFeedDownEfficiencyRecoBin(hFeedDownSimulationpt,hFeedDownReconstructionpt);
 205   }
 206   else {
 207     spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt);
 208   }
 209   // option specific histos (theory)
 210   cout << " the theoretical spectra";
 211   if(option==1){
 212     spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt);
 213     if(asym) spectra->SetMCptDistributionsBounds(hDirectMCptMax,hDirectMCptMin,hFeedDownMCptMax,hFeedDownMCptMin);
 214   }
 215   else if(option==2){
 216     spectra->SetFeedDownMCptSpectra(hFeedDownMCpt);
 217     if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCptMax,hFeedDownMCptMin);
 218   }
 219
 220   cout << " and the normalization" <<endl;
 221   // Set normalization factors (uncertainties set to 0. as example)
 222   double lumiUnc = 0.10*lumi; // 10% uncertainty on the luminosity
 223   spectra->SetLuminosity(lumi,lumiUnc);
 224   spectra->SetTriggerEfficiency(effTrig,0.);
 225
 226   // Set the global uncertainties on the efficiencies (in percent)
 227   double globalEffUnc = 0.15;
 228   double globalBCEffRatioUnc = 0.15;
 229   //  double globalEffUnc = 0.;
 230   //  double globalBCEffRatioUnc = 0.;
 231   spectra->SetAccEffPercentageUncertainty(globalEffUnc,globalBCEffRatioUnc);
 232
 233   // Do the calculations
 234   cout << " Doing the calculation... "<< endl;
 235   double deltaY = 1.0;
 236   double branchingRatioC = 1.0;
 237   double branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction
 238   spectra->ComputeHFPtSpectrum(deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
 239   cout << "   ended the calculation, getting the histograms back " << endl;
 240
 241   // Set the systematics externally
 242   AliHFSystErr *systematics = new AliHFSystErr();
 243   systematics->Init(decay);
 244   bool combineFeedDown = true;
 245   spectra->ComputeSystUncertainties(systematics,combineFeedDown);
 246
 247
 248   //
 249   // Get the output histograms
 250   //
 251   // the corrected yield and cross-section
 252   histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum();
 253   histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum();
 254   histoYieldCorrMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum();
 255   histoYieldCorrMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum();
 256   histoSigmaCorrMax = (TH1D*)spectra->GetHistoUpperLimitCrossSectionFromYieldSpectrum();
 257   histoSigmaCorrMin = (TH1D*)spectra->GetHistoLowerLimitCrossSectionFromYieldSpectrum();
 258   histoYieldCorr->SetNameTitle("histoYieldCorr","corrected yield");
 259   histoYieldCorrMax->SetNameTitle("histoYieldCorrMax","max corrected yield");
 260   histoYieldCorrMin->SetNameTitle("histoYieldCorrMin","min corrected yield");
 261   histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section");
 262   histoSigmaCorrMax->SetNameTitle("histoSigmaCorrMax","max corrected invariant cross-section");
 263   histoSigmaCorrMin->SetNameTitle("histoSigmaCorrMin","min corrected invariant cross-section");
 264   // the efficiencies
 265   if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection();
 266   if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection();
 267
 268   // Get & Rename the TGraphs
 269   if (asym) {
 270     gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
 271     gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
 272     gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme();
 273     gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme();
 274     gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative();
 275     gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative();
 276   }
 277
 278   // Get & Rename the TGraphs
 279   if (option==0 && asym){
 280     gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)");
 281     gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)");
 282   }
 283   if (option==1){
 284     // fc histos
 285     histofc = (TH1D*)spectra->GetHistoFeedDownCorrectionFc();
 286     histofcMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectionFc();
 287     histofcMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectionFc();
 288     histofc->SetNameTitle("histofc","fc correction factor");
 289     histofcMax->SetNameTitle("histofcMax","max fc correction factor");
 290     histofcMin->SetNameTitle("histofcMin","min fc correction factor");
 291     if (asym) {
 292       gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)");
 293       gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)");
 294       gFcExtreme = spectra->GetFeedDownCorrectionFcExtreme();
 295       gFcExtreme->SetNameTitle("gFcExtreme","gFcExtreme");
 296       gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by fc)");
 297       gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by fc)");
 298       gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
 299       gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
 300       gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by fc)");
 301       gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by fc)");
 302     }
 303   }
 304   if (option==2 && asym) {
 305     gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)");
 306     gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)");
 307     gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by Nb)");
 308     gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by Nb)");
 309     gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by Nb)");
 310     gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by Nb)");
 311     gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
 312     gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
 313   }
 314
 315   //
 316   // Now, plot the results ! :)
 317   //
 318
 319   cout << " Drawing the results ! " << endl;
 320
 321   // control plots
 322   if (option==1) {
 323
 324     TCanvas *ceff = new TCanvas("ceff","efficiency drawing");
 325     ceff->Divide(1,2);
 326     ceff->cd(1);
 327     hDirectEffpt->Draw();
 328     ceff->cd(2);
 329     hFeedDownEffpt->Draw();
 330     ceff->Update();
 331
 332     TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
 333     cTheoryRebin->Divide(1,2);
 334     cTheoryRebin->cd(1);
 335     hDirectMCpt->Draw("");
 336     TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum();
 337     hDirectMCptRebin->SetLineColor(2);
 338     hDirectMCptRebin->Draw("same");
 339     cTheoryRebin->cd(2);
 340     hFeedDownMCpt->Draw("");
 341     TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum();
 342     hFeedDownRebin->SetLineColor(2);
 343     hFeedDownRebin->Draw("same");
 344     cTheoryRebin->Update();
 345
 346     TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
 347     cTheoryRebinLimits->Divide(1,2);
 348     cTheoryRebinLimits->cd(1);
 349     hDirectMCptMax->Draw("");
 350     TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum();
 351     hDirectMCptMaxRebin->SetLineColor(2);
 352     hDirectMCptMaxRebin->Draw("same");
 353     hDirectMCptMin->Draw("same");
 354     TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum();
 355     hDirectMCptMinRebin->SetLineColor(2);
 356     hDirectMCptMinRebin->Draw("same");
 357     cTheoryRebinLimits->cd(2);
 358     hFeedDownMCptMax->Draw("");
 359     TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum();
 360     hFeedDownMaxRebin->SetLineColor(2);
 361     hFeedDownMaxRebin->Draw("same");
 362     hFeedDownMCptMin->Draw("same");
 363     TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum();
 364     hFeedDownMinRebin->SetLineColor(2);
 365     hFeedDownMinRebin->Draw("same");
 366     cTheoryRebinLimits->Update();
 367   }
 368
 369   if (option==1) {
 370
 371     TCanvas * cfc = new TCanvas("cfc","Fc");
 372     histofcMax->Draw("c");
 373     histofc->Draw("csame");
 374     histofcMin->Draw("csame");
 375     cfc->Update();
 376
 377     if (asym) {
 378       TH2F *histofcDraw= new TH2F("histofcDraw","histofc (for drawing)",100,0,33.25,100,0.01,1.25);
 379       histofcDraw->SetStats(0);
 380       histofcDraw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 381       histofcDraw ->GetXaxis()->SetTitleSize(0.05);
 382       histofcDraw->GetXaxis()->SetTitleOffset(0.95);
 383       histofcDraw->GetYaxis()->SetTitle(" fc ");
 384       histofcDraw->GetYaxis()->SetTitleSize(0.05);
 385
 386       if (gFcExtreme){
 387
 388 //      for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 389 //        Double_t center=0., value=0.;
 390 //        gFcExtreme->GetPoint(item,center,value);
 391 //        Double_t highunc = gFcExtreme->GetErrorYhigh(item) / value ;
 392 //        Double_t lowunc = gFcExtreme->GetErrorYlow(item) / value ;
 393 //        cout << "Fc extreme: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 394 //      }
 395 //      for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 396 //        Double_t center=0., value=0.;
 397 //        gFcConservative->GetPoint(item,center,value);
 398 //        Double_t highunc = gFcConservative->GetErrorYhigh(item) / value ;
 399 //        Double_t lowunc = gFcConservative->GetErrorYlow(item) / value ;
 400 //        cout << "Fc conservative: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 401 //      }
 402         TCanvas *cfcExtreme = new TCanvas("cfcExtreme","Extreme Asymmetric fc (TGraphAsymmErr)");
 403         gFcExtreme->SetFillStyle(3006);
 404         gFcExtreme->SetLineWidth(3);
 405         gFcExtreme->SetMarkerStyle(20);
 406         gFcExtreme->SetFillColor(2);
 407         histofcDraw->Draw();
 408         gFcExtreme->Draw("3same");
 409
 410         if(gFcConservative){
 411           gFcConservative->SetFillStyle(3007);
 412           gFcConservative->SetFillColor(4);
 413           gFcConservative->Draw("3same");
 414         }
 415
 416         cfcExtreme->Update();
 417       }
 418     }
 419
 420   }
 421
 422   //
 423   // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra)
 424   //
 425   TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma");
 426   hDirectMCpt->SetMarkerStyle(20);
 427   hDirectMCpt->SetMarkerColor(4);
 428   hDirectMCpt->Draw("p");
 429   histoSigmaCorr->SetMarkerStyle(21);
 430   histoSigmaCorr->SetMarkerColor(2);
 431   histoSigmaCorr->Draw("psame");
 432   histoYieldCorr->SetMarkerStyle(22);
 433   histoYieldCorr->SetMarkerColor(6);
 434   histoYieldCorr->Draw("psame");
 435   hRECpt->SetMarkerStyle(23);
 436   hRECpt->SetMarkerColor(3);
 437   hRECpt->Draw("psame");
 438   cresult->SetLogy();
 439   cresult->Update();
 440
 441   TCanvas * cresult2 = new TCanvas("cresult2","corrected yield & sigma");
 442   histoSigmaCorr->SetMarkerStyle(21);
 443   histoSigmaCorr->SetMarkerColor(2);
 444   histoSigmaCorr->Draw("p");
 445   histoYieldCorr->SetMarkerStyle(22);
 446   histoYieldCorr->SetMarkerColor(6);
 447   histoYieldCorr->Draw("psame");
 448   hRECpt->SetMarkerStyle(23);
 449   hRECpt->SetMarkerColor(3);
 450   hRECpt->Draw("psame");
 451   cresult2->SetLogy();
 452   cresult2->Update();
 453
 454
 455   if (asym) {
 456
 457     TH2F *histoDraw = new TH2F("histoDraw","histo (for drawing)",100,0,33.25,100,50.,1e7);
 458     float max = 1.1*gYieldCorr->GetMaximum();
 459     histoDraw->SetAxisRange(0.1,max,"Y");
 460     histoDraw->SetStats(0);
 461     histoDraw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 462     histoDraw->GetXaxis()->SetTitleSize(0.05);
 463     histoDraw->GetXaxis()->SetTitleOffset(0.95);
 464     histoDraw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]");
 465     histoDraw->GetYaxis()->SetTitleSize(0.05);
 466     TCanvas * cyieldAsym = new TCanvas("cyieldAsym","Asymmetric corrected yield (TGraphAsymmErr)");
 467     gYieldCorr->SetFillStyle(3001);
 468     gYieldCorr->SetLineWidth(3);
 469     gYieldCorr->SetMarkerStyle(20);
 470     gYieldCorr->SetFillColor(3);
 471     histoDraw->Draw();
 472     gYieldCorr->Draw("3LPsame");
 473     gYieldCorr->Draw("Xsame");
 474     cyieldAsym->SetLogy();
 475     cyieldAsym->Update();
 476
 477     TCanvas * cyieldExtreme = new TCanvas("cyieldExtreme","Extreme Asymmetric corrected yield (TGraphAsymmErr)");
 478     histoYieldCorr->Draw();
 479     gYieldCorrExtreme->SetFillStyle(3002);
 480     gYieldCorrExtreme->SetLineWidth(3);
 481     gYieldCorrExtreme->SetMarkerStyle(20);
 482     gYieldCorrExtreme->SetFillColor(2);
 483     histoYieldCorr->Draw();
 484     gYieldCorr->Draw("3same");
 485     gYieldCorrExtreme->Draw("3same");
 486     cyieldExtreme->SetLogy();
 487     cyieldExtreme->Update();
 488
 489     TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9);
 490     max = 1.1*gSigmaCorr->GetMaximum();
 491     histo2Draw->SetAxisRange(0.1,max,"Y");
 492     histo2Draw->SetStats(0);
 493     histo2Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 494     histo2Draw->GetXaxis()->SetTitleSize(0.05);
 495     histo2Draw->GetXaxis()->SetTitleOffset(0.95);
 496     histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}");
 497     histo2Draw->GetYaxis()->SetTitleSize(0.05);
 498     TCanvas * csigmaAsym = new TCanvas("csigmaAsym","Asymmetric corrected sigma (TGraphAsymmErr)");
 499     gSigmaCorr->SetFillStyle(3001);
 500     gSigmaCorr->SetLineWidth(3);
 501     gSigmaCorr->SetMarkerStyle(21);
 502     gSigmaCorr->SetFillColor(3);
 503     histo2Draw->Draw();
 504     gSigmaCorr->Draw("3LPsame");
 505     gSigmaCorr->Draw("Xsame");
 506     csigmaAsym->SetLogy();
 507     csigmaAsym->Update();
 508
 509 //     cout << endl <<" Sytematics (stat approach) " <<endl;
 510 //     for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
 511 //       Double_t center=0., value=0.;
 512 //       gSigmaCorr->GetPoint(item,center,value);
 513 //       Double_t highunc = gSigmaCorr->GetErrorYhigh(item) / value ;
 514 //       Double_t lowunc = gSigmaCorr->GetErrorYlow(item) / value ;
 515 //       cout << "Sigma syst (stat), i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 516 //     }
 517
 518     TCanvas * csigmaExtreme = new TCanvas("csigmaExtreme","Asymmetric extreme corrected sigma (TGraphAsymmErr)");
 519     histoSigmaCorr->Draw();
 520     gSigmaCorr->Draw("3Psame");
 521     gSigmaCorrExtreme->SetFillStyle(3002);
 522     gSigmaCorrExtreme->SetLineWidth(3);
 523     gSigmaCorrExtreme->SetMarkerStyle(21);
 524     gSigmaCorrExtreme->SetFillColor(2);
 525     gSigmaCorrExtreme->Draw("3Psame");
 526     csigmaExtreme->SetLogy();
 527     csigmaExtreme->Update();
 528
 529 //     cout << endl << " Sytematics (Extreme approach)" <<endl;
 530 //     for(Int_t item=0; item<gSigmaCorrExtreme->GetN(); item++){
 531 //       Double_t center=0., value=0.;
 532 //       gSigmaCorrExtreme->GetPoint(item,center,value);
 533 //       Double_t highunc = gSigmaCorrExtreme->GetErrorYhigh(item) / value ;
 534 //       Double_t lowunc = gSigmaCorrExtreme->GetErrorYlow(item) / value ;
 535 //       cout << "Sigma syst (extreme) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 536 //     }
 537
 538 //     cout << endl << " Sytematics (Conservative approach)" <<endl;
 539 //     for(Int_t item=0; item<gSigmaCorrConservative->GetN(); item++){
 540 //       Double_t center=0., value=0.;
 541 //       gSigmaCorrConservative->GetPoint(item,center,value);
 542 //       Double_t highunc = gSigmaCorrConservative->GetErrorYhigh(item) / value ;
 543 //       Double_t lowunc = gSigmaCorrConservative->GetErrorYlow(item) / value ;
 544 //       cout << "Sigma syst (conservative) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
 545 //     }
 546
 547  }
 548
 549
 550
 551   //
 552   // Write the histograms to the output file
 553   //
 554   cout << " Saving the results ! " << endl<< endl;
 555
 556   out->cd();
 557   //
 558   hDirectMCpt->Write();        hFeedDownMCpt->Write();
 559   hDirectMCptMax->Write();     hDirectMCptMin->Write();
 560   hFeedDownMCptMax->Write();   hFeedDownMCptMin->Write();
 561   if(hDirectEffpt) hDirectEffpt->Write();        if(hFeedDownEffpt) hFeedDownEffpt->Write();
 562   hRECpt->Write();
 563   //
 564   histoYieldCorr->Write();
 565   histoYieldCorrMax->Write();     histoYieldCorrMin->Write();
 566   histoSigmaCorr->Write();
 567   histoSigmaCorrMax->Write();     histoSigmaCorrMin->Write();
 568
 569   if(asym){
 570     gYieldCorr->Write();
 571     gSigmaCorr->Write();
 572     if(gYieldCorrExtreme) gYieldCorrExtreme->Write();
 573     if(gSigmaCorrExtreme) gSigmaCorrExtreme->Write();
 574     if(gYieldCorrConservative) gYieldCorrConservative->Write();
 575     if(gSigmaCorrConservative) gSigmaCorrConservative->Write();
 576     if(asym && gFcConservative) gFcConservative->Write();
 577   }
 578
 579   if(option==1){
 580     histofc->Write();
 581     histofcMax->Write();     histofcMin->Write();
 582     if(asym && gFcExtreme) gFcExtreme->Write();
 583   }
 584
 585
 586   // Draw the cross-section
 587   //  spectra->DrawSpectrum(gPrediction);
 588
 589   //  out->Close();
 590
 591 }