integrate mc jets
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 May 2012 19:38:00 +0000 (19:38 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 May 2012 19:38:00 +0000 (19:38 +0000)
PWGGA/EMCALJetTasks/AliEmcalJetTask.cxx
PWGGA/EMCALJetTasks/AliEmcalJetTask.h

index c5f86d3..2765c53 100644 (file)
@@ -2,23 +2,24 @@
 //
 // Emcal jet finder task.
 //
-// Author: C.Loizides
+// Authors: C.Loizides, S.Aiola
 
-#include "AliEmcalJetTask.h"
 #include <TChain.h>
 #include <TClonesArray.h>
-#include <TH1F.h>
-#include <TH2F.h>
 #include <TList.h>
 #include <TLorentzVector.h>
-#include <TParticle.h>
+
 #include "AliAnalysisManager.h"
 #include "AliCentrality.h"
 #include "AliEMCALGeometry.h"
 #include "AliEmcalJet.h"
 #include "AliFJWrapper.h"
 #include "AliVCluster.h"
-#include "AliVTrack.h"
+#include "AliVParticle.h"
+#include "AliVEvent.h"
+#include "AliMCEvent.h"
+
+#include "AliEmcalJetTask.h"
 
 ClassImp(AliEmcalJetTask)
 
@@ -28,6 +29,7 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
+  fMC(kFALSE),
   fAlgo(1),
   fRadius(0.4),
   fType(0),
@@ -35,7 +37,8 @@ AliEmcalJetTask::AliEmcalJetTask() :
   fMinJetClusPt(0.15),
   fMinJetArea(0.01),
   fMinJetPt(1.0),
-  fJets(0)
+  fJets(0),
+  fEvent(0)
 {
   // Default constructor.
 }
@@ -46,6 +49,7 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fTracksName("Tracks"),
   fCaloName("CaloClusters"),
   fJetsName("Jets"),
+  fMC(kFALSE),
   fAlgo(1),
   fRadius(0.4),
   fType(0),
@@ -53,7 +57,8 @@ AliEmcalJetTask::AliEmcalJetTask(const char *name) :
   fMinJetClusPt(0.15),
   fMinJetArea(0.01),
   fMinJetPt(1.0),
-  fJets(0)
+  fJets(0),
+  fEvent(0)
 {
   // Standard constructor.
 
@@ -80,32 +85,42 @@ void AliEmcalJetTask::UserExec(Option_t *)
 {
   // Main loop, called for each event.
 
+  // get the event
+  if (fMC) 
+    fEvent = MCEvent();
+  else
+    fEvent = InputEvent();
+
   // add jets to event if not yet there
-  if (!(InputEvent()->FindListObject(fJetsName)))
-    InputEvent()->AddObject(fJets);
+  if (!(fEvent->FindListObject(fJetsName)))
+    fEvent->AddObject(fJets);
 
   // delete jet output
   fJets->Delete();
 
+  if (!fEvent) {
+    AliError(Form("Could not get the event! fMC = %d", fMC));
+    return;
+  }
+
   // get input collections
   TClonesArray *tracks = 0;
   TClonesArray *clus   = 0;
-  TList *l = InputEvent()->GetList();
   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 
-  if ((fType==0)||(fType==1)) {
+  if ((fType == 0) || (fType == 1)) {
     if (fTracksName == "Tracks")
       am->LoadBranch("Tracks");
-    tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
+    tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
     if (!tracks) {
       AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
       return;
     }
   }
-  if ((fType==0)||(fType==2)) {
+  if ((fType == 0) || (fType == 2)) {
     if (fCaloName == "CaloClusters")
       am->LoadBranch("CaloClusters");
-    clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
+    clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
     if (!clus) {
       AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
       return;
@@ -114,12 +129,12 @@ void AliEmcalJetTask::UserExec(Option_t *)
 
   // get centrality
   Float_t cent = -1; 
-  AliCentrality *centrality = InputEvent()->GetCentrality();
+  AliCentrality *centrality = fEvent->GetCentrality();
   if (centrality)
     cent = centrality->GetCentralityPercentile("V0M");
   else
     cent=99; // probably pp data
-  if (cent<0) {
+  if (cent < 0) {
     AliError(Form("Centrality negative: %f", cent));
     return;
   }
@@ -158,18 +173,23 @@ void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, D
   if (tracks) {
     const Int_t Ntracks = tracks->GetEntries();
     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
-      AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
+      AliVParticle *t = static_cast<AliVParticle*>(tracks->At(iTracks));
       if (!t)
         continue;
       if (t->Pt()<fMinJetTrackPt) 
         continue;
-      fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
+
+      Int_t index = 1;
+      if(fMC && t->Charge() == 0)
+       index =- 1;
+
+      fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), (iTracks + 123) * index);
     }
   }
 
   if (clus) {
     Double_t vertex[3] = {0, 0, 0};
-    InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+    fEvent->GetPrimaryVertex()->GetXYZ(vertex);
     const Int_t Nclus = clus->GetEntries();
     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
       AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
@@ -180,9 +200,9 @@ void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, D
       TLorentzVector nPart;
       c->GetMomentum(nPart, vertex);
       Double_t et = nPart.Pt();
-      if (et<fMinJetClusPt) 
+      if (et < fMinJetClusPt) 
         continue;
-      fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus-123);
+      fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 123);
     }
   }
 
@@ -197,7 +217,7 @@ void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, D
   }
   
   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
-  for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
+  for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
     if (jets_incl[ij].perp()<fMinJetPt) 
       continue;
     if (fjw.GetJetArea(ij)<fMinJetArea)
@@ -205,48 +225,53 @@ void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, D
     AliEmcalJet *jet = new ((*fJets)[jetCount]) 
       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
     Double_t vertex[3] = {0, 0, 0};
-    InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
+    fEvent->GetPrimaryVertex()->GetXYZ(vertex);
     vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
-    Int_t nt = 0;
-    Int_t nc = 0;
+    Int_t nt            = 0;
+    Int_t nc            = 0;
     Double_t neutralE   = 0;
     Double_t maxTrack   = 0;
     Double_t maxCluster = 0;
     Int_t gall          = 0;
     Int_t gemc          = 0;
+
     jet->SetNumberOfTracks(constituents.size());
     jet->SetNumberOfClusters(constituents.size());
-    for(UInt_t ic=0; ic<constituents.size(); ++ic) {
+
+    for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
       Int_t uid = constituents[ic].user_index();
-      if ((uid==-1) && (constituents[ic].kt2()<1e-25)) { //ghost particle
+
+      if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
         ++gall;
-        Double_t gphi = constituents[ic].phi()*TMath::RadToDeg();
+        Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
         Double_t geta = constituents[ic].eta();
-        if ((gphi>geom->GetArm1PhiMin()) && (gphi<geom->GetArm1PhiMax()) &&
-            (geta>geom->GetArm1EtaMin()) && (geta<geom->GetArm1EtaMax()))
+        if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
+            (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
           ++gemc;
         continue;
-      }
-      if (uid>=0){
-        AliVTrack *t = static_cast<AliVTrack*>(tracks->At(uid));
+      }        else if (fMC || uid >= 0) {
+       Int_t tid = TMath::Abs(uid) - 123;
+        AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
         if (t) {
-          if (t->Pt()>maxTrack)
-            maxTrack=t->Pt();
-          jet->AddTrackAt(uid, nt);
+         if (uid >= 0)
+           neutralE += t->P();
+          if (t->Pt() > maxTrack)
+            maxTrack = t->Pt();
+          jet->AddTrackAt(tid, nt);
           nt++;
         }
       } else {
-        Int_t cid = -uid - 123;
-        AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
-        if (c) {
-          TLorentzVector nP;
-          c->GetMomentum(nP, vertex);
-          neutralE += nP.P();
-          if (nP.Pt()>maxCluster)
-            maxCluster=nP.Pt();
-          jet->AddClusterAt(cid,nc);
-          nc++;
-        }
+       Int_t cid = TMath::Abs(uid) - 123;
+       AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
+       if (c) {
+         TLorentzVector nP;
+         c->GetMomentum(nP, vertex);
+         neutralE += nP.P();
+         if (nP.Pt() > maxCluster)
+           maxCluster = nP.Pt();
+         jet->AddClusterAt(cid, nc);
+         nc++;
+       }
       }
     }
     jet->SetNumberOfTracks(nt);
@@ -254,14 +279,16 @@ void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, D
     jet->SortConstituents();
     jet->SetMaxTrackPt(maxTrack);
     jet->SetMaxClusterPt(maxCluster);
-    jet->SetNEF(neutralE/jet->E());
+    jet->SetNEF(neutralE / jet->E());
     jet->SetArea(fjw.GetJetArea(ij));
-    if (gall>0)
-      jet->SetAreaEmc(fjw.GetJetArea(ij)*gemc/gall);
+    if (gall > 0)
+      jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
     else 
       jet->SetAreaEmc(-1);
-    if ((jet->Phi()>geom->GetArm1PhiMin()*TMath::DegToRad()) && (jet->Phi()<geom->GetArm1PhiMax()*TMath::DegToRad()) &&
-        (jet->Eta()>geom->GetArm1EtaMin()) && (jet->Eta()<geom->GetArm1EtaMax()))
+    if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) && 
+        (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
+        (jet->Eta() > geom->GetArm1EtaMin()) && 
+        (jet->Eta() < geom->GetArm1EtaMax()))
       jet->SetAxisInEmcal(kTRUE);
     jetCount++;
   }
index 178de16..932cd87 100644 (file)
@@ -4,9 +4,7 @@
 // $Id$
 
 class TClonesArray;
-class TList;
-class TH1;
-class TH2;
+class AliVEvent;
 
 #include "AliAnalysisTaskSE.h"
 
@@ -16,27 +14,30 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   AliEmcalJetTask(const char *name);
   virtual ~AliEmcalJetTask();
 
-  void         UserCreateOutputObjects();
-  void         UserExec(Option_t *option);
-  void         Terminate(Option_t *option);
-
-  void         SetAlgo(Int_t a)                 { fAlgo          = a;  }
-  void         SetClusName(const char *n)       { fCaloName      = n;  }
-  void         SetJetsName(const char *n)       { fJetsName      = n;  }
-  void         SetMinJetArea(Double_t a)        { fMinJetArea    = a;  }
-  void         SetMinJetClusPt(Double_t min)    { fMinJetClusPt  = min;}
-  void         SetMinJetPt(Double_t j)          { fMinJetPt      = j;  }
-  void         SetMinJetTrackPt(Double_t min)   { fMinJetTrackPt = min;}
-  void         SetRadius(Double_t r)            { fRadius        = r;  }
-  void         SetTracksName(const char *n)     { fTracksName    = n;  }
-  void         SetType(Int_t t)                 { fType          = t;  }
+  void                   UserCreateOutputObjects();
+  void                   UserExec(Option_t *option);
+  void                   Terminate(Option_t *option);
+
+  void                   SetTracksName(const char *n)     { fTracksName    = n     ; }
+  void                   SetClusName(const char *n)       { fCaloName      = n     ; }
+  void                   SetJetsName(const char *n)       { fJetsName      = n     ; }
+  void                   SetMC(Bool_t mc = kTRUE)         { fMC            = mc    ; }
+  void                   SetAlgo(Int_t a)                 { fAlgo          = a     ; }
+  void                   SetRadius(Double_t r)            { fRadius        = r     ; }
+  void                   SetType(Int_t t)                 { fType          = t     ; }
+  void                   SetMinJetClusPt(Double_t min)    { fMinJetClusPt  = min   ; }
+  void                   SetMinJetTrackPt(Double_t min)   { fMinJetTrackPt = min   ; }
+  void                   SetMinJetArea(Double_t a)        { fMinJetArea    = a     ; }
+  void                   SetMinJetPt(Double_t j)          { fMinJetPt      = j     ; }
+
 
  protected:
-  void         FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/);
+  void                   FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/);
 
   TString                fTracksName;             // name of track collection
   TString                fCaloName;               // name of calo cluster collection
   TString                fJetsName;               // name of jet collection
+  Bool_t                 fMC;                     // true = MC analysis
   Int_t                  fAlgo;                   // algo (0==kt, 1==antikt)
   Double_t               fRadius;                 // jet radius
   Int_t                  fType;                   // jet type (0=all, 1=ch, 2=neutral)
@@ -45,11 +46,12 @@ class AliEmcalJetTask : public AliAnalysisTaskSE {
   Double_t               fMinJetArea;             // min area to keep jet in output
   Double_t               fMinJetPt;               // min jet pt to keep jet in output
   TClonesArray          *fJets;                   //!jet collection
+  AliVEvent             *fEvent;                  //!current event
 
  private:
   AliEmcalJetTask(const AliEmcalJetTask&);            // not implemented
   AliEmcalJetTask &operator=(const AliEmcalJetTask&); // not implemented
 
-  ClassDef(AliEmcalJetTask, 2) // Jet producing task
+  ClassDef(AliEmcalJetTask, 3) // Jet producing task
 };
 #endif