Merge branch 'feature-movesplit'

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

diff --git a/PWG/muon/AliAnalysisTaskSingleMu.cxx b/PWG/muon/AliAnalysisTaskSingleMu.cxx
index 23d15c103c1697c6f181ba6d40b4845b8010b324..cd02c8f2098bd1b106b20d3bd044056e92446dae 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);
@@ -249,15 +248,27 @@ void AliAnalysisTaskSingleMu::ProcessEvent(TString physSel, const TObjArray& sel
       // (before ISR, FSR and kt kick) and the final state one.
       // The first muon is of course there only for information and should be rejected.
       // The Pythia code for initial state particles is 21
-      Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 && AliAnalysisMuonUtility::GetStatusCode(track) != 21 );
+      // When running with POWHEG, Pythia puts the hard process input of POWHEG in the stack
+      // with state 21, and then re-add it to stack before applying ISR, FSR and kt kick.
+      // This muon produces the final state muon, and its status code is 11
+      // To avoid all problems, keep only final state muon (status code <10)
+      // FIXME: is the convention valid for other generators as well?
+      Bool_t isSelected = ( istep == kStepReconstructed ) ? fMuonTrackCuts->IsSelected(track) : ( TMath::Abs(track->PdgCode()) == 13 && AliAnalysisMuonUtility::GetStatusCode(track) < 10 );
       if ( ! isSelected ) continue;
       
       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()));
@@ -450,10 +461,12 @@ void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
         }
         TH1* histo = gridSparseArray[igrid]->Project(kHvarPt, kHvarEta);
         histo->SetName(Form("hPtEta_%s_%s_%s", gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
+        histo->SetDirectory(0);
         if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
         for ( Int_t iproj=0; iproj<4; ++iproj ) {
           histo = gridSparseArray[igrid]->Project(iproj);
           histo->SetName(Form("proj%i_%s_%s_%s", iproj, gridSparseName[igrid].Data(), srcAxis->GetBinLabel(isrc), chargeNames[icharge].Data()));
+          histo->SetDirectory(0);
           if ( histo->Integral() > 0 ) histoList[icharge].Add(histo);
         } // loop on projections
       } // loop on grid sparse
@@ -487,6 +500,7 @@ void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
     axisTitle.ReplaceAll("MuMinus","#mu^{-}");
     histo->GetYaxis()->SetTitle(axisTitle.Data());
     histo->SetStats(kFALSE);
+    histo->SetDirectory(0);
     histoListRatio.Add(histo);
   }
   
@@ -568,16 +582,14 @@ void AliAnalysisTaskSingleMu::Terminate(Option_t *) {
   evtSel = Form("trigger:%s/%s", trigClassName.Data(), countPhysSel.Data());
   fEventCounters->Print("centrality",evtSel.Data(),kTRUE);
   
-  
   TString outFilename = Form("/tmp/out%s.root", GetName());
+  printf("\nCreating file %s\n", outFilename.Data());
   TFile* outFile = new TFile(outFilename.Data(),"RECREATE");
   for ( Int_t icharge=0; icharge<3; icharge++ ) {
     histoList[icharge].Write();
   }
   histoListRatio.Write();
   outFile->Close();
-  printf("\nCreating file %s\n", outFilename.Data(
-         ));
   
   ///////////////////
   // Vertex method //