PWGHF/vertexingHF/macros/DrawFpromptVsRaaElossHypoCombined.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 "TMath.h"
  14 #endif
  15
  16 /******************************************************
  17  *
  18  *  Macro to compute fprompt for pPb and PbPb data
  19  *    using the outputs of the HFPtSpectrumRaa maco
  20  *  Questions and complains to Z.Conesa del Valle
  21  *
  22  *******************************************************/
  23
  24 Bool_t elossFDQuadSum = true;
  25
  26 TGraphAsymmErrors * DrawFpromptVsRaaElossHypo(const char *infile="", //const char *outfile="",
  27                                               Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t CentralHypo = 1.0, Bool_t isYieldUnc=true);
  28
  29
  30 //_______________________________________________________________________
  31 void DrawFpromptVsRaaElossHypoCombined(const char *infileNb="", const char *infilefc="", const char *outfile="", Double_t MinHypo=1./3., Double_t MaxHypo=3.0, Double_t CentralHypo = 1.0, Bool_t isYieldUnc=true)
  32 {
  33
  34   TH2F *hdraw = new TH2F("hdraw","fprompt vs pt; p_{T} (GeV/c); f_{prompt}",36,0.,36.,10.,0.,1);
  35   hdraw->Draw();
  36
  37   cout<<endl<<endl<<" ********* Checking Nb method ********* "<<endl<<endl;
  38   TGraphAsymmErrors *gfPromptNb = DrawFpromptVsRaaElossHypo(infileNb,MinHypo,MaxHypo,CentralHypo,isYieldUnc);
  39   gfPromptNb->SetNameTitle("gfPromptNb","fprompt Nb method, Eloss range considered");
  40   gfPromptNb->SetLineColor(kBlack);
  41   gfPromptNb->SetMarkerColor(kGreen+2);
  42   gfPromptNb->SetMarkerStyle(25);
  43   gfPromptNb->Draw("P");
  44
  45   cout<<endl<<endl<<" ********* Checking fc method ********* "<<endl<<endl;
  46   TGraphAsymmErrors *gfPromptfc = DrawFpromptVsRaaElossHypo(infilefc,MinHypo,MaxHypo,CentralHypo,isYieldUnc);
  47   gfPromptfc->SetNameTitle("gfPromptfc","fprompt fc method, Eloss range considered");
  48   gfPromptfc->SetLineColor(kBlack);
  49   gfPromptfc->SetMarkerColor(kRed);
  50   gfPromptfc->SetMarkerStyle(21);
  51   gfPromptfc->Draw("P");
  52
  53   cout<<endl<<endl<<" ********* Computing the envelope ******** "<<endl<<endl;
  54   TGraphAsymmErrors *gFPromptCombined = new TGraphAsymmErrors(0);
  55   gFPromptCombined->SetNameTitle("gFPromptCombined","fprompt Nb + fc, Eloss range considered");
  56
  57   Int_t nbins=gfPromptNb->GetN();
  58   for(Int_t i=0; i<nbins; i++) {
  59     Double_t pt=0.,value=0.,ptfc=0.,valuefc=0.,exlow=0.,eylow=0.,eyNblow=0.,eyfclow=0.,eyhigh=0.,eyNbhigh=0.,eyfchigh=0.;
  60     gfPromptNb->GetPoint(i,pt,value);
  61     gfPromptfc->GetPoint(i,ptfc,valuefc);
  62     exlow = gfPromptNb->GetErrorXlow(i);
  63     eyNblow = gfPromptNb->GetErrorYlow(i);
  64     eyNbhigh = gfPromptNb->GetErrorYhigh(i);
  65     eyfclow = gfPromptfc->GetErrorYlow(i);
  66     eyfchigh = gfPromptfc->GetErrorYhigh(i);
  67     Double_t uncertainties[5]={ value, value-eyNblow, value+eyNbhigh, valuefc-eyfclow, valuefc+eyfchigh };
  68     eylow = value - TMath::MinElement(5,uncertainties);
  69     eyhigh = TMath::MaxElement(5,uncertainties) - value;
  70     gFPromptCombined->SetPoint(i,pt,value);
  71     gFPromptCombined->SetPointError(i,exlow,exlow,eylow,eyhigh);
  72   }
  73
  74   gFPromptCombined->SetFillStyle(0);
  75   gFPromptCombined->SetLineColor(kBlack);
  76   gFPromptCombined->SetMarkerColor(kBlack);
  77   gFPromptCombined->SetMarkerStyle(20);
  78   gFPromptCombined->Draw("2P");
  79
  80   cout<<endl<<endl;
  81
  82   TFile* fout= new TFile(outfile,"recreate");
  83   gfPromptNb->Write();
  84   gfPromptfc->Write();
  85   gFPromptCombined->Write();
  86   fout->Write();
  87 }
  88
  89 //_______________________________________________________________________
  90 TGraphAsymmErrors * DrawFpromptVsRaaElossHypo(const char *infile,
  91                                               //const char *outfile="",
  92                                               Double_t MinHypo, Double_t MaxHypo, Double_t CentralHypo, Bool_t isYieldUnc)
  93 {
  94
  95   Float_t pt=0.,TAB=0.,sigmaPP=0.,invyieldAB=0.,RaaCharm=0.,RaaBeauty=0., fc=0.;
  96   Float_t yieldFdHigh=0., yieldFdLow=0.,RaaCharmFdHigh=0.,RaaCharmFdLow=0.;
  97
  98   TFile *fRaa = new TFile(infile,"read");
  99   TH1D *hRABvsPt = (TH1D*)fRaa->Get("hRABvsPt");
 100   TNtuple *ntRaa = (TNtuple*)fRaa->Get("ntupleRAB");
 101   ntRaa->SetBranchAddress("pt",&pt);
 102   ntRaa->SetBranchAddress("TAB",&TAB);
 103   ntRaa->SetBranchAddress("sigmaPP",&sigmaPP);
 104   ntRaa->SetBranchAddress("invyieldAB",&invyieldAB);
 105   ntRaa->SetBranchAddress("RABCharm",&RaaCharm);
 106   ntRaa->SetBranchAddress("RABBeauty",&RaaBeauty);
 107   ntRaa->SetBranchAddress("RABCharmFDHigh",&RaaCharmFdHigh);
 108   ntRaa->SetBranchAddress("RABCharmFDLow",&RaaCharmFdLow);
 109   ntRaa->SetBranchAddress("invyieldABFDHigh",&yieldFdHigh);
 110   ntRaa->SetBranchAddress("invyieldABFDLow",&yieldFdLow);
 111   ntRaa->SetBranchAddress("fc",&fc);
 112
 113   // define the binning
 114   Int_t nbins = hRABvsPt->GetNbinsX();
 115
 116   //
 117   //
 118   // Search the central value of the energy loss hypothesis Rb = Rc (bin)
 119   //
 120   Double_t ElossCentral[nbins+1];
 121   Double_t fcV[nbins+1], fcVmax[nbins+1], fcVmin[nbins+1];
 122   Double_t fdMax[nbins+1], fdMin[nbins+1];
 123   for(Int_t i=0; i<=nbins; i++) {
 124     ElossCentral[i]=0.;
 125     fcV[i]=0.; fcVmax[i]=0.; fcVmin[i]=6.;
 126     fdMax[i]=0.; fdMin[i]=6.;
 127   }
 128   //
 129   Int_t netries = ntRaa->GetEntries();
 130   //
 131   for(Int_t ientry=0; ientry<=netries; ientry++){
 132     ntRaa->GetEntry(ientry);
 133     Double_t ElossHypo =  1. / (RaaCharm / RaaBeauty) ;
 134     //
 135     // Find the bin for the central Eloss hypo
 136     //
 137     if( TMath::Abs( ElossHypo - CentralHypo ) < 0.075 ){
 138       Int_t hABbin = hRABvsPt->FindBin( pt );
 139       Double_t DeltaIni = TMath::Abs( ElossCentral[ hABbin ] - CentralHypo );
 140       Double_t DeltaV = TMath::Abs( ElossHypo - CentralHypo );
 141       //      cout << " pt " << ptTot << " ECentral Tot " << ElossCentralTot[ hABbin ] << " Ehypo "<< ElossHypo ;
 142       if ( DeltaV < DeltaIni ) ElossCentral[ hABbin ] = ElossHypo;
 143       //      cout << " final ECentral " << ElossCentralTot[ hABbin ] << endl;
 144     }
 145   }
 146
 147   //
 148   for(Int_t ientry=0; ientry<=netries; ientry++){
 149
 150     ntRaa->GetEntry(ientry);
 151     Double_t ElossHypo =  1. / (RaaCharm / RaaBeauty) ;
 152     Int_t hABbin = hRABvsPt->FindBin( pt );
 153
 154     // Get the central value of Raa
 155     if ( ElossHypo == ElossCentral[ hABbin ] ) {
 156
 157       fcV[hABbin] = fc;
 158
 159       // Now look at FONLL band for the central value
 160       Double_t elems[3]= {0., 0., 0.};
 161       //        fdMax[hABbin]= (RaaCharmFdHigh-RaaCharm)/RaaCharm;
 162       //        fdMin[hABbin]= (RaaCharm-RaaCharmFdLow)/RaaCharm;
 163       if(isYieldUnc) {
 164         elems[0]=invyieldAB;
 165         elems[1]=yieldFdHigh;
 166         elems[2]=yieldFdLow;
 167       } else {
 168         elems[0]=RaaCharm;
 169         elems[1]=RaaCharmFdHigh;
 170         elems[2]=RaaCharmFdLow;
 171       }
 172       fdMax[hABbin]= (TMath::MaxElement(3,elems)-elems[0])/elems[0];
 173       fdMin[hABbin]= (elems[0]-TMath::MinElement(3,elems))/elems[0];
 174       //      cout<<" Check FONLL band giving for pt="<< pt <<" +"<< fdMax[hABbin]<< " -"<< fdMin[hABbin]<<" (relative variation)"<<endl;
 175     }
 176
 177     // Here check the Eloss band
 178     if( RaaCharm>0 && ElossHypo >= MinHypo && ElossHypo <=MaxHypo ) {
 179       if( fc < fcVmin[hABbin] ) fcVmin[hABbin] = fc;
 180       if( fc > fcVmax[hABbin] ) fcVmax[hABbin] = fc;
 181     }
 182   }
 183
 184   TGraphAsymmErrors *gFPrompt = new TGraphAsymmErrors(0);
 185   Int_t j=0;
 186   for(Int_t i=1; i<=nbins; i++) {
 187     Double_t xpt = hRABvsPt->GetBinCenter(i);
 188     Double_t width  = hRABvsPt->GetBinWidth(i)/2;
 189     Double_t elossLow = (fcV[i]-fcVmin[i])/fcV[i];
 190     Double_t elossHigh = (fcVmax[i]-fcV[i])/fcV[i];
 191     Double_t fcLow = fcV[i]-fcVmin[i];
 192     Double_t fcHigh = fcVmax[i]-fcV[i];
 193     //    cout<<" pt="<<xpt<<" , fprompt="<<fcV[i]<<" ["<< fcVmin[i] <<","<<fcVmax[i] <<"] ( Eloss) "<<endl;
 194
 195     if(elossFDQuadSum){ // add quadratically
 196       fcLow = TMath::Sqrt( fdMin[i]*fdMin[i] + elossLow*elossLow )*fcV[i];
 197       fcHigh = TMath::Sqrt( fdMax[i]*fdMax[i] + elossHigh*elossHigh )*fcV[i];
 198     } else { // add linearly
 199       fcLow = ( fdMin[i] + elossLow )*fcV[i];
 200       fcHigh = ( fdMax[i] + elossHigh )*fcV[i];
 201     }
 202
 203     gFPrompt->SetPoint(j,xpt,fcV[i]);
 204     //    gFPrompt->SetPointError(j,width,width,fcV[i]-fcVmin[i],fcVmax[i]-fcV[i]);
 205     gFPrompt->SetPointError(j,width,width,fcLow,fcHigh);
 206     j++;
 207
 208     cout<<" pt="<<xpt<<" , fprompt="<<fcV[i]<<" +"<<fcHigh<<" -"<<fcLow<<" ( splitted it is +"<<elossLow*fcV[i]<<" -"<<elossHigh*fcV[i]<<" Eloss, +"<<fdMax[i]*fcV[i]<<" -"<<fdMin[i]*fcV[i]<<" FONLL-band )"<<endl;
 209   }
 210   //  gFPrompt->Draw("AP");
 211
 212   // if(strcmp(outfile,"")!=0) {
 213   //   TFile *out = new TFile(outfile,"recreate");
 214   //   gFPrompt->Write();
 215   // }
 216
 217   fRaa->Close();
 218
 219   return gFPrompt;
 220
 221 }