PWGGA/EMCALJetTasks/AliEmcalJetTask.cxx

   1 // $Id$
   2 //
   3 // Emcal jet finder task.
   4 //
   5 // Author: C.Loizides
   6
   7 #include "AliEmcalJetTask.h"
   8 #include <TChain.h>
   9 #include <TClonesArray.h>
  10 #include <TH1F.h>
  11 #include <TH2F.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TParticle.h>
  15 #include "AliAnalysisManager.h"
  16 #include "AliCentrality.h"
  17 #include "AliFJWrapper.h"
  18 #include "AliVCluster.h"
  19 #include "AliVTrack.h"
  20 #include "AliEmcalJet.h"
  21
  22 ClassImp(AliEmcalJetTask)
  23
  24 //________________________________________________________________________
  25 AliEmcalJetTask::AliEmcalJetTask() :
  26   AliAnalysisTaskSE("AliEmcalJetTask"),
  27   fTracksName("Tracks"),
  28   fCaloName("CaloClusters"),
  29   fJetsName("Jets"),
  30   fAlgo(1),
  31   fRadius(0.4),
  32   fType(0),
  33   fMinJetTrackPt(0.15),
  34   fMinJetClusPt(0.15),
  35   fJets(0)
  36 {
  37   // Default constructor.
  38 }
  39
  40 //________________________________________________________________________
  41 AliEmcalJetTask::AliEmcalJetTask(const char *name) :
  42   AliAnalysisTaskSE(name),
  43   fTracksName("Tracks"),
  44   fCaloName("CaloClusters"),
  45   fJetsName("Jets"),
  46   fAlgo(1),
  47   fRadius(0.4),
  48   fType(0),
  49   fMinJetTrackPt(0.15),
  50   fMinJetClusPt(0.15),
  51   fJets(0)
  52 {
  53   // Standard constructor.
  54
  55   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
  56 }
  57
  58 //________________________________________________________________________
  59 AliEmcalJetTask::~AliEmcalJetTask()
  60 {
  61   // Destructor
  62 }
  63
  64 //________________________________________________________________________
  65 void AliEmcalJetTask::UserCreateOutputObjects()
  66 {
  67   // Create user objects.
  68
  69   fJets = new TClonesArray("AliEmcalJet");
  70   fJets->SetName(fJetsName);
  71 }
  72
  73 //________________________________________________________________________
  74 void AliEmcalJetTask::UserExec(Option_t *)
  75 {
  76   // Main loop, called for each event.
  77
  78   // add jets to event if not yet there
  79   if (!(InputEvent()->FindListObject(fJetsName)))
  80     InputEvent()->AddObject(fJets);
  81
  82   // delete jet output
  83   fJets->Delete();
  84
  85
  86   // get input collections
  87   TClonesArray *tracks = 0;
  88   TClonesArray *clus   = 0;
  89   TList *l = InputEvent()->GetList();
  90   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
  91
  92   if ((fType==0)||(fType==1)) {
  93     if (fTracksName == "Tracks")
  94       am->LoadBranch("Tracks");
  95     tracks = dynamic_cast<TClonesArray*>(l->FindObject(fTracksName));
  96     if (!tracks) {
  97       AliError(Form("Pointer to tracks %s == 0", fTracksName.Data() ));
  98       return;
  99     }
 100   }
 101   if ((fType==0)||(fType==2)) {
 102     if (fCaloName == "CaloClusters")
 103       am->LoadBranch("CaloClusters");
 104     clus = dynamic_cast<TClonesArray*>(l->FindObject(fCaloName));
 105     if (!clus) {
 106       AliError(Form("Pointer to clus %s == 0", fCaloName.Data() ));
 107       return;
 108     }
 109   }
 110
 111   // get centrality
 112   Float_t cent = -1;
 113   AliCentrality *centrality = InputEvent()->GetCentrality() ;
 114   if (centrality)
 115     cent = centrality->GetCentralityPercentile("V0M");
 116   else
 117     cent=99; // probably pp data
 118   if (cent<0) {
 119     AliError(Form("Centrality negative: %f", cent));
 120     return;
 121   }
 122
 123   FindJets(tracks, clus, fAlgo, fRadius, cent);
 124 }
 125
 126 //________________________________________________________________________
 127 void AliEmcalJetTask::Terminate(Option_t *)
 128 {
 129   // Called once at the end of the analysis.
 130 }
 131
 132 //________________________________________________________________________
 133 void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
 134 {
 135   // Find jets.
 136
 137   if (!tracks && !clus)
 138     return;
 139
 140   TString name("kt");
 141   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
 142   if (algo>=1) {
 143     name  = "antikt";
 144     jalgo = fastjet::antikt_algorithm;
 145   }
 146
 147   AliFJWrapper fjw(name, name);
 148   fjw.SetR(radius);
 149   fjw.SetAlgorithm(jalgo);
 150   fjw.SetMaxRap(0.9);
 151   fjw.Clear();
 152
 153   if (tracks) {
 154     const Int_t Ntracks = tracks->GetEntries();
 155     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 156       AliVTrack *t = static_cast<AliVTrack*>(tracks->At(iTracks));
 157       if (!t)
 158         continue;
 159       if (t->Pt()<fMinJetTrackPt)
 160         continue;
 161       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks);
 162     }
 163   }
 164
 165   if (clus) {
 166     Double_t vertex[3] = {0, 0, 0};
 167     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 168     const Int_t Nclus = clus->GetEntries();
 169     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 170       AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
 171       if (!c)
 172         continue;
 173       if (!c->IsEMCAL())
 174         continue;
 175       TLorentzVector nPart;
 176       c->GetMomentum(nPart, vertex);
 177       Double_t et = nPart.Pt();
 178       if (et<fMinJetClusPt)
 179         continue;
 180       fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus-1);
 181     }
 182   }
 183
 184   // run jet finder
 185   fjw.Run();
 186
 187   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
 188   for(UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
 189     if (jets_incl[ij].perp()<1)
 190       continue;
 191     AliEmcalJet *jet = new ((*fJets)[jetCount])
 192       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
 193     jet->SetArea(fjw.GetJetArea(ij));
 194     Double_t vertex[3] = {0, 0, 0};
 195     InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
 196     vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
 197     Double_t neutralE = 0;Double_t maxTrack = 0;Double_t maxCluster=0;
 198     jet->SetNumberOfTracks(constituents.size());
 199     jet->SetNumberOfClusters(constituents.size());
 200     Int_t nt = 0;
 201     Int_t nc = 0;
 202     for(UInt_t ic=0; ic<constituents.size(); ++ic) {
 203       Int_t uid = constituents[ic].user_index();
 204       if (uid>=0){
 205         jet->AddTrackAt(uid, nt);
 206         AliVTrack *t = static_cast<AliVTrack*>(tracks->At(uid));
 207         if (t) {
 208           if (t->Pt()>maxTrack)
 209             maxTrack=t->Pt();
 210           nt++;
 211         }
 212       } else {
 213         jet->AddClusterAt(-(uid+1),nc);
 214         AliVCluster *c = static_cast<AliVCluster*>(clus->At(-(uid+1)));
 215         if (c) {
 216           TLorentzVector nP;
 217           c->GetMomentum(nP, vertex);
 218           neutralE += nP.P();
 219           if (nP.Pt()>maxCluster)
 220             maxCluster=nP.Pt();
 221         }
 222         nc++;
 223       }
 224     }
 225     jet->SetNumberOfTracks(nt);
 226     jet->SetNumberOfClusters(nc);
 227     jet->SetMaxTrackPt(maxTrack);
 228     jet->SetMaxClusterPt(maxCluster);
 229     jet->SetNEF(neutralE/jet->E());
 230     jetCount++;
 231   }
 232 }