Enable writing of a AOD cluster branch (of tpe AliAODJet) with track references,...

diff --git a/JETAN/AliAnalysisTaskJetCluster.cxx b/JETAN/AliAnalysisTaskJetCluster.cxx
index 3b693ad3dc72255fcf6dd3c67204d0420f95abfa..d737ca52866bf06897c02f2fa72cdaa4f47985b9 100644 (file)
--- a/JETAN/AliAnalysisTaskJetCluster.cxx
+++ b/JETAN/AliAnalysisTaskJetCluster.cxx
@@ -25,7 +25,7 @@
 #include <TSystem.h>
 #include <TInterpreter.h>
 #include <TChain.h>
-#include <TRandom.h>
+#include <TRefArray.h>
 #include <TFile.h>
 #include <TKey.h>
 #include <TH1F.h>
@@ -42,8 +42,6 @@
 #include "AliJetFinder.h"
 #include "AliJetHeader.h"
 #include "AliJetReader.h"
-#include "AliJetReaderHeader.h"
-#include "AliUA1JetHeaderV1.h"
 #include "AliESDEvent.h"
 #include "AliAODEvent.h"
 #include "AliAODHandler.h"
@@ -69,17 +67,26 @@
 
 ClassImp(AliAnalysisTaskJetCluster)
 
+AliAnalysisTaskJetCluster::~AliAnalysisTaskJetCluster(){
+  delete fRef;
+}
+
 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
   fAOD(0x0),
+  fAODExtension(0x0),
+  fRef(new TRefArray),
   fUseAODTrackInput(kFALSE),
   fUseAODMCInput(kFALSE),
   fUseGlobalSelection(kFALSE),
   fFilterMask(0),
   fTrackTypeRec(kTrackUndef),
-  fTrackTypeGen(kTrackUndef),
-  fAvgTrials(1),
+  fTrackTypeGen(kTrackUndef),  fAvgTrials(1),
   fExternalWeight(1),    
   fRecEtaWindow(0.5),
+  fTrackPtCut(0.),                                                     
+  fJetOutputMinPt(1),
+  fNonStdBranch(""),
+  fNonStdFile(""),
   fRparam(1.0), 
   fAlgorithm(fastjet::kt_algorithm),
   fStrategy(fastjet::Best),
@@ -141,6 +148,8 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(): AliAnalysisTaskSE(),
 AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
   AliAnalysisTaskSE(name),
   fAOD(0x0),
+  fAODExtension(0x0),
+  fRef(new TRefArray),
   fUseAODTrackInput(kFALSE),
   fUseAODMCInput(kFALSE),
   fUseGlobalSelection(kFALSE),
@@ -150,6 +159,10 @@ AliAnalysisTaskJetCluster::AliAnalysisTaskJetCluster(const char* name):
   fAvgTrials(1),
   fExternalWeight(1),    
   fRecEtaWindow(0.5),
+  fTrackPtCut(0.),                                                     
+  fJetOutputMinPt(1),
+  fNonStdBranch(""),
+  fNonStdFile(""),
   fRparam(1.0), 
   fAlgorithm(fastjet::kt_algorithm),
   fStrategy(fastjet::Best),
@@ -227,11 +240,50 @@ void AliAnalysisTaskJetCluster::UserCreateOutputObjects()
   //
 
 
+
+
   // Connect the AOD
 
 
   if (fDebug > 1) printf("AnalysisTaskJetCluster::UserCreateOutputObjects() \n");
 
+
+
+  if(fNonStdBranch.Length()!=0)
+    {
+      // only create the output branch if we have a name
+      // Create a new branch for jets...
+      //  -> cleared in the UserExec....
+      // here we can also have the case that the brnaches are written to a separate file
+
+      TClonesArray *tca = new TClonesArray("AliAODJet", 0);
+      tca->SetName(fNonStdBranch.Data());
+      AddAODBranch("TClonesArray",&tca,fNonStdFile.Data());
+      if(fNonStdFile.Length()!=0){
+       // 
+       // case that we have an AOD extension we need to fetch the jets from the extended output
+       // we identifay the extension aod event by looking for the branchname
+       AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
+       TObjArray* extArray = aodH->GetExtensions();
+       if (extArray) {
+         TIter next(extArray);
+         while ((fAODExtension=(AliAODExtension*)next())){
+           TObject *obj = fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data());
+           if(fDebug>10){
+             Printf("%s:%d Dumping..",(char*)__FILE__,__LINE__);
+             fAODExtension->GetAOD()->Dump();
+           }
+           if(obj){
+             if(fDebug>1)Printf("AODExtension found for %s",fNonStdBranch.Data());
+             break;
+           }
+           fAODExtension = 0;
+         }
+       }
+      }
+    }
+
+
   OpenFile(1);
   if(!fHistList)fHistList = new TList();
 
@@ -451,6 +503,15 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
     return;
   }
 
+  // handle and reset the output jet branch 
+  // only need this once
+  TClonesArray* jarray = 0;  
+  if(fNonStdBranch.Length()!=0) {
+    if(AODEvent())jarray = (TClonesArray*)(AODEvent()->FindListObject(fNonStdBranch.Data()));
+    if(!jarray)jarray = (TClonesArray*)(fAODExtension->GetAOD()->FindListObject(fNonStdBranch.Data()));
+    if(jarray)jarray->Delete();    // this is our responsibility, clear before filling again
+  }
+
   //
   // Execute analysis for current event
   //
@@ -476,11 +537,13 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
   }
   
   Bool_t selectEvent =  false;
-  Bool_t physicsSelection = (fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
+  Bool_t physicsSelection = true;// handled by the framework(fInputHandler->IsEventSelected()&AliVEvent::kMB)==AliVEvent::kMB;
 
   if(fAOD){
     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
-    if(vtxAOD->GetNContributors()>0){
+    TString vtxTitle(vtxAOD->GetTitle());
+
+    if(vtxAOD->GetNContributors()>2&&!vtxTitle.Contains("TPCVertex")){
        Float_t zvtx = vtxAOD->GetZ();
        Float_t yvtx = vtxAOD->GetY();
        Float_t xvtx = vtxAOD->GetX();
@@ -545,6 +608,15 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
     jInpRan.set_user_index(i);
     inputParticlesRecRan.push_back(jInpRan);
 
+
+    // fill the tref array, only needed when we write out jets
+    if(jarray){
+      if(i == 0){
+       fRef->Delete(); // make sure to delete before placement new...
+       new(fRef) TRefArray(TProcessID::GetProcessWithUID(vp));
+      }
+      fRef->Add(vp);
+    }
   }
 
   if(inputParticlesRec.size()==0){
@@ -565,11 +637,14 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
 
  // loop over all jets an fill information, first on is the leading jet
 
- Int_t nRecOver = inclusiveJets.size();
- Int_t nRec     = inclusiveJets.size();
- if(inclusiveJets.size()>0){
-   AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
+  Int_t nRecOver = inclusiveJets.size();
+  Int_t nRec     = inclusiveJets.size();
+  if(inclusiveJets.size()>0){
+    AliAODJet leadingJet (sortedJets[0].px(), sortedJets[0].py(), sortedJets[0].pz(), sortedJets[0].E());
     Float_t pt = leadingJet.Pt();
+    Int_t nAodOutJets = 0;
+    Int_t nAodOutTracks = 0;
+    AliAODJet *aodOutJet = 0;
 
     Int_t iCount = 0;
     for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
@@ -588,11 +663,11 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
     if(phi<0)phi+=TMath::Pi()*2.;    
     Float_t eta = leadingJet.Eta();
     pt = leadingJet.Pt();
-
+    
     // correlation of leading jet with tracks
     TIterator *recIter = recParticles.MakeIterator();
     AliVParticle *tmpRecTrack = (AliVParticle*)(recIter->Next());  
-
+    
     recIter->Reset();
     while((tmpRecTrack = (AliVParticle*)(recIter->Next()))){
       Float_t tmpPt = tmpRecTrack->Pt();
@@ -605,54 +680,73 @@ void AliAnalysisTaskJetCluster::UserExec(Option_t */*option*/)
       //      Float_t dEta = eta - tmpRecTrack->Eta();
       fh2TracksLeadingJetPhiPt->Fill(dPhi,pt);
       fh2TracksLeadingJetPhiPtW->Fill(dPhi,pt,tmpPt);
-   }  
+    }  
+    
+   
+    
+   for(int j = 0; j < nRec;j++){
+     AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
+     aodOutJet = 0;
+     nAodOutTracks = 0;
+     Float_t tmpPt = tmpRec.Pt();  
+     fh1PtJetsRecIn->Fill(tmpPt);
+     // Fill Spectra with constituents
+     vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
+     fh1NConstRec->Fill(constituents.size());
+     fh2PtNch->Fill(nCh,tmpPt);
+     fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
+     fh2NConstPt->Fill(tmpPt,constituents.size());
+     // loop over constiutents and fill spectrum
+
+     // Add the jet information and the track references to the output AOD
+     
+     if(tmpPt>fJetOutputMinPt&&jarray){
+       aodOutJet =  new ((*jarray)[nAodOutJets++]) AliAODJet(tmpRec);
+     }
 
+     
+     for(unsigned int ic = 0; ic < constituents.size();ic++){
+       AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
+       fh1PtJetConstRec->Fill(part->Pt());
+       if(aodOutJet){
+        aodOutJet->AddTrack(fRef->At(constituents[ic].user_index()));
+       }
+       if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
+     }
 
 
-    for(int j = 0; j < nRec;j++){
-      AliAODJet tmpRec (sortedJets[j].px(), sortedJets[j].py(), sortedJets[j].pz(), sortedJets[j].E());
-      Float_t tmpPt = tmpRec.Pt();
-      fh1PtJetsRecIn->Fill(tmpPt);
-      // Fill Spectra with constituents
-      vector<fastjet::PseudoJet> constituents = clustSeq.constituents(sortedJets[j]);
-      fh1NConstRec->Fill(constituents.size());
-      fh2PtNch->Fill(nCh,tmpPt);
-      fh2PtNchN->Fill(nCh,tmpPt,constituents.size());
-      fh2NConstPt->Fill(tmpPt,constituents.size());
-      // loop over constiutents and fill spectrum
-      for(unsigned int ic = 0; ic < constituents.size();ic++){
-       AliVParticle *part = (AliVParticle*)recParticles.At(constituents[ic].user_index());
-       fh1PtJetConstRec->Fill(part->Pt());
-       if(j==0)fh1PtJetConstLeadingRec->Fill(part->Pt());
-      }
 
-      // correlation
-      Float_t tmpPhi =  tmpRec.Phi();
-      Float_t tmpEta =  tmpRec.Eta();
-      if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
+     
+     
 
-      if(j==0){
-       fh1PtJetsLeadingRecIn->Fill(tmpPt);
-       fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
-       fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
-       fh1NConstLeadingRec->Fill(constituents.size());
-       fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
-       continue;
-      }
-      fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
-      fh2JetEtaPt->Fill(tmpEta,tmpPt);
-      Float_t dPhi = phi - tmpPhi;
-      if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
-      if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
-      Float_t dEta = eta - tmpRec.Eta();
-      fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
-      fh2JetsLeadingPhiPt->Fill(dPhi,pt);
-      fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
-    }
-    delete recIter;
- }
- 
 
+     
+     // correlation
+     Float_t tmpPhi =  tmpRec.Phi();
+     Float_t tmpEta =  tmpRec.Eta();
+     if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;    
+     
+     if(j==0){
+       fh1PtJetsLeadingRecIn->Fill(tmpPt);
+       fh2LeadingJetPhiEta->Fill(tmpPhi,tmpEta);
+       fh2LeadingJetEtaPt->Fill(tmpEta,tmpPt);
+       fh1NConstLeadingRec->Fill(constituents.size());
+       fh2NConstLeadingPt->Fill(tmpPt,constituents.size());
+       continue;
+     }
+     fh2JetPhiEta->Fill(tmpRec.Phi(),tmpEta);
+     fh2JetEtaPt->Fill(tmpEta,tmpPt);
+     Float_t dPhi = phi - tmpPhi;
+     if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
+     if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();      
+     Float_t dEta = eta - tmpRec.Eta();
+     fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
+     fh2JetsLeadingPhiPt->Fill(dPhi,pt);
+     fh2JetsLeadingPhiPtW->Fill(dPhi,pt,tmpPt);
+   }
+   delete recIter;
+  }
+ 
  // fill track information
  Int_t nTrackOver = recParticles.GetSize();
   // do the same for tracks and jets

diff --git a/JETAN/AliAnalysisTaskJetCluster.h b/JETAN/AliAnalysisTaskJetCluster.h
index 01db107d1f3c3e331037b0a06a4b51fc792938ba..f3b4f7148bccf4b2df20819ded0dda12734e37e8 100644 (file)
--- a/JETAN/AliAnalysisTaskJetCluster.h
+++ b/JETAN/AliAnalysisTaskJetCluster.h
@@ -19,6 +19,7 @@
 class AliJetHeader;
 class AliESDEvent;
 class AliAODEvent;
+class AliAODExtension;
 class AliAODJet;
 class AliGenPythiaEventHeader;
 class AliCFManager;
@@ -29,7 +30,7 @@ class TH2F;
 class TH1F;
 class TH3F;
 class TProfile;
-
+class TRefArray;
 
 
 class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
@@ -37,7 +38,7 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
  public:
     AliAnalysisTaskJetCluster();
     AliAnalysisTaskJetCluster(const char* name);
-    virtual ~AliAnalysisTaskJetCluster() {;}
+    virtual ~AliAnalysisTaskJetCluster();
     // Implementation of interface methods
     virtual void UserCreateOutputObjects();
     virtual void Init();
@@ -52,8 +53,13 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
     virtual void SetRecEtaWindow(Float_t f){fRecEtaWindow = f;}
     virtual void SetTrackTypeGen(Int_t i){fTrackTypeGen = i;}
     virtual void SetTrackTypeRec(Int_t i){fTrackTypeRec = i;}
+    virtual void SetTrackPtCut(Float_t x){fTrackPtCut = x;}
     virtual void SetFilterMask(UInt_t i){fFilterMask = i;}
 
+    virtual void SetJetOutputBranch(const char *c){fNonStdBranch = c;}
+    virtual void SetJetOutputFile(const char *c){fNonStdFile = c;}
+    virtual void SetJetOutputMinPt(Float_t x){fJetOutputMinPt = x;}
+
     // for Fast Jet
     fastjet::JetAlgorithm        GetAlgorithm()         const {return fAlgorithm;}
     fastjet::Strategy            GetStrategy()          const {return fStrategy;}
@@ -86,17 +92,26 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
 
     Int_t GetListOfTracks(TList *list,Int_t type);
 
-    AliAODEvent  *fAOD; // where we take the jets from can be input or output AOD
-
+    AliAODEvent     *fAOD;                // ! where we take the jets from can be input or output AOD
+    AliAODExtension *fAODExtension;       // ! AOD extension in case we write a non-sdt branch to a separate file and the aod is standard
+    TRefArray       *fRef;               // ! trefarray for track references within the jet
     Bool_t        fUseAODTrackInput;      // take track from input AOD not from ouptu AOD
     Bool_t        fUseAODMCInput;         // take MC from input AOD not from ouptu AOD
     Bool_t        fUseGlobalSelection;    // Limit the eta of the generated jets
-    UInt_t        fFilterMask;             // filter bit for slecected tracks
+    UInt_t        fFilterMask;            // filter bit for slecected tracks
     Int_t         fTrackTypeRec;          // type of tracks used for FF 
     Int_t         fTrackTypeGen;          // type of tracks used for FF 
     Float_t       fAvgTrials;             // Average nimber of trials
     Float_t       fExternalWeight;        // external weight
     Float_t       fRecEtaWindow;          // eta window used for corraltion plots between rec and gen 
+    Float_t       fTrackPtCut;            // minimum track pt to be accepted
+    Float_t       fJetOutputMinPt;        // minimum p_t for jets to be written out
+
+    // output configurartion
+    TString       fNonStdBranch;      // the name of the non-std branch name, if empty no branch is filled
+    TString       fNonStdFile;        // The optional name of the output file the non-std brnach is written to
+    
+
     // Fast jet
     Double_t fRparam;
     fastjet::JetAlgorithm fAlgorithm; //fastjet::kt_algorithm
@@ -160,7 +175,7 @@ class AliAnalysisTaskJetCluster : public AliAnalysisTaskSE
     TList *fHistList; // Output list
    
 
-    ClassDef(AliAnalysisTaskJetCluster, 3) // Analysis task for standard jet analysis
+    ClassDef(AliAnalysisTaskJetCluster, 4) 
 };
  
 #endif

diff --git a/PWG4/macros/AddTaskJetCluster.C b/PWG4/macros/AddTaskJetCluster.C
index c29aaf4fd79ae93484dc8021e3d7ccb01a59abd1..6e924ccabeec75b79df2acf3a19580c33fa1b235 100644 (file)
--- a/PWG4/macros/AddTaskJetCluster.C
+++ b/PWG4/macros/AddTaskJetCluster.C
@@ -84,8 +84,10 @@ AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec,char* bGen ,UInt_t filte
    }\r
    else if (typeRec.Contains("AOD")) {\r
      pwg4spec->SetTrackTypeRec(AliAnalysisTaskJetCluster::kTrackAOD);\r
+     pwg4spec->SetTrackPtCut(0.15);\r
    }\r
 \r
+\r
    if(typeGen.Contains("AODMC2b")){// work down from the top AODMC2b -> AODMC2 -> AODMC -> AOD\r
      pwg4spec->SetTrackTypeGen(AliAnalysisTaskJetCluster::kTrackAODMCChargedAcceptance);\r
    }\r
@@ -119,6 +121,11 @@ AliAnalysisTaskJetCluster *AddTaskJetCluster(char* bRec,char* bGen ,UInt_t filte
    }\r
 \r
    \r
+   if(TMath::Abs(radius-0.4)<0.01){\r
+     pwg4spec->SetJetOutputBranch(Form("clusters%s_%s%s",bRec,jf,cRadius));\r
+     pwg4spec->SetJetOutputMinPt(1); // store only jets / clusters above a certain threshold\r
+   }\r
+\r
    if(iPhysicsSelection)pwg4spec->SelectCollisionCandidates();\r
 \r
    mgr->AddTask(pwg4spec);\r