PWGGA/EMCALJetTasks/AliEsdJetTask.cxx

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