PWGGA/EMCALJetTasks/AliEsdJetTask.cxx

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