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
  16 #include "AliHFPtSpectrum.h"
  17
  18 //
  19 // Macro execution parameters:
  20 //  0) filename with the theoretical predictions  (direct & feed-down)
  21 //  1) acceptance and reconstruction efficiencies file name (direct & feed-down)
  22 //  2) reconstructed spectra file name
  23 //  3) output file name
  24 //  4) Set the feed-down calculation option flag: 0=none, 1=fc only, 2=Nb only
  25 //  5) Set the luminosity
  26 //  6) Set the trigger efficiency
  27 //
  28 void HFPtSpectrum(const char *mcfilename="FeedDownCorrectionMC.root",
  29                   const char *efffilename="Efficiencies.root",
  30                   const char *recofilename="Reconstructed.root",
  31                   const char *outfilename="HFPtSpectrum.root",
  32                   int option=1, double lumi=1.0, double eff_trig=1.0){
  33
  34   //  Set if calculation considers asymmetric uncertainties or not
  35   bool asym = true;
  36
  37   // Set the meson and decay
  38   // (only D0 -> K pi & D+-->K pi pi implemented here)
  39   bool isD0Kpi = true;
  40   bool isDplusKpipi = false;
  41   if (isD0Kpi && isDplusKpipi) {
  42     cout << "Sorry, can not deal with the two corrections at the same time"<<endl;
  43     return;
  44   }
  45
  46   if (option>2) {
  47     cout<< "Bad calculation option, should be <=2"<<endl;
  48     return;
  49   }
  50   if (option==0) asym = false;
  51
  52   //
  53   // Get the histograms from the files
  54   //
  55   TH1D *hDirectMCpt;            // Input MC c-->D spectra
  56   TH1D *hFeedDownMCpt;          // Input MC b-->D spectra
  57   TH1D *hDirectMCpt_max;        // Input MC maximum c-->D spectra
  58   TH1D *hDirectMCpt_min;        // Input MC minimum c-->D spectra
  59   TH1D *hFeedDownMCpt_max;      // Input MC maximum b-->D spectra
  60   TH1D *hFeedDownMCpt_min;      // Input MC minimum b-->D spectra
  61   TH1D *hDirectEffpt;           // c-->D Acceptance and efficiency correction
  62   TH1D *hFeedDownEffpt;         // b-->D Acceptance and efficiency correction
  63   TH1D *hRECpt;                 // all reconstructed D
  64
  65   //
  66   // Define/Get the input histograms
  67   //
  68   TFile * mcfile = new TFile(mcfilename,"read");
  69   if (isD0Kpi){
  70     hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central");
  71     hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central");
  72     hDirectMCpt_max = (TH1D*)mcfile->Get("hD0Kpipred_max");
  73     hDirectMCpt_min = (TH1D*)mcfile->Get("hD0Kpipred_min");
  74     hFeedDownMCpt_max = (TH1D*)mcfile->Get("hD0KpifromBpred_max");
  75     hFeedDownMCpt_min = (TH1D*)mcfile->Get("hD0KpifromBpred_min");
  76   }
  77   else if (isDplusKpipi){
  78     hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central");
  79     hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central");
  80     hDirectMCpt_max = (TH1D*)mcfile->Get("hDpluskpipipred_max");
  81     hDirectMCpt_min = (TH1D*)mcfile->Get("hDpluskpipipred_min");
  82     hFeedDownMCpt_max = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max");
  83     hFeedDownMCpt_min = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min");
  84   }
  85   //
  86   hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra");
  87   hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra");
  88   hDirectMCpt_max->SetNameTitle("hDirectMCpt_max","max direct MC spectra");
  89   hDirectMCpt_min->SetNameTitle("hDirectMCpt_min","min direct MC spectra");
  90   hFeedDownMCpt_max->SetNameTitle("hFeedDownMCpt_max","max feed-down MC spectra");
  91   hFeedDownMCpt_min->SetNameTitle("hFeedDownMCpt_min","min feed-down MC spectra");
  92   //
  93   TFile * efffile = new TFile(efffilename,"read");
  94   hDirectEffpt = (TH1D*)efffile->Get("hDirectEffpt");
  95   hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
  96   hFeedDownEffpt = (TH1D*)efffile->Get("hFeedDownEffpt");
  97   hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff");
  98   //
  99   TFile * recofile = new TFile(recofilename,"read");
 100   hRECpt = (TH1D*)recofile->Get("hRecoAll");
 101   hRECpt->SetNameTitle("hRECpt","Reconstructed spectra");
 102
 103   //
 104   // Define the output histograms
 105   //
 106   TFile *out = new TFile(outfilename,"recreate");
 107   //
 108   //
 109   TH1D *histofc = (TH1D*)hRECpt->Clone();
 110   histofc->SetNameTitle("histofc","fc correction factor");
 111   histofc->Reset();
 112   TH1D *histofc_max = (TH1D*)histofc->Clone();
 113   TH1D *histofc_min = (TH1D*)histofc->Clone();
 114   TH1D *histoYieldCorr = (TH1D*)hRECpt->Clone();
 115   histoYieldCorr->SetNameTitle("histoYieldCorrFc","corrected yield");
 116   histoYieldCorr->Reset();
 117   TH1D *histoYieldCorr_max = (TH1D*)histoYieldCorr->Clone();
 118   histoYieldCorr_max->SetNameTitle("histoYieldCorr_max","max corrected yield");
 119   TH1D *histoYieldCorr_min = (TH1D*)histoYieldCorr->Clone();
 120   histoYieldCorr_min->SetNameTitle("histoYieldCorr_min","min corrected yield");
 121   TH1D *histoSigmaCorr = (TH1D*)hRECpt->Clone();
 122   histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section");
 123   histoSigmaCorr->Reset();
 124   //
 125   int nbins = hRECpt->GetNbinsX();
 126   TGraphAsymmErrors * gFc = new TGraphAsymmErrors(nbins);
 127   TGraphAsymmErrors * gYieldCorr = new TGraphAsymmErrors(nbins);
 128   TGraphAsymmErrors * gSigmaCorr = new TGraphAsymmErrors(nbins);
 129
 130
 131   //
 132   // Main functionalities for the calculation
 133   //
 134
 135   // Define and set the basic option flags
 136   AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option);
 137   spectra->SetFeedDownCalculationOption(option);
 138   spectra->SetComputeAsymmetricUncertainties(asym);
 139
 140   // Feed the input histograms
 141   spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt);
 142   spectra->SetReconstructedSpectrum(hRECpt);
 143   // option specific histos
 144   if(option==1){
 145     spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt);
 146     if(asym) spectra->SetMCptDistributionsBounds(hDirectMCpt_max,hDirectMCpt_min,hFeedDownMCpt_max,hFeedDownMCpt_min);
 147   }
 148   else if(option==2){
 149     spectra->SetFeedDownMCptSpectra(hFeedDownMCpt);
 150     if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCpt_max,hFeedDownMCpt_min);
 151   }
 152
 153   // Set normalization factors (uncertainties set to 0. as example)
 154   spectra->SetLuminosity(lumi,0.);
 155   spectra->SetTriggerEfficiency(eff_trig,0.);
 156
 157   // Do the calculations
 158   cout << " Doing the calculation... "<< endl;
 159   double delta_y = 1.0;
 160   double BR_c=1.0, BR_b_decay=1.0;
 161   spectra->ComputeHFPtSpectrum(delta_y,BR_c,BR_b_decay);
 162   cout << "   ended the calculation, getting the histograms back " << endl;
 163
 164   //
 165   // Get the output histograms
 166   //
 167   // the corrected yield and cross-section
 168   histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum();
 169   histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum();
 170   histoYieldCorr_max = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum();
 171   histoYieldCorr_min = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum();
 172   if (asym) {
 173     gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
 174     gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
 175   }
 176
 177   if (option==0 && asym){
 178     gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)");
 179     gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)");
 180   }
 181   if (option==1){
 182     // fc histos
 183     histofc = (TH1D*)spectra->GetHistoFeedDownCorrection_fc();
 184     histofc_max = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrection_fc();
 185     histofc_min = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrection_fc();
 186     histofc_max->SetNameTitle("hfc_max","max fc correction factor");
 187     histofc_min->SetNameTitle("histofc_min","min fc correction factor");
 188     if (asym) {
 189       gFc = spectra->GetFeedDownCorrection_fc();
 190       gFc->SetNameTitle("gFc","gFc");
 191       gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)");
 192       gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)");
 193     }
 194   }
 195   if (option==2 && asym) {
 196     gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)");
 197     gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)");
 198   }
 199
 200   //
 201   // Now, plot the results ! :)
 202   //
 203
 204   cout << " Drawing the results ! " << endl;
 205
 206   if (option==1) {
 207
 208     TCanvas * cfc = new TCanvas("cfc","Fc");
 209     histofc_max->Draw("c");
 210     histofc->Draw("csame");
 211     histofc_min->Draw("csame");
 212     cfc->Update();
 213
 214     if (asym) {
 215       TH2F *histofc_draw= new TH2F("histofc_draw","histofc (for drawing)",100,0,33.25,100,0.01,1.25);
 216       histofc_draw->SetStats(0);
 217       histofc_draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 218       histofc_draw ->GetXaxis()->SetTitleSize(0.05);
 219       histofc_draw->GetXaxis()->SetTitleOffset(0.95);
 220       histofc_draw->GetYaxis()->SetTitle(" fc ");
 221       histofc_draw->GetYaxis()->SetTitleSize(0.05);
 222       TCanvas *cfc_asym = new TCanvas("cfc_asym","Asymmetric fc (TGraphAsymmErr)");
 223       gFc->SetFillStyle(3001);
 224       gFc->SetLineWidth(3);
 225       gFc->SetMarkerStyle(20);
 226       gFc->SetFillColor(3);
 227       histofc_draw->Draw();
 228       gFc->Draw("3Csame");
 229       gFc->Draw("Xsame");
 230       cfc_asym->Update();
 231     }
 232
 233   }
 234
 235   //
 236   // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra)
 237   //
 238   TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma");
 239   hDirectMCpt->SetMarkerStyle(20);
 240   hDirectMCpt->SetMarkerColor(4);
 241   hDirectMCpt->Draw("p");
 242   histoSigmaCorr->SetMarkerStyle(21);
 243   histoSigmaCorr->SetMarkerColor(2);
 244   histoSigmaCorr->Draw("psame");
 245   histoYieldCorr->SetMarkerStyle(22);
 246   histoYieldCorr->SetMarkerColor(6);
 247   histoYieldCorr->Draw("psame");
 248   hRECpt->SetMarkerStyle(23);
 249   hRECpt->SetMarkerColor(3);
 250   hRECpt->Draw("psame");
 251   cresult->SetLogy();
 252   cresult->Update();
 253
 254   if (asym) {
 255
 256       TH2F *histo_draw = new TH2F("histo_draw","histo (for drawing)",100,0,33.25,100,50.,1e7);
 257       float max = 1.1*gYieldCorr->GetMaximum();
 258       histo_draw->SetAxisRange(0.1,max,"Y");
 259       histo_draw->SetStats(0);
 260       histo_draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 261       histo_draw->GetXaxis()->SetTitleSize(0.05);
 262       histo_draw->GetXaxis()->SetTitleOffset(0.95);
 263       histo_draw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]");
 264       histo_draw->GetYaxis()->SetTitleSize(0.05);
 265       TCanvas * cyield_asym = new TCanvas("cyield_asym","Asymmetric corrected yield (TGraphAsymmErr)");
 266       gYieldCorr->SetFillStyle(3001);
 267       gYieldCorr->SetLineWidth(3);
 268       gYieldCorr->SetMarkerStyle(20);
 269       gYieldCorr->SetFillColor(3);
 270       histo_draw->Draw();
 271       gYieldCorr->Draw("3Csame");
 272       gYieldCorr->Draw("Xsame");
 273       cyield_asym->SetLogy();
 274       cyield_asym->Update();
 275
 276       TH2F *histo2_draw = new TH2F("histo2_draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9);
 277       max = 1.1*gSigmaCorr->GetMaximum();
 278       histo2_draw->SetAxisRange(0.1,max,"Y");
 279       histo2_draw->SetStats(0);
 280       histo2_draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 281       histo2_draw->GetXaxis()->SetTitleSize(0.05);
 282       histo2_draw->GetXaxis()->SetTitleOffset(0.95);
 283       histo2_draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}");
 284       histo2_draw->GetYaxis()->SetTitleSize(0.05);
 285       TCanvas * csigma_asym = new TCanvas("csigma_asym","Asymmetric corrected sigma (TGraphAsymmErr)");
 286       gSigmaCorr->SetFillStyle(3001);
 287       gSigmaCorr->SetLineWidth(3);
 288       gSigmaCorr->SetMarkerStyle(21);
 289       gSigmaCorr->SetFillColor(4);
 290       histo2_draw->Draw();
 291       gSigmaCorr->Draw("3Csame");
 292       gSigmaCorr->Draw("Xsame");
 293       csigma_asym->SetLogy();
 294       csigma_asym->Update();
 295   }
 296
 297
 298
 299   //
 300   // Write the histograms to the output file
 301   //
 302   out->cd();
 303   //
 304   hDirectMCpt->Write();         hFeedDownMCpt->Write();
 305   hDirectMCpt_max->Write();     hDirectMCpt_min->Write();
 306   hFeedDownMCpt_max->Write();   hFeedDownMCpt_min->Write();
 307   hDirectEffpt->Write();        hFeedDownEffpt->Write();
 308   hRECpt->Write();
 309   //
 310   histoYieldCorr->Write();
 311   histoYieldCorr_max->Write();     histoYieldCorr_min->Write();
 312   histoSigmaCorr->Write();
 313
 314   if(asym){
 315     gYieldCorr->Write();
 316     gSigmaCorr->Write();
 317   }
 318
 319   if(option==1){
 320     histofc->Write();
 321     histofc_max->Write();     histofc_min->Write();
 322     if(asym) gFc->Write();
 323   }
 324
 325   // out->Close();
 326
 327 }