Use fixed Ecut instead of pT-dependent

[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetResponseV2.cxx

diff --git a/PWGJE/AliAnalysisTaskJetResponseV2.cxx b/PWGJE/AliAnalysisTaskJetResponseV2.cxx
index 7d704aa1387db4749a914fe0a9f3402e86779643..b7a016eef42a14db42c7b3c4e53f69e315dabdf9 100644 (file)
--- a/PWGJE/AliAnalysisTaskJetResponseV2.cxx
+++ b/PWGJE/AliAnalysisTaskJetResponseV2.cxx
@@ -45,6 +45,7 @@
 
 #include "AliAODEvent.h"
 #include "AliAODJet.h"
+#include "AliAODHandler.h"
 
 #include "AliAnalysisTaskJetResponseV2.h"
 
@@ -57,6 +58,9 @@ AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2() :
 AliAnalysisTaskSE(),
   fESD(0x0),
   fAOD(0x0),
+  fAODOut(0x0),
+  fAODExtension(0x0),
+  fNonStdFile(""),
   fBackgroundBranch(""),
   fIsPbPb(kTRUE),
   fOfflineTrgMask(AliVEvent::kAny),
@@ -121,6 +125,9 @@ AliAnalysisTaskJetResponseV2::AliAnalysisTaskJetResponseV2(const char *name) :
   AliAnalysisTaskSE(name),
   fESD(0x0),
   fAOD(0x0),
+  fAODOut(0x0),
+  fAODExtension(0x0),
+  fNonStdFile(""),
   fBackgroundBranch(""),
   fIsPbPb(kTRUE),
   fOfflineTrgMask(AliVEvent::kAny),
@@ -196,8 +203,10 @@ void AliAnalysisTaskJetResponseV2::SetBranchNames(const TString &branch1, const
   fJetBranchName[1] = branch2;
   fJetBranchName[2] = branch3;
 
-  if(strlen(fJetBranchName[2].Data()) ) 
+  if(strlen(fJetBranchName[2].Data()) ) {
     fkNbranches = 3;
+    fbJets3Branches = kTRUE;
+  }
 }
 
 void AliAnalysisTaskJetResponseV2::Init()
@@ -315,10 +324,10 @@ void AliAnalysisTaskJetResponseV2::UserCreateOutputObjects()
     fhnJetsArea = NewTHnSparseF("fhnJetsArea", entries, opt);
   }
 
-  // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : pT hard bin
+  // cent : nInpTrks : jetPt(3x) : deltaPt : delta : jetArea(3x) : fraction(2x) : deltaR(1x) : pT hard bin
   if(fbJets3Branches){
-    entries = 1<<0 | 1<<1 | 1<<6 | 1<<7 | 1<<7 | 1<<14 | 1<<14 | 1<<12 | 1<<13 | 1<<13 | 1<<19 | 1<<19 | 1<<26;
-    opt = 1<<6 | 1<<7 | 1<<14;
+    entries = 1<<0 | 1<<1 | 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28 | 1<<12 | 1<<13 | 1<<29 | 1<<19 | 1<<30 | 1<<17 | 1<<26;
+    opt = 1<<6 | 1<<7 | 1<<27 | 1<<14 | 1<<28;
     fhnJets3Branches = NewTHnSparseF("fhnJets3Branches", entries, opt);
   }
 
@@ -385,7 +394,10 @@ void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
   fESD=dynamic_cast<AliESDEvent*>(InputEvent());
   if (!fESD) {
     AliError("ESD not available");
+
     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
+    //  assume that the AOD is in the general output...
+    fAODOut  = AODEvent();
   } else {
     fAOD = dynamic_cast<AliAODEvent*>(AODEvent());
   }
@@ -394,6 +406,15 @@ void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
     return;
   }
 
+  if(fNonStdFile.Length()!=0){
+    // case that we have an AOD extension we need can fetch the jets from the extended output
+    AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+    fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);    
+    if(!fAODExtension){
+      if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
+    }
+  }
+
   // -- event selection --
   fHistEvtSelection->Fill(1); // number of events before event selection
 
@@ -470,23 +491,37 @@ void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
   AliAODJetEventBackground* externalBackground = 0;
   if(!externalBackground&&fBackgroundBranch.Length()){
     externalBackground =  (AliAODJetEventBackground*)(fAOD->FindListObject(fBackgroundBranch.Data()));
+    if(!externalBackground && fAODOut) externalBackground =  (AliAODJetEventBackground*)(fAODOut->FindListObject(fBackgroundBranch.Data()));
+    if(!externalBackground && fAODExtension)  externalBackground =  (AliAODJetEventBackground*)(fAODExtension->GetAOD()->FindListObject(fBackgroundBranch.Data()));
     //if(!externalBackground)Printf("%s:%d Background branch not found %s",(char*)__FILE__,__LINE__,fBackgroundBranch.Data());;
   }
   Float_t rho = 0;
-  if(externalBackground)rho = externalBackground->GetBackground(0);
+  if(externalBackground) rho = externalBackground->GetBackground(0);
 
 
   // fetch jets
   TClonesArray *aodJets[3];
   aodJets[0] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
+  if(!aodJets[0] && fAODOut) aodJets[0] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[0].Data())); // in general: embedded jet
+  if(!aodJets[0] && fAODExtension) aodJets[0] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[0].Data()));
   aodJets[1] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet + UE version1
-  if( strlen(fJetBranchName[2].Data()) )
+  if(!aodJets[1] && fAODOut) aodJets[1] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[1].Data())); // in general: embedded jet
+  if(!aodJets[1] && fAODExtension) aodJets[1] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[1].Data()));
+  if( strlen(fJetBranchName[2].Data()) ) {
     aodJets[2] = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet + UE version2
+    if(!aodJets[2] && fAODOut) aodJets[2] = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName[2].Data())); // in general: embedded jet
+    if(!aodJets[2] && fAODExtension) aodJets[2] = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName[2].Data()));
+    if(aodJets[2]) fkNbranches=3;
+    if(fDebug>10) printf("3rd branch: %s\n",fJetBranchName[2].Data());
+  }
 
+  if(fDebug>10)  printf("fkNbranches %d\n",fkNbranches);
   for (Int_t iJetType = 0; iJetType < fkNbranches; iJetType++) {
     fListJets[iJetType]->Clear();
-    if (!aodJets[iJetType]) continue;
-      
+    if (!aodJets[iJetType]) {
+        if(fDebug) Printf("%s: no jet branch\n",fJetBranchName[iJetType].Data());
+      continue;
+    }
     if(fDebug) Printf("%s: %d jets",fJetBranchName[iJetType].Data(),aodJets[iJetType]->GetEntriesFast());
       
     for (Int_t iJet = 0; iJet < aodJets[iJetType]->GetEntriesFast(); iJet++) {
@@ -509,10 +544,12 @@ void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
                                            aMatchIndex, aPtFraction, fDebug, fMatchMaxDist, fIsPbPb?1:2);
   static TArrayI aMatchIndexv2(fListJets[0]->GetEntries());
   static TArrayF aPtFractionv2(fListJets[0]->GetEntries());
+  if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
+  if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
   if( strlen(fJetBranchName[2].Data()) ) {
-    if(aMatchIndexv2.GetSize()<fListJets[0]->GetEntries()) aMatchIndexv2.Set(fListJets[0]->GetEntries());
-    if(aPtFractionv2.GetSize()<fListJets[0]->GetEntries()) aPtFractionv2.Set(fListJets[0]->GetEntries());
-    AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()), fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()), aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2);
+    AliAnalysisHelperJetTasks::GetJetMatching(fListJets[0], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[0]->GetEntries()), 
+                                             fListJets[2], TMath::Min((Int_t)fNMatchJets,(Int_t)fListJets[2]->GetEntries()), 
+                                             aMatchIndexv2, aPtFractionv2, fDebug, fMatchMaxDist, fIsPbPb?1:2);
   }
    
   // loop over matched jets
@@ -530,7 +567,8 @@ void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
    
   for(Int_t ig=0; ig<fListJets[0]->GetEntries(); ++ig){
     ir = aMatchIndex[ig];
-    ir2 = aMatchIndexv2[ig];
+    if(aMatchIndexv2.GetSize()>ig) ir2 = aMatchIndexv2[ig];
+    else ir2 =0;
 
     //fetch jets
     jet[0] = (AliAODJet*)(fListJets[0]->At(ig));
@@ -558,7 +596,8 @@ void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
       }
     }
     fraction  = aPtFraction[ig];
-    fraction2 = aPtFractionv2[ig];
+    if(aPtFractionv2.GetSize()>ig) fraction2 = aPtFractionv2[ig];
+    else fraction2 = 0.;
 
     // jet statistics
     fHistJetSelection->Fill(1); // all probe jets
@@ -765,11 +804,10 @@ void AliAnalysisTaskJetResponseV2::UserExec(Option_t *)
     if(fbJets3Branches){
       Double_t jetEntries3Branches[14] = {
        (Double_t)centValue, (Double_t)nInputTracks, 
-       (Double_t)jetPt[0], (Double_t)jetPt[1], (Double_t)jetPt[2], 
-       (Double_t)deltaPt, (Double_t)delta,
-       (Double_t)jetArea[0], (Double_t)jetArea[1],(Double_t)jetArea[3], 
-       (Double_t)fraction, (Double_t)fraction2,
-       (Double_t)pthardbin
+       (Double_t)jetPt[0], (Double_t)jetPt[1],
+       (Double_t)jetArea[0], (Double_t)jetArea[1],
+       (Double_t)deltaPt, (Double_t)fraction, (Double_t)pthardbin,
+       (Double_t)jetPt[2],(Double_t)delta,(Double_t)jetArea[2], (Double_t)fraction2, (Double_t)deltaR
       };                                
       fhnJets3Branches->Fill(jetEntries3Branches);
     }
@@ -801,10 +839,14 @@ Int_t AliAnalysisTaskJetResponseV2::GetNInputTracks()
   }
   // use only HI event
   if(jbname.Contains("AODextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
+  else if(jbname.Contains("AODMCextraonly")) jbname.ReplaceAll("AODextraonly","AOD");
   if(jbname.Contains("AODextra")) jbname.ReplaceAll("AODextra","AOD");
+  else if(jbname.Contains("AODMCextra")) jbname.ReplaceAll("AODextra","AOD");
 
   if(fDebug) Printf("Multiplicity from jet branch %s", jbname.Data());
   TClonesArray *tmpAODjets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(jbname.Data()));
+  if(!tmpAODjets && fAODOut) tmpAODjets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(jbname.Data()));
+  if(!tmpAODjets && fAODExtension) tmpAODjets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(jbname.Data()));
   if(!tmpAODjets){
     Printf("Jet branch %s not found", jbname.Data());
     Printf("AliAnalysisTaskJetResponseV2::GetNInputTracks FAILED");
@@ -935,8 +977,9 @@ void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString
       
   case 6:
   case 7:
+  case 27:
     if(iEntry==6)label = "probe p_{T} (GeV/c)";
-    if(iEntry==7)label = "rec p_{T} (GeV/c)";
+    if(iEntry==7 || iEntry==27)label = "rec p_{T} (GeV/c)";
     if(hr){
       nbins = 300;
       xmin = -50.;
@@ -983,8 +1026,9 @@ void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString
       
   case 12:
   case 13:
+  case 29:
     if(iEntry==12)label = "probe area";
-    if(iEntry==13)label = "rec area";
+    if(iEntry==13 || iEntry==29)label = "rec area";
     if(hr){
       nbins = 100;
       xmin = 0.;
@@ -997,7 +1041,9 @@ void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString
     break;
       
   case 14:
-    label = "#Delta p_{T}";
+  case 28:
+    if(iEntry==14) label = "#delta p_{T}";
+    if(iEntry==28) label = "#delta";
     if(hr){
       nbins = 241;
       xmin = -120.5;
@@ -1066,6 +1112,7 @@ void AliAnalysisTaskJetResponseV2::GetDimParams(Int_t iEntry, Bool_t hr, TString
       
       
   case 19:
+  case 30:
     label = "fraction";
     if(hr){
       nbins = 52;