Fixes for not filled histograms and calculation of Dijet binning
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskFragmentationFunction.cxx
index 708c35648fa7c64df0094e3fac8e22711152efa9..ce780956c8062941eeeeada39b97bba5d0d6b280 100644 (file)
@@ -27,6 +27,9 @@
 #include "TH2F.h"
 #include "TString.h"
 #include "THnSparse.h"
+#include "TProfile.h"
+#include "TFile.h"
+#include "TKey.h"
 
 #include "AliAODInputHandler.h" 
 #include "AliAODHandler.h" 
@@ -57,6 +60,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction()
    ,fJetTypeGen(0)
    ,fJetTypeRecEff(0)
    ,fFilterMask(0)
+   ,fUsePhysicsSelection(kTRUE)
    ,fTrackPtCut(0)
    ,fTrackEtaMin(0)
    ,fTrackEtaMax(0)
@@ -195,6 +199,10 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction()
    ,fh1VertexNContributors(0)
    ,fh1VertexZ(0)
    ,fh1EvtMult(0)
+   ,fh1Xsec(0)
+   ,fh1Trials(0)
+   ,fh1PtHard(0)
+   ,fh1PtHardTrials(0)
    ,fh1nRecJetsCuts(0)
    ,fh1nGenJets(0)
    ,fh1nRecEffJets(0)
@@ -216,6 +224,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fJetTypeGen(0)
   ,fJetTypeRecEff(0)
   ,fFilterMask(0)
+  ,fUsePhysicsSelection(kTRUE)
   ,fTrackPtCut(0)
   ,fTrackEtaMin(0)
   ,fTrackEtaMax(0)
@@ -354,6 +363,10 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fh1VertexNContributors(0)
   ,fh1VertexZ(0)
   ,fh1EvtMult(0)
+  ,fh1Xsec(0)
+  ,fh1Trials(0)
+  ,fh1PtHard(0)
+  ,fh1PtHardTrials(0)
   ,fh1nRecJetsCuts(0)
   ,fh1nGenJets(0)
   ,fh1nRecEffJets(0)
@@ -379,6 +392,7 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fJetTypeGen(copy.fJetTypeGen)
   ,fJetTypeRecEff(copy.fJetTypeRecEff)
   ,fFilterMask(copy.fFilterMask)
+  ,fUsePhysicsSelection(copy.fUsePhysicsSelection)
   ,fTrackPtCut(copy.fTrackPtCut)
   ,fTrackEtaMin(copy.fTrackEtaMin)
   ,fTrackEtaMax(copy.fTrackEtaMax)
@@ -517,6 +531,10 @@ AliAnalysisTaskFragmentationFunction::AliAnalysisTaskFragmentationFunction(const
   ,fh1VertexNContributors(copy.fh1VertexNContributors)
   ,fh1VertexZ(copy.fh1VertexZ)
   ,fh1EvtMult(copy.fh1EvtMult)
+  ,fh1Xsec(copy.fh1Xsec)
+  ,fh1Trials(copy.fh1Trials)
+  ,fh1PtHard(copy.fh1PtHard)  
+  ,fh1PtHardTrials(copy.fh1PtHardTrials)  
   ,fh1nRecJetsCuts(copy.fh1nRecJetsCuts)
   ,fh1nGenJets(copy.fh1nGenJets)
   ,fh1nRecEffJets(copy.fh1nRecEffJets)
@@ -544,6 +562,7 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper
     fJetTypeGen                   = o.fJetTypeGen;
     fJetTypeRecEff                = o.fJetTypeRecEff;
     fFilterMask                   = o.fFilterMask;
+    fUsePhysicsSelection          = o.fUsePhysicsSelection;
     fTrackPtCut                   = o.fTrackPtCut;
     fTrackEtaMin                  = o.fTrackEtaMin;
     fTrackEtaMax                  = o.fTrackEtaMax;
@@ -682,6 +701,10 @@ AliAnalysisTaskFragmentationFunction& AliAnalysisTaskFragmentationFunction::oper
     fh1VertexNContributors        = o.fh1VertexNContributors;
     fh1VertexZ                    = o.fh1VertexZ;
     fh1EvtMult                    = o.fh1EvtMult;
+    fh1Xsec                       = o.fh1Xsec;
+    fh1Trials                     = o.fh1Trials;
+    fh1PtHard                     = o.fh1PtHard;
+    fh1PtHardTrials               = o.fh1PtHardTrials;
     fh1nRecJetsCuts               = o.fh1nRecJetsCuts;
     fh1nGenJets                   = o.fh1nGenJets; 
     fh1nRecEffJets                = o.fh1nRecEffJets;
@@ -1734,6 +1757,88 @@ void AliAnalysisTaskFragmentationFunction::AliFragFuncIntraJetHistos::AddToOutpu
 
 }
 
+
+Bool_t AliAnalysisTaskFragmentationFunction::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // (taken from AliAnalysisTaskJetSpectrum)
+  // 
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  Double_t xsection = 0;
+  UInt_t   ntrials  = 0;
+  Float_t   ftrials  = 0;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+    if(!fh1Xsec||!fh1Trials){
+      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+      return kFALSE;
+    }
+
+    TString fileName(curfile->GetName());
+    if(fileName.Contains("AliESDs.root")){
+        fileName.ReplaceAll("AliESDs.root", "");
+    }
+    else if(fileName.Contains("AliAOD.root")){
+        fileName.ReplaceAll("AliAOD.root", "");
+    }
+    else if(fileName.Contains("AliAODs.root")){
+        fileName.ReplaceAll("AliAODs.root", "");
+    }
+    else if(fileName.Contains("galice.root")){
+        // for running with galice and kinematics alone...                      
+        fileName.ReplaceAll("galice.root", "");
+    }
+    TFile *fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec.root"));
+    if(!fxsec){
+      if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec.root"));
+      // next trial fetch the histgram file
+      fxsec = TFile::Open(Form("%s%s",fileName.Data(),"pyxsec_hists.root"));
+      if(!fxsec){
+       // not a severe condition
+       if(fDebug>0)Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,Form("%s%s",fileName.Data(),"pyxsec_hists.root"));        
+       return kTRUE;
+      }
+      else{
+       // find the tlist we want to be independtent of the name so use the Tkey
+       TKey* key = (TKey*)fxsec->GetListOfKeys()->At(0); 
+       if(!key){
+         if(fDebug>0)Printf("%s:%d key not found in the file",(char*)__FILE__,__LINE__);       
+         return kTRUE;
+       }
+       TList *list = dynamic_cast<TList*>(key->ReadObj());
+       if(!list){
+         if(fDebug>0)Printf("%s:%d key is not a tlist",(char*)__FILE__,__LINE__);      
+         return kTRUE;
+       }
+       xsection = ((TProfile*)list->FindObject("h1Xsec"))->GetBinContent(1);
+       ftrials  = ((TH1F*)list->FindObject("h1Trials"))->GetBinContent(1);
+      }
+    }
+    else{
+      TTree *xtree = (TTree*)fxsec->Get("Xsection");
+      if(!xtree){
+       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+      }
+      xtree->SetBranchAddress("xsection",&xsection);
+      xtree->SetBranchAddress("ntrials",&ntrials);
+      ftrials = ntrials;
+      xtree->GetEntry(0);
+    }
+    fh1Xsec->Fill("<#sigma>",xsection);
+    fh1Trials->Fill("#sum{ntrials}",ftrials);
+  }
+  
+  return kTRUE;
+}
+
+
+
 //__________________________________________________________________
 void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
 {
@@ -1789,7 +1894,15 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
   fh1EvtSelection            = new TH1F("fh1EvtSelection", "Event Selection", 6, -0.5, 5.5);
   fh1VertexNContributors     = new TH1F("fh1VertexNContributors", "Vertex N contributors", 11,-.5, 10.5);
   fh1VertexZ                 = new TH1F("fh1VertexZ", "Vertex z distribution", 30, -15., 15.);
-  fh1EvtMult                = new TH1F("fh1EvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",100,0.,100.);
+  fh1EvtMult                = new TH1F("fh1EvtMult","Event multiplicity, track pT cut > 150 MeV/c, |#eta| < 0.9",120,0.,120.);
+
+  fh1Xsec                    = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+  fh1Trials                  = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+  fh1PtHard                  = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",350,-.5,349.5);
+  fh1PtHardTrials            = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",350,-.5,349.5);
+
   fh1nRecJetsCuts            = new TH1F("fh1nRecJetsCuts","reconstructed jets per event",10,-0.5,9.5);
   fh1nGenJets                = new TH1F("fh1nGenJets","generated jets per event",10,-0.5,9.5);
   fh1nRecEffJets             = new TH1F("fh1nRecEffJets","reconstruction effiency: jets per event",10,-0.5,9.5);
@@ -2001,31 +2114,47 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
   fQADiJetHistosRecCuts->DefineQADiJetHistos();
   fQADiJetHistosGen->DefineQADiJetHistos();
 
+  Bool_t genJets    = (fJetTypeGen != kJetsUndef) ? kTRUE : kFALSE;
+  Bool_t genTracks  = (fTrackTypeGen != kTrackUndef) ? kTRUE : kFALSE;
+  Bool_t recJetsEff = (fJetTypeRecEff != kJetsUndef) ? kTRUE : kFALSE;
+
+
   const Int_t saveLevel = 5;
   if(saveLevel>0){
     fCommonHistList->Add(fh1EvtSelection);
     fFFHistosRecCuts->AddToOutput(fCommonHistList);
     fFFHistosRecLeading->AddToOutput(fCommonHistList);
     fFFHistosRecLeadingTrack->AddToOutput(fCommonHistList);
-    fFFHistosGen->AddToOutput(fCommonHistList);
-    fFFHistosGenLeading->AddToOutput(fCommonHistList);
-    fFFHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+
+    if(genJets && genTracks){
+       fFFHistosGen->AddToOutput(fCommonHistList);
+       fFFHistosGenLeading->AddToOutput(fCommonHistList);
+       fFFHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+
+       fCommonHistList->Add(fh1Xsec);
+       fCommonHistList->Add(fh1Trials);
+       fCommonHistList->Add(fh1PtHard);
+       fCommonHistList->Add(fh1PtHardTrials);
+    }
   }
   if(saveLevel>1){
     fQATrackHistosRec->AddToOutput(fCommonHistList);
     fQATrackHistosRecCuts->AddToOutput(fCommonHistList);
-    fQATrackHistosGen->AddToOutput(fCommonHistList);
+    if(genTracks) fQATrackHistosGen->AddToOutput(fCommonHistList);
     
     fQAJetHistosRec->AddToOutput(fCommonHistList);
     fQAJetHistosRecCuts->AddToOutput(fCommonHistList);
     fQAJetHistosRecCutsLeading->AddToOutput(fCommonHistList);
-    fQAJetHistosGen->AddToOutput(fCommonHistList);
-    fQAJetHistosGenLeading->AddToOutput(fCommonHistList);
-    fQAJetHistosRecEffLeading->AddToOutput(fCommonHistList); 
+    if(recJetsEff) fQAJetHistosRecEffLeading->AddToOutput(fCommonHistList); 
+
+    if(genJets){
+       fQAJetHistosGen->AddToOutput(fCommonHistList);
+       fQAJetHistosGenLeading->AddToOutput(fCommonHistList);
+    }
 
     fCommonHistList->Add(fh1EvtMult);
     fCommonHistList->Add(fh1nRecJetsCuts);
-    fCommonHistList->Add(fh1nGenJets);
+    if(genJets) fCommonHistList->Add(fh1nGenJets);
   }
   if(saveLevel>2){
     fCommonHistList->Add(fh1VertexNContributors);
@@ -2035,22 +2164,31 @@ void AliAnalysisTaskFragmentationFunction::UserCreateOutputObjects()
     fIJHistosRecCuts->AddToOutput(fCommonHistList);
     fIJHistosRecLeading->AddToOutput(fCommonHistList);
     fIJHistosRecLeadingTrack->AddToOutput(fCommonHistList);
-    fIJHistosGen->AddToOutput(fCommonHistList);
-    fIJHistosGenLeading->AddToOutput(fCommonHistList);
-    fIJHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+
+    if(genJets && genTracks){
+       fIJHistosGen->AddToOutput(fCommonHistList);
+       fIJHistosGenLeading->AddToOutput(fCommonHistList);
+       fIJHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+    }
   }
   if(saveLevel>4){
     fFFDiJetHistosRecCuts->AddToOutput(fCommonHistList);
     fFFDiJetHistosRecLeading->AddToOutput(fCommonHistList);
     fFFDiJetHistosRecLeadingTrack->AddToOutput(fCommonHistList);
-    fFFDiJetHistosGen->AddToOutput(fCommonHistList);
-    fFFDiJetHistosGenLeading->AddToOutput(fCommonHistList);
-    fFFDiJetHistosGenLeadingTrack->AddToOutput(fCommonHistList);
     fQADiJetHistosRecCuts->AddToOutput(fCommonHistList);
-    fQADiJetHistosGen->AddToOutput(fCommonHistList); 
-    fCommonHistList->Add(fhnSingleTrackRecEffHisto);
-    fCommonHistList->Add(fhnJetTrackRecEffHisto);
-    fCommonHistList->Add(fh1nRecEffJets);
+
+    if(genJets && genTracks){
+        fFFDiJetHistosGen->AddToOutput(fCommonHistList);
+        fFFDiJetHistosGenLeading->AddToOutput(fCommonHistList);
+        fFFDiJetHistosGenLeadingTrack->AddToOutput(fCommonHistList);
+        fQADiJetHistosGen->AddToOutput(fCommonHistList);
+    }
+
+    if(recJetsEff && genTracks){
+       fCommonHistList->Add(fhnSingleTrackRecEffHisto);
+       fCommonHistList->Add(fhnJetTrackRecEffHisto);
+       fCommonHistList->Add(fh1nRecEffJets);
+    }
   }
 
   // =========== Switch on Sumw2 for all histos ===========
@@ -2092,7 +2230,7 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
     fh1EvtSelection->Fill(1.);
   } else {
     fh1EvtSelection->Fill(0.);
-    if(inputHandler->InheritsFrom("AliESDInputHandler")){ // PhysicsSelection only with ESD input
+    if(inputHandler->InheritsFrom("AliESDInputHandler") && fUsePhysicsSelection){ // PhysicsSelection only with ESD input
       if (fDebug > 1 ) Printf(" Trigger Selection: event REJECTED ... ");
       PostData(1, fCommonHistList);
       return;
@@ -2162,6 +2300,26 @@ void AliAnalysisTaskFragmentationFunction::UserExec(Option_t *)
   }
   if (fDebug > 1) Printf("%s:%d primary vertex selection: event ACCEPTED ...",(char*)__FILE__,__LINE__); 
   fh1EvtSelection->Fill(5.);
+
+
+  //___ get MC information __________________________________________________________________
+
+  Double_t ptHard = 0.;
+  Double_t nTrials = 1; // trials for MC trigger weight for real data
+  AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(fMCEvent);
+  if(!pythiaGenHeader){
+     if(fJetTypeGen != kJetsUndef && fTrackTypeGen != kTrackUndef){
+        Printf("%s:%d no pythiaGenHeader found", (char*)__FILE__,__LINE__);
+        return;
+     }
+  } else {
+     nTrials = pythiaGenHeader->Trials();
+     ptHard  = pythiaGenHeader->GetPtHard();
+
+     fh1PtHard->Fill(ptHard);
+     fh1PtHardTrials->Fill(ptHard,nTrials);
+  }
   
   
   //___ fetch jets __________________________________________________________________________
@@ -2729,9 +2887,9 @@ Double_t AliAnalysisTaskFragmentationFunction::GetDiJetBin(Double_t invMass, Dou
        {
          jetBin = fDiJetJetInvMassMin + i*stepInvMass/2.;
          if(((fDiJetJetInvMassMin+i*stepInvMass) <= invMass) &&
-            (fDiJetJetInvMassMin + (i+1)*stepInvMass) > invMass) jetBinOk = jetBin;
+            (fDiJetJetInvMassMin + (i+1)*stepInvMass) > invMass) {jetBinOk = jetBin; break;}
+          else jetBinOk = -1.;
        }
-      jetBinOk = -1.;
     }
   else if (kindBins == 3)
     {
@@ -2739,9 +2897,9 @@ Double_t AliAnalysisTaskFragmentationFunction::GetDiJetBin(Double_t invMass, Dou
        {
          jetBin = fDiJetJetPtMin + i*stepPt/2.;
          if(((fDiJetJetPtMin+i*stepPt) <= EtMean) &&
-            (fDiJetJetPtMin + (i+1)*stepPt) > EtMean) jetBinOk = jetBin;
+            (fDiJetJetPtMin + (i+1)*stepPt) > EtMean) {jetBinOk = jetBin; break;}
+          else jetBinOk = -1.;
        }
-      jetBinOk = -1.;
     }
   else if (kindBins == 2)
     {
@@ -2749,9 +2907,9 @@ Double_t AliAnalysisTaskFragmentationFunction::GetDiJetBin(Double_t invMass, Dou
        {
          jetBin = fDiJetJetPtMin + i*stepPt/2.;
          if(((fDiJetJetPtMin+i*stepPt) <= leadingJetPt) &&
-            (fDiJetJetPtMin + (i+1)*stepPt) > leadingJetPt) jetBinOk = jetBin;
+            (fDiJetJetPtMin + (i+1)*stepPt) > leadingJetPt) {jetBinOk = jetBin; break;}
+          else jetBinOk = -1.;
        }
-      jetBinOk = -1.;
     }
   else {Printf("WARNING: kindBins wrongly set ! Please make sure to call SetKindSlices() and set the kind parameter to 1, 2 or 3.\n");}
 
@@ -2783,7 +2941,7 @@ Int_t AliAnalysisTaskFragmentationFunction::GetListOfTracks(TList *list, Int_t t
   if(type==kTrackUndef) return 0;
   
   Int_t iCount = 0;
-  if(type==kTrackAODCuts || type==kTrackAOD){
+  if(type==kTrackAODCuts || type==kTrackAODQualityCuts || type==kTrackAOD){
 
     // all rec. tracks, esd filter mask, eta range
     if(!fAOD) return -1;
@@ -2791,11 +2949,13 @@ Int_t AliAnalysisTaskFragmentationFunction::GetListOfTracks(TList *list, Int_t t
     for(Int_t it=0; it<fAOD->GetNumberOfTracks(); ++it){
       AliAODTrack *tr = fAOD->GetTrack(it);
       
-      if(type == kTrackAODCuts){
+      if(type == kTrackAODCuts || type==kTrackAODQualityCuts ){
        if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))   continue;
-       if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
-       if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
-       if(tr->Pt()  < fTrackPtCut) continue;
+       if(type == kTrackAODCuts){
+         if(tr->Eta() < fTrackEtaMin || tr->Eta() > fTrackEtaMax) continue;
+         if(tr->Phi() < fTrackPhiMin || tr->Phi() > fTrackPhiMax) continue;
+         if(tr->Pt()  < fTrackPtCut) continue;
+       }
       }
       list->Add(tr);
       iCount++;