Update (Zaida)

diff --git a/PWG3/vertexingHF/AliHFPtSpectrum.cxx b/PWG3/vertexingHF/AliHFPtSpectrum.cxx
index d883112b9b29cfedc21d3369168353e18c8a0a9e..2b983da645095836574d51779cc3caaa96aff240 100644 (file)
--- a/PWG3/vertexingHF/AliHFPtSpectrum.cxx
+++ b/PWG3/vertexingHF/AliHFPtSpectrum.cxx
@@ -33,7 +33,6 @@
 #include "TMath.h"
 #include "TH1.h"
 #include "TH1D.h"
-#include "TCanvas.h"
 #include "TGraphAsymmErrors.h"
 #include "TNamed.h"
 
@@ -169,6 +168,56 @@ AliHFPtSpectrum::~AliHFPtSpectrum(){
 }
   
 
+//_________________________________________________________________________________________________________
+TH1 * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1 *hTheory) {
+  //
+  // Function to rebin the theoretical spectrum 
+  //  with respect to the real-data reconstructed spectrum binning 
+  //
+  
+  if (!hTheory || !fhRECpt) {
+    AliError("Feed-down or reconstructed spectra don't exist");
+    return NULL;
+  }
+
+  //
+  // Get the reconstructed spectra bins & limits
+  Int_t nbins = fhRECpt->GetNbinsX();
+  Int_t nbinsMC = hTheory->GetNbinsX();
+  Double_t *limits = new Double_t[nbins+1];
+  Double_t xlow=0., binwidth=0.;
+  for (Int_t i=0; i<=nbins; i++) {
+    binwidth = fhRECpt->GetBinWidth(i);
+    xlow = fhRECpt->GetBinLowEdge(i);
+    limits[i-1] = xlow;
+  }
+  limits[nbins] = xlow + binwidth;
+
+  //
+  // Define a new histogram with the real-data reconstructed spectrum binning 
+  TH1D * hTheoryRebin = new TH1D("hTheoryRebin"," theoretical rebinned prediction",nbins,limits);
+
+  Double_t sum[nbins], items[nbins];
+  for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
+    
+    for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
+      if( hTheory->GetBinCenter(ibin)>limits[ibinrec] && 
+         hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
+       sum[ibinrec]+=hTheory->GetBinContent(ibin);
+       items[ibinrec]+=1.;
+      }
+    }
+    
+  }
+
+  // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
+  for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
+    hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
+  }
+  
+  return (TH1*)hTheoryRebin;
+}
+
 //_________________________________________________________________________________________________________
 void AliHFPtSpectrum::SetMCptSpectra(TH1 *hDirect, TH1 *hFeedDown){
   //
@@ -188,23 +237,9 @@ void AliHFPtSpectrum::SetMCptSpectra(TH1 *hDirect, TH1 *hFeedDown){
     return;
   }
 
-  //
-  // Get the reconstructed spectra bins & limits
-  Int_t nbins = fhRECpt->GetNbinsX();
-  Double_t *limits = new Double_t[nbins+1];
-  Double_t xlow=0., binwidth=0.;
-  for (Int_t i=0; i<=nbins; i++) {
-    binwidth = fhRECpt->GetBinWidth(i);
-    xlow = fhRECpt->GetBinLowEdge(i);
-    limits[i-1] = xlow;
-    //    cout <<" bin = " << i <<", bin width = "<< binwidth << ", low edge = " << xlow << ", high edge = " << xlow+binwidth <<endl;
-  }
-  limits[nbins] = xlow + binwidth;
-
   //
   // If the predictions shape do not have the same
   //  bin width (check via number of bins) as the reconstructed spectra, change them 
-  Int_t nbinsMC = hDirect->GetNbinsX();
   if (hDirect->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
 
     fhDirectMCpt = hDirect;
@@ -213,51 +248,11 @@ void AliHFPtSpectrum::SetMCptSpectra(TH1 *hDirect, TH1 *hFeedDown){
   } 
   else {
 
-    fhDirectMCpt = new TH1D("fhDirectMCpt"," direct theoretical prediction",nbins,limits);
-    fhFeedDownMCpt = new TH1D("fhFeedDownMCpt"," feed-down theoretical prediction",nbins,limits);
-
-    Double_t sumd[nbins], sumf[nbins], itemsd[nbins], itemsf[nbins];
-    for (Int_t ibinrec=0; ibinrec<=nbins; ibinrec++){
-      sumd[ibinrec]=0.; sumf[ibinrec]=0.; itemsd[ibinrec]=0.; itemsf[ibinrec]=0.;
-    }
-    for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
-
-      for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
-       if( hDirect->GetBinCenter(ibin)>limits[ibinrec] && 
-           hDirect->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumd[ibinrec]+=hDirect->GetBinContent(ibin);
-         itemsd[ibinrec]+=1.;
-         //      cout << " binrec="<< ibinrec << " sumd="<< sumd[ibinrec] << endl;
-       }
-       if( hFeedDown->GetBinCenter(ibin)>limits[ibinrec] && 
-           hFeedDown->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumf[ibinrec]+=hFeedDown->GetBinContent(ibin);
-         itemsf[ibinrec]+=1.;
-         //      cout << " binrec="<< ibinrec << ", sumf=" << sumf[ibinrec] << endl;
-       }
-      }
-
-    }
-
-    // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
-    for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
-      fhDirectMCpt->SetBinContent(ibinrec+1,sumd[ibinrec]/itemsd[ibinrec]);
-      fhFeedDownMCpt->SetBinContent(ibinrec+1,sumf[ibinrec]/itemsf[ibinrec]);
-      //      cout << " Setting: binrec="<< ibinrec <<", at x="<<fhDirectMCpt->GetBinCenter(ibinrec)<< " sumd="<< sumd[ibinrec] << ", sumf=" << sumf[ibinrec] << endl;
-    }
+    fhDirectMCpt = RebinTheoreticalSpectra(hDirect);
+    fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
+    fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown);
+    fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
 
-    // control plot
-    TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
-    cTheoryRebin->Divide(1,2);
-    cTheoryRebin->cd(1);
-    hDirect->Draw("");
-    fhDirectMCpt->SetLineColor(2);
-    fhDirectMCpt->Draw("same");
-    cTheoryRebin->cd(2);
-    hFeedDown->Draw("");
-    fhFeedDownMCpt->SetLineColor(2);
-    fhFeedDownMCpt->Draw("same");
-    cTheoryRebin->Update();
   }
 
 }
@@ -274,22 +269,9 @@ void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1 *hFeedDown){
     return;
   }
 
-  //
-  // Get the reconstructed spectra bins & limits
-  Int_t nbins = fhRECpt->GetNbinsX();
-  Double_t *limits = new Double_t[nbins+1];
-  Double_t xlow=0., binwidth=0.;
-  for (Int_t i=0; i<=nbins; i++) {
-    binwidth = fhRECpt->GetBinWidth(i);
-    xlow = fhRECpt->GetBinLowEdge(i);
-    limits[i-1] = xlow;
-  }
-  limits[nbins] = xlow + binwidth;
-
   //
   // If the predictions shape do not have the same
   //  bin width (check via number of bins) as the reconstructed spectra, change them 
-  Int_t nbinsMC = hFeedDown->GetNbinsX();
   if (hFeedDown->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
 
     fhFeedDownMCpt = hFeedDown;
@@ -297,32 +279,9 @@ void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1 *hFeedDown){
   } 
   else {
 
-    fhFeedDownMCpt = new TH1D("fhFeedDownMCpt"," feed-down theoretical prediction",nbins,limits);
-
-    Double_t sumf[nbins], itemsf[nbins];
-    for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
-
-      for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
-       if( hFeedDown->GetBinCenter(ibin)>limits[ibinrec] && 
-           hFeedDown->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumf[ibinrec]+=hFeedDown->GetBinContent(ibin);
-         itemsf[ibinrec]+=1.;
-       }
-      }
+    fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown);
+    fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
 
-    }
-
-    // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
-    for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
-      fhFeedDownMCpt->SetBinContent(ibinrec+1,sumf[ibinrec]/itemsf[ibinrec]);
-    }
-
-    // control plot
-    TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
-    hFeedDown->Draw("");
-    fhFeedDownMCpt->SetLineColor(2);
-    fhFeedDownMCpt->Draw("same");
-    cTheoryRebin->Update();
   }
 
 }
@@ -349,22 +308,10 @@ void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1 *hDirectMax, TH1 *hDirectMi
     return;
   }
 
-  //
-  // Get the reconstructed spectra bins & limits
-  Int_t nbins = fhRECpt->GetNbinsX();
-  Double_t *limits = new Double_t[nbins+1];
-  Double_t xlow=0., binwidth=0.;
-  for (Int_t i=0; i<=nbins; i++) {
-    binwidth = fhRECpt->GetBinWidth(i);
-    xlow = fhRECpt->GetBinLowEdge(i);
-    limits[i-1] = xlow;
-  }
-  limits[nbins] = xlow + binwidth;
 
   //
   // If the predictions shape do not have the same
   //  bin width (check via number of bins) as the reconstructed spectra, change them 
-  Int_t nbinsMC =  hDirectMax->GetNbinsX();
   if (hFeedDownMax->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
     
     fhDirectMCptMax = hDirectMax;
@@ -375,66 +322,15 @@ void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1 *hDirectMax, TH1 *hDirectMi
   } 
   else {
 
-    fhDirectMCptMax = new TH1D("fhDirectMCptMax"," maximum direct theoretical prediction",nbins,limits);
-    fhDirectMCptMin = new TH1D("fhDirectMCptMin"," minimum direct theoretical prediction",nbins,limits);
-    fhFeedDownMCptMax = new TH1D("fhFeedDownMCptMax"," maximum feed-down theoretical prediction",nbins,limits);
-    fhFeedDownMCptMin = new TH1D("fhFeedDownMCptMin"," minimum feed-down theoretical prediction",nbins,limits);
+    fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax);
+    fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
+    fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin);
+    fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
+    fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax);
+    fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
+    fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin);
+    fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
 
-    Double_t sumdmax[nbins], sumfmax[nbins], itemsdmax[nbins], itemsfmax[nbins];
-    Double_t sumdmin[nbins], sumfmin[nbins], itemsdmin[nbins], itemsfmin[nbins];
-    for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
-
-      for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
-       if( hDirectMax->GetBinCenter(ibin)>limits[ibinrec] && 
-           hDirectMax->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumdmax[ibinrec]+=hDirectMax->GetBinContent(ibin);
-         itemsdmax[ibinrec]+=1.;
-       }
-       if( hDirectMin->GetBinCenter(ibin)>limits[ibinrec] && 
-           hDirectMin->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumdmin[ibinrec]+=hDirectMin->GetBinContent(ibin);
-         itemsdmin[ibinrec]+=1.;
-       }
-       if( hFeedDownMax->GetBinCenter(ibin)>limits[ibinrec] && 
-           hFeedDownMax->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumfmax[ibinrec]+=hFeedDownMax->GetBinContent(ibin);
-         itemsfmax[ibinrec]+=1.;
-       }
-       if( hFeedDownMin->GetBinCenter(ibin)>limits[ibinrec] && 
-           hFeedDownMin->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumfmin[ibinrec]+=hFeedDownMin->GetBinContent(ibin);
-         itemsfmin[ibinrec]+=1.;
-       }
-      }
-
-    }
-
-    // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
-    for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
-      fhDirectMCptMax->SetBinContent(ibinrec+1,sumdmax[ibinrec]/itemsdmax[ibinrec]);
-      fhDirectMCptMin->SetBinContent(ibinrec+1,sumdmin[ibinrec]/itemsdmin[ibinrec]);
-      fhFeedDownMCptMax->SetBinContent(ibinrec+1,sumfmax[ibinrec]/itemsfmax[ibinrec]);
-      fhFeedDownMCptMin->SetBinContent(ibinrec+1,sumfmin[ibinrec]/itemsfmin[ibinrec]);
-    }
-
-    // control plot
-    TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
-    cTheoryRebinLimits->Divide(1,2);
-    cTheoryRebinLimits->cd(1);
-    hDirectMax->Draw("");
-    fhDirectMCptMax->SetLineColor(2);
-    fhDirectMCptMax->Draw("same");
-    hDirectMin->Draw("same");
-    fhDirectMCptMin->SetLineColor(2);
-    fhDirectMCptMin->Draw("same");
-    cTheoryRebinLimits->cd(2);
-    hFeedDownMax->Draw("");
-    fhFeedDownMCptMax->SetLineColor(2);
-    fhFeedDownMCptMax->Draw("same");
-    hFeedDownMin->Draw("same");
-    fhFeedDownMCptMin->SetLineColor(2);
-    fhFeedDownMCptMin->Draw("same");
-    cTheoryRebinLimits->Update();
   }
 
 }
@@ -459,22 +355,10 @@ void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1 *hFeedDownMax, TH1
     return;
   }
 
-  //
-  // Get the reconstructed spectra bins & limits
-  Int_t nbins = fhRECpt->GetNbinsX();
-  Double_t *limits = new Double_t[nbins+1];
-  Double_t xlow=0., binwidth=0.;
-  for (Int_t i=0; i<=nbins; i++) {
-    binwidth = fhRECpt->GetBinWidth(i);
-    xlow = fhRECpt->GetBinLowEdge(i);
-    limits[i-1] = xlow;
-  }
-  limits[nbins] = xlow + binwidth;
 
   //
   // If the predictions shape do not have the same
   //  bin width (check via number of bins) as the reconstructed spectra, change them 
-  Int_t nbinsMC = hFeedDownMax->GetNbinsX();
   if (hFeedDownMax->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
     
     fhFeedDownMCptMax = hFeedDownMax;
@@ -483,43 +367,11 @@ void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1 *hFeedDownMax, TH1
   } 
   else {
 
-    fhFeedDownMCptMax = new TH1D("fhFeedDownMCptMax"," maximum feed-down theoretical prediction",nbins,limits);
-    fhFeedDownMCptMin = new TH1D("fhFeedDownMCptMin"," minimum feed-down theoretical prediction",nbins,limits);
+    fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax);
+    fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
+    fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin);
+    fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
 
-    Double_t sumfmax[nbins], itemsfmax[nbins];
-    Double_t sumfmin[nbins], itemsfmin[nbins];
-    for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
-
-      for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
-       if( hFeedDownMax->GetBinCenter(ibin)>limits[ibinrec] && 
-           hFeedDownMax->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumfmax[ibinrec]+=hFeedDownMax->GetBinContent(ibin);
-         itemsfmax[ibinrec]+=1.;
-       }
-       if( hFeedDownMin->GetBinCenter(ibin)>limits[ibinrec] && 
-           hFeedDownMin->GetBinCenter(ibin)<limits[ibinrec+1]){
-         sumfmin[ibinrec]+=hFeedDownMin->GetBinContent(ibin);
-         itemsfmin[ibinrec]+=1.;
-       }
-      }
-
-    }
-
-    // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
-    for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
-      fhFeedDownMCptMax->SetBinContent(ibinrec+1,sumfmax[ibinrec]/itemsfmax[ibinrec]);
-      fhFeedDownMCptMin->SetBinContent(ibinrec+1,sumfmin[ibinrec]/itemsfmin[ibinrec]);
-    }
-
-    // control plot
-    TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
-    hFeedDownMax->Draw("");
-    fhFeedDownMCptMax->SetLineColor(2);
-    fhFeedDownMCptMax->Draw("same");
-    hFeedDownMin->Draw("same");
-    fhFeedDownMCptMin->SetLineColor(2);
-    fhFeedDownMCptMin->Draw("same");
-    cTheoryRebinLimits->Update();
   }
 
 }
@@ -703,9 +555,9 @@ void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRat
 }
 
 //_________________________________________________________________________________________________________
-TH1 * AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco) {
+TH1 * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco) {
   //
-  // Function that computes the Direct  acceptance and efficiency correction
+  // Function that computes the acceptance and efficiency correction
   //  based on the simulated and reconstructed spectra
   //  and using the reconstructed spectra bin width
   //
@@ -727,7 +579,7 @@ TH1 * AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1 *hSimu, TH1 *hR
   }
   limits[nbins] = xlow + binwidth;
 
-  fhDirectEffpt = new TH1D("fhDirectEffpt"," direct acceptance #times efficiency",nbins,limits);
+  TH1D * hEfficiency = new TH1D("hEfficiency"," acceptance #times efficiency",nbins,limits);
   
   Double_t sumSimu[nbins], sumReco[nbins];
   for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
@@ -751,66 +603,52 @@ TH1 * AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1 *hSimu, TH1 *hR
   for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
     eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
     erreff = TMath::Sqrt( eff * (1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
-    fhDirectEffpt->SetBinContent(ibinrec+1,eff);
-    fhDirectEffpt->SetBinError(ibinrec+1,erreff);
+    hEfficiency->SetBinContent(ibinrec+1,eff);
+    hEfficiency->SetBinError(ibinrec+1,erreff);
   }
 
-  return (TH1*)fhDirectEffpt;
+  return (TH1*)hEfficiency;
 }
 
 //_________________________________________________________________________________________________________
-TH1 * AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco) {
+TH1 * AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco) {
   //
-  // Function that computes the Feed-Down acceptance and efficiency correction
+  // Function that computes the Direct  acceptance and efficiency correction
   //  based on the simulated and reconstructed spectra
   //  and using the reconstructed spectra bin width
   //
   //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
   // 
-  
+
   if(!fhRECpt){
     AliInfo("Hey, the reconstructed histogram was not set yet !"); 
     return NULL;
   }
 
-  Int_t nbins = fhRECpt->GetNbinsX();
-  Double_t *limits = new Double_t[nbins+1];
-  Double_t xlow=0.,binwidth=0.;
-  for (Int_t i=0; i<=nbins; i++) {
-    binwidth = fhRECpt->GetBinWidth(i);
-    xlow = fhRECpt->GetBinLowEdge(i);
-    limits[i-1] = xlow;
-  }
-  limits[nbins] = xlow + binwidth;
+  fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco);
+  fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
 
-  fhFeedDownEffpt = new TH1D("fhFeedDownEffpt"," feed-down acceptance #times efficiency",nbins,limits);
-  
-  Double_t sumSimu[nbins], sumReco[nbins];
-  for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
+  return (TH1*)fhDirectEffpt;
+}
 
-    for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
-      if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] && 
-          hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
-       sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
-      }
-      if ( hReco->GetBinCenter(ibin)>limits[ibinrec] && 
-          hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
-       sumReco[ibinrec]+=hReco->GetBinContent(ibin);
-      }
-    }
-    
-  }
+//_________________________________________________________________________________________________________
+TH1 * AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco) {
+  //
+  // Function that computes the Feed-Down acceptance and efficiency correction
+  //  based on the simulated and reconstructed spectra
+  //  and using the reconstructed spectra bin width
+  //
+  //  eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
+  // 
   
-  // the efficiency is computed as reco/sim (in each bin)
-  //  its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
-  Double_t eff=0., erreff=0.;
-  for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
-    eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
-    erreff = TMath::Sqrt( eff * (1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
-    fhFeedDownEffpt->SetBinContent(ibinrec+1,eff);
-    fhFeedDownEffpt->SetBinError(ibinrec+1,erreff);
+  if(!fhRECpt){
+    AliInfo("Hey, the reconstructed histogram was not set yet !"); 
+    return NULL;
   }
 
+  fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco);
+  fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
+
   return (TH1*)fhFeedDownEffpt;
 }
 

diff --git a/PWG3/vertexingHF/AliHFPtSpectrum.h b/PWG3/vertexingHF/AliHFPtSpectrum.h
index 485f10797baad4e0f10cba2b6b863afe6caa2142..da85742c2a505f607b3feffe2368db78b8468d90 100644 (file)
--- a/PWG3/vertexingHF/AliHFPtSpectrum.h
+++ b/PWG3/vertexingHF/AliHFPtSpectrum.h
@@ -70,6 +70,16 @@ class AliHFPtSpectrum: public TNamed
   //
   // Getters
   //
+  // Return the theoretical predictions used for the calculation (rebinned if needed)
+  TH1 * GetDirectTheoreticalSpectrum() { return (fhDirectMCpt ? fhDirectMCpt : NULL); }
+  TH1 * GetDirectTheoreticalUpperLimitSpectrum() { return (fhDirectMCptMax ? fhDirectMCptMax : NULL); }
+  TH1 * GetDirectTheoreticalLowerLimitSpectrum() { return (fhDirectMCptMin ? fhDirectMCptMin : NULL); }
+  TH1 * GetFeedDownTheoreticalSpectrum() { return (fhFeedDownMCpt ? fhFeedDownMCpt : NULL); }
+  TH1 * GetFeedDownTheoreticalUpperLimitSpectrum() { return (fhFeedDownMCptMax ? fhFeedDownMCptMax : NULL); }
+  TH1 * GetFeedDownTheoreticalLowerLimitSpectrum() { return (fhFeedDownMCptMin ? fhFeedDownMCptMin : NULL); }
+  // Return the acceptance and efficiency corrections (rebinned if needed)
+  TH1 * GetDirectAccEffCorrection() { return (fhDirectEffpt ? fhDirectEffpt : NULL); }
+  TH1 * GetFeedDownAccEffCorrection() { return (fhFeedDownEffpt ? fhFeedDownEffpt : NULL); }
   // Return the TGraphAsymmErrors of the feed-down correction
   TGraphAsymmErrors * GetFeedDownCorrectionFc() { return (fgFc ?  fgFc : NULL); }
   // Return the histogram of the feed-down correction
@@ -128,6 +138,10 @@ class AliHFPtSpectrum: public TNamed
 
   // Check histograms consistency function
   Bool_t CheckHistosConsistency(TH1 *h1, TH1 *h2);
+  // Function to rebin the theoretical spectra in the data-reconstructed spectra binning
+  TH1 * RebinTheoreticalSpectra(TH1 *hTheory);
+  // Function to estimate the efficiency in the data-reconstructed spectra binning
+  TH1 * EstimateEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco);
 
   //
   // Input spectra

diff --git a/PWG3/vertexingHF/macros/ComputeEfficiencyInputFromTwoSteps.C b/PWG3/vertexingHF/macros/ComputeEfficiencyInputFromTwoSteps.C
new file mode 100644 (file)
index 0000000..8e17131
--- /dev/null
+++ b/PWG3/vertexingHF/macros/ComputeEfficiencyInputFromTwoSteps.C
@@ -0,0 +1,86 @@
+#include "TFile.h"
+#include "TH1.h"
+#include "TH1D.h"
+#include "TMath.h"
+
+#include <Riostream.h>
+
+//
+// Multiply two histos of variable bin width
+//
+TH1 * MultiplyHistos(TH1 *h1, TH1 *h2) {
+  
+  Int_t nbins1 = h1->GetNbinsX();
+  Int_t nbins2 = h2->GetNbinsX();
+  Double_t binwidth1 = h1->GetBinWidth(1);
+  Double_t binwidth2 = h2->GetBinWidth(1);
+
+  if ( nbins1!=nbins2 || binwidth1!=binwidth2 ) {
+    cout << "Histos do not have the same binning, do not seem compatible" << endl;
+    return NULL;
+  }
+
+  //
+  // Get the bins & limits
+  Double_t *limits = new Double_t[nbins1+1];
+  Double_t xlow=0., binwidth=0.;
+  for (Int_t i=0; i<=nbins1; i++) {
+    binwidth = h1->GetBinWidth(i);
+    xlow = h1->GetBinLowEdge(i);
+    limits[i-1] = xlow;
+  }
+  limits[nbins1] = xlow + binwidth;
+
+  TH1D *hMultiply = new TH1D("hMultiply","hMultiply",nbins1,limits);
+
+  Double_t value=0., err=0.;
+  for (Int_t ibin=0; ibin<nbins1; ibin++) {
+    value = h1->GetBinContent(ibin) * h2->GetBinContent(ibin);
+    err = value * TMath::Sqrt(  (h1->GetBinError(ibin)/h1->GetBinContent(ibin)) * (h1->GetBinError(ibin)/h1->GetBinContent(ibin))  +
+                               (h2->GetBinError(ibin)/h2->GetBinContent(ibin)) * (h2->GetBinError(ibin)/h2->GetBinContent(ibin))   );
+    hMultiply->SetBinContent(ibin,value);
+    hMultiply->SetBinError(ibin,err);
+  }
+  
+  return (TH1*)hMultiply;
+}
+
+//
+// Main function
+// 
+void ComputeEfficiencyInputFromTwoSteps (const char* recolhc10d3filename="Distributions.root",
+                                        const char* recolhc10d3histoname="RECPIDpt",
+                                        const char* simuAcclhc10d3filename="Distributions.root",
+                                        const char* simuAcclhc10d3histoname="MCAccpt",
+                                        const char* simuLimAcclhc10d4filename="",
+                                        const char* simuLimAcclhc10d4histoname="MCLimAccpt",
+                                        const char* simuAcclhc10d4filename="",
+                                        const char* simuAcclhc10d4histoname="MCAccpt",
+                                        const char* outfilename="ComputeEfficiencyInputFromTwoSteps.root") 
+{
+
+  TFile *recolhc10d3file = new TFile(recolhc10d3filename,"read");
+  TH1D *hrecolhc10d3 = (TH1D*)recolhc10d3file->Get(recolhc10d3histoname);
+
+  TFile *simuAcclhc10d3file = new TFile(simuAcclhc10d3filename,"read");
+  TH1D *hsimuAcclhc10d3 = (TH1D*)simuAcclhc10d3file->Get(simuAcclhc10d3histoname);
+
+  TFile *simuLimAcclhc10d4file = new TFile(simuLimAcclhc10d4filename,"read");
+  TH1D *hsimuLimAcclhc10d4 = (TH1D*)simuLimAcclhc10d4file->Get(simuLimAcclhc10d4histoname);
+
+  TFile *simuAcclhc10d4file = new TFile(simuAcclhc10d4filename,"read");
+  TH1D *hsimuAcclhc10d4 = (TH1D*)simuAcclhc10d4file->Get(simuAcclhc10d4histoname);
+
+
+  TFile *out = new TFile(outfilename,"recreate");
+  TH1D *hRecoPIDCorr = (TH1D*)MultiplyHistos(hrecolhc10d3,hsimuAcclhc10d4);
+  hRecoPIDCorr->SetNameTitle("hRecoPIDCorr","hRecoPIDCorr");
+  TH1D *hSimuCorr = (TH1D*)MultiplyHistos(hsimuAcclhc10d3,hsimuLimAcclhc10d4);
+  hSimuCorr->SetNameTitle("hSimuCorr","hSimuCorr");
+
+  out->cd(); 
+  hRecoPIDCorr->Write();
+  hSimuCorr->Write();
+  out->Close();
+
+}

diff --git a/PWG3/vertexingHF/macros/HFPtSpectrum.C b/PWG3/vertexingHF/macros/HFPtSpectrum.C
index 37639d0893c57f781d962736aad578ac15987b27..cd6846761e9cf45d8203fa5f3d6fe44daa0d53cf 100644 (file)
--- a/PWG3/vertexingHF/macros/HFPtSpectrum.C
+++ b/PWG3/vertexingHF/macros/HFPtSpectrum.C
@@ -28,7 +28,7 @@
 //
 void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
                    const char *efffilename="Efficiencies.root",
-                   const char *recofilename="Reconstructed.root",
+                   const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0",
                    const char *outfilename="HFPtSpectrum.root",
                    int option=1, double lumi=1.0, double effTrig=1.0,
                    const char *directsimufilename="", const char *directsimuhistoname="CFHFccontainer0_New_3Prong_SelStep0_proj-pt", 
@@ -135,7 +135,7 @@ void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
   //
   //
   TFile * recofile = new TFile(recofilename,"read");
-  hRECpt = (TH1D*)recofile->Get("hRecoAll");
+  hRECpt = (TH1D*)recofile->Get(recohistoname);
   hRECpt->SetNameTitle("hRECpt","Reconstructed spectra");
 
   //
@@ -268,6 +268,49 @@ void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
 
   cout << " Drawing the results ! " << endl;
 
+  // control plots
+  TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
+  cTheoryRebin->Divide(1,2);
+  cTheoryRebin->cd(1);
+  hDirectMCpt->Draw("");
+  TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum();
+  hDirectMCptRebin->SetNameTitle("hDirectMCptRebin","hDirectMCptRebin");
+  hDirectMCptRebin->SetLineColor(2);
+  hDirectMCptRebin->Draw("same");
+  cTheoryRebin->cd(2);
+  hFeedDownMCpt->Draw("");
+  TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum();
+  hFeedDownRebin->SetNameTitle("hFeedDownRebin","hFeedDownRebin");
+  hFeedDownRebin->SetLineColor(2);
+  hFeedDownRebin->Draw("same");
+  cTheoryRebin->Update();
+
+  TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
+  cTheoryRebinLimits->Divide(1,2);
+  cTheoryRebinLimits->cd(1);
+  hDirectMCptMax->Draw("");
+  TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum();
+  hDirectMCptMaxRebin->SetNameTitle("hDirectMCptMaxRebin","hDirectMCptMaxRebin");
+  hDirectMCptMaxRebin->SetLineColor(2);
+  hDirectMCptMaxRebin->Draw("same");
+  hDirectMCptMin->Draw("same");
+  TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum();
+  hDirectMCptMinRebin->SetNameTitle("hDirectMCptMinRebin","hDirectMCptMinRebin");
+  hDirectMCptMinRebin->SetLineColor(2);
+  hDirectMCptMinRebin->Draw("same");
+  cTheoryRebinLimits->cd(2);
+  hFeedDownMCptMax->Draw("");
+  TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum();
+  hFeedDownMaxRebin->SetNameTitle("hFeedDownMaxRebin","hFeedDownMaxRebin");
+  hFeedDownMaxRebin->SetLineColor(2);
+  hFeedDownMaxRebin->Draw("same");
+  hFeedDownMCptMin->Draw("same");
+  TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum();
+  hFeedDownMinRebin->SetNameTitle("hFeedDownMinRebin","hFeedDownMinRebin");
+  hFeedDownMinRebin->SetLineColor(2);
+  hFeedDownMinRebin->Draw("same");
+  cTheoryRebinLimits->Update();
+  
   if (option==1) {
 
     TCanvas * cfc = new TCanvas("cfc","Fc");
@@ -378,6 +421,8 @@ void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
   //
   // Write the histograms to the output file
   //
+  cout << " Saving the results ! " << endl<< endl;
+
   out->cd();
   //
   hDirectMCpt->Write();         hFeedDownMCpt->Write();
@@ -402,6 +447,6 @@ void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
     if(asym) gFc->Write();
   }
 
-  // out->Close();
+  //  out->Close();
 
 }