Use TF1 or TH1 as weight for beauty (Diego)

[u/mrichter/AliRoot.git] / PWG / muon / AliAnalysisTaskSingleMu.cxx

diff --git a/PWG/muon/AliAnalysisTaskSingleMu.cxx b/PWG/muon/AliAnalysisTaskSingleMu.cxx
index e2c028500ff732e2fb0e5d3ac0c6e00eb316136d..74b0964b7d0c8a776502bda5023df23281a7cec6 100644 (file)
--- a/PWG/muon/AliAnalysisTaskSingleMu.cxx
+++ b/PWG/muon/AliAnalysisTaskSingleMu.cxx
@@ -43,6 +43,7 @@
 #include "TStyle.h"
 //#include "TMCProcess.h"
 #include "TArrayI.h"
+#include "TArrayD.h"
 #include "TPaveStats.h"
 #include "TFitResultPtr.h"
 #include "TFile.h"
@@ -225,8 +226,6 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
     ((TH1*)GetMergeableObject(physSel, trigClassName, centrality, "hIpVtx"))->Fill(ipVz);
   }
   
-  TF1* beautyMuWgt = static_cast<TF1*>(GetWeight("beautyMu"));
-
   // Bool_t isPileupFromSPD = ( fAODEvent && ! fAODEvent->GetTracklets() ) ? InputEvent()->IsPileupFromSPD(3, 0.8, 3., 2., 5.) : InputEvent()->IsPileupFromSPDInMultBins(); // Avoid break when reading Muon AODs (tracklet info is not present and IsPileupFromSPDInMultBins crashes // UNCOMMENT TO REJECT PILEUP
   // if ( isPileupFromSPD ) return; // UNCOMMENT TO REJECT PILEUP
   
@@ -239,7 +238,7 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
     
     TObjArray selectedTracks(nTracks);
     TArrayI trackSources(nTracks);
-    TArrayF trackWgt(nTracks);
+    TArrayD trackWgt(nTracks);
     trackWgt.Reset(1.);
     for (Int_t itrack = 0; itrack < nTracks; itrack++) {
       track = ( istep == kStepReconstructed ) ? AliAnalysisMuonUtility::GetTrack(itrack,InputEvent()) : MCEvent()->GetTrack(itrack);
@@ -255,9 +254,16 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
       selectedTracks.AddAt(track,itrack);
       trackSources[itrack] = GetParticleType(track);
       
-      if ( trackSources[itrack] == kBeautyMu && beautyMuWgt ) {
-        AliVParticle* mcTrack = ( kStepReconstructed ) ? MCEvent()->GetTrack(track->GetLabel()) : track;
-        trackWgt[itrack] = beautyMuWgt->Eval(mcTrack->Pt());
+      TObject* wgtObj = GetWeight(fSrcKeys->At(trackSources[itrack])->GetName());
+      
+      if ( wgtObj  ) {
+        AliVParticle* mcTrack = ( istep == kStepReconstructed ) ? MCEvent()->GetTrack(track->GetLabel()) : track;
+        if ( wgtObj->IsA() == TF1::Class() ) trackWgt[itrack] = static_cast<TF1*>(wgtObj)->Eval(mcTrack->Pt());
+        else if ( wgtObj->IsA()->InheritsFrom(TH1::Class()) ) {
+          TH1* wgtHisto = static_cast<TH1*>(wgtObj);
+          trackWgt[itrack] = wgtHisto->GetBinContent(wgtHisto->GetXaxis()->FindBin(mcTrack->Pt()));
+        }
+        AliDebug(3,Form("Apply weights %s:  pt %g  gen pt %g  weight %g",wgtObj->GetName(),track->Pt(),mcTrack->Pt(),trackWgt[itrack]));
 //        Int_t iancestor = fUtilityMuonAncestor->GetAncestor(track,MCEvent());
 //        AliVParticle* motherPart = MCEvent()->GetTrack(iancestor);
 //        trackWgt[itrack] = beautyMuWgt->GetBinContent(beautyMuWgt->GetXaxis()->FindBin(motherPart->Pt()));
@@ -422,11 +428,8 @@ void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
   Bool_t isMC = furtherOpt.Contains("MC");
   
   TAxis* srcAxis = gridSparseArray[0]->GetAxis(kHvarMotherType);
-  Int_t unIdBin = srcAxis->GetNbins();
-  for ( Int_t ibin=1; ibin<=srcAxis->GetNbins(); ibin++ ) {
-    currName = srcAxis->GetBinLabel(ibin);
-    if ( currName.Contains("Unidentified") ) unIdBin = ibin;
-  }
+  Int_t unIdBin = srcAxis->FindBin(fSrcKeys->At(kUnidentified)->GetName());
+  if ( unIdBin < 1 ) unIdBin = srcAxis->GetNbins();
   
   Int_t firstSrcBin = ( isMC ) ? 1 : unIdBin;
   Int_t lastSrcBin  = ( isMC ) ? srcAxis->GetNbins() - 1 : unIdBin;
@@ -618,14 +621,15 @@ void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
   sumMothers[kRecoHF].Set(0);
   sumMothers[kRecoBkg].Set(0);
   sumMothers[kInputHF].Set(3);
-  sumMothers[kInputHF][0] = kCharmMu;
-  sumMothers[kInputHF][1] = kBeautyMu;
-  sumMothers[kInputHF][2] = kQuarkoniumMu;
+  
+  sumMothers[kInputHF][0] = srcAxis->FindBin(fSrcKeys->At(kCharmMu)->GetName());
+  sumMothers[kInputHF][1] = srcAxis->FindBin(fSrcKeys->At(kBeautyMu)->GetName());
+  sumMothers[kInputHF][2] = srcAxis->FindBin(fSrcKeys->At(kQuarkoniumMu)->GetName());
   sumMothers[kInputDecay].Set(1);
-  sumMothers[kInputDecay][0] = kDecayMu;
+  sumMothers[kInputDecay][0] = srcAxis->FindBin(fSrcKeys->At(kDecayMu)->GetName());
   sumMothers[kRecoAll].Set(srcAxis->GetNbins());
-  for ( Int_t isrc=1; isrc<srcAxis->GetNbins(); ++isrc ) {
-    sumMothers[kRecoAll][isrc-1] = isrc;
+  for ( Int_t isrc=0; isrc<srcAxis->GetNbins(); ++isrc ) {
+    sumMothers[kRecoAll][isrc] = isrc+1;
   }
   
   meanZ = vtxFit->GetParameter(1);
@@ -683,6 +687,7 @@ void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
       histo->Divide(eventVertex);
       Double_t norm = histo->GetBinContent(histo->FindBin(0.));
       histo->GetYaxis()->SetTitle("#frac{dN_{#mu}}{dv_{z}} / #left(#frac{1}{N_{MB}}#frac{dN_{MB}}{dv_{z}}#right)");
+      histo->SetTitle(Form("%g < p_{T} (GeV/c) < %g",ptAxis->GetBinLowEdge(firstPtBin), ptAxis->GetBinUpEdge(lastPtBin)));
       slope = ( histo->GetBinContent(histo->FindBin(meanZ+sigmaZ)) -
                histo->GetBinContent(histo->FindBin(meanZ-sigmaZ)) ) / ( 2. * sigmaZ );