PWGGA/EMCALJetTasks/AliEmcalJetTask.cxx

   1 // $Id$
   2 //
   3 // Emcal jet finder task.
   4 //
   5 // Authors: C.Loizides, S.Aiola
   6
   7 #include <TChain.h>
   8 #include <TClonesArray.h>
   9 #include <TList.h>
  10 #include <TLorentzVector.h>
  11
  12 #include "AliAnalysisManager.h"
  13 #include "AliCentrality.h"
  14 #include "AliEMCALGeometry.h"
  15 #include "AliEmcalJet.h"
  16 #include "AliFJWrapper.h"
  17 #include "AliVCluster.h"
  18 #include "AliVParticle.h"
  19 #include "AliVEvent.h"
  20 #include "AliESDEvent.h"
  21 #include "AliMCEvent.h"
  22
  23 #include "AliEmcalJetTask.h"
  24
  25 ClassImp(AliEmcalJetTask)
  26
  27 //________________________________________________________________________
  28 AliEmcalJetTask::AliEmcalJetTask() :
  29   AliAnalysisTaskSE("AliEmcalJetTask"),
  30   fTracksName("Tracks"),
  31   fCaloName("CaloClusters"),
  32   fJetsName("Jets"),
  33   fAlgo(1),
  34   fRadius(0.4),
  35   fType(0),
  36   fMinJetTrackPt(0.15),
  37   fMinJetClusPt(0.15),
  38   fMinJetArea(0.01),
  39   fMinJetPt(1.0),
  40   fJets(0),
  41   fEvent(0)
  42 {
  43   // Default constructor.
  44 }
  45
  46 //________________________________________________________________________
  47 AliEmcalJetTask::AliEmcalJetTask(const char *name) :
  48   AliAnalysisTaskSE(name),
  49   fTracksName("Tracks"),
  50   fCaloName("CaloClusters"),
  51   fJetsName("Jets"),
  52   fAlgo(1),
  53   fRadius(0.4),
  54   fType(0),
  55   fMinJetTrackPt(0.15),
  56   fMinJetClusPt(0.15),
  57   fMinJetArea(0.01),
  58   fMinJetPt(1.0),
  59   fJets(0),
  60   fEvent(0)
  61 {
  62   // Standard constructor.
  63
  64   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
  65 }
  66
  67 //________________________________________________________________________
  68 AliEmcalJetTask::~AliEmcalJetTask()
  69 {
  70   // Destructor
  71 }
  72
  73 //________________________________________________________________________
  74 void AliEmcalJetTask::UserCreateOutputObjects()
  75 {
  76   // Create user objects.
  77
  78   fJets = new TClonesArray("AliEmcalJet");
  79   fJets->SetName(fJetsName);
  80 }
  81
  82 //_____________________________________________________
  83 TString AliEmcalJetTask::GetBeamType()
  84 {
  85   // Get beam type : pp-AA-pA
  86   // ESDs have it directly, AODs get it from hardcoded run number ranges
  87
  88   TString beamType;
  89
  90   AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fEvent);
  91   if (esd) {
  92     const AliESDRun *run = esd->GetESDRun();
  93     beamType = run->GetBeamType();
  94   }
  95   else
  96   {
  97     Int_t runNumber = fEvent->GetRunNumber();
  98     if ((runNumber >= 136851 && runNumber <= 139517) ||  // LHC10h
  99         (runNumber >= 166529 && runNumber <= 170593))    // LHC11h
 100     {
 101       beamType = "A-A";
 102     }
 103     else
 104     {
 105       beamType = "p-p";
 106     }
 107   }
 108
 109   return beamType;
 110 }
 111
 112
 113 //________________________________________________________________________
 114 void AliEmcalJetTask::UserExec(Option_t *)
 115 {
 116   // Main loop, called for each event.
 117
 118   // get the event
 119   fEvent = InputEvent();
 120
 121   if (!fEvent) {
 122     AliError("Could not retrieve event! Returning");
 123     return;
 124   }
 125
 126   // add jets to event if not yet there
 127   if (!(fEvent->FindListObject(fJetsName)))
 128     fEvent->AddObject(fJets);
 129
 130   // delete jet output
 131   fJets->Delete();
 132
 133   // get input collections
 134   TClonesArray *tracks = 0;
 135   TClonesArray *clus   = 0;
 136   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 137
 138   if ((fType == 0) || (fType == 1)) {
 139     if (fTracksName == "Tracks")
 140       am->LoadBranch("Tracks");
 141     tracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
 142     if (!tracks) {
 143       AliError(Form("Pointer to tracks %s == 0", fTracksName.Data()));
 144       return;
 145     }
 146   }
 147   if ((fType == 0) || (fType == 2)) {
 148     if (fCaloName == "CaloClusters")
 149       am->LoadBranch("CaloClusters");
 150     clus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
 151     if (!clus) {
 152       AliError(Form("Pointer to clus %s == 0", fCaloName.Data()));
 153       return;
 154     }
 155   }
 156
 157
 158   // get centrality
 159   Double_t cent = 99;
 160
 161   if (GetBeamType() == "A-A") {
 162     AliCentrality *centrality = InputEvent()->GetCentrality();
 163
 164     if (centrality)
 165       cent = centrality->GetCentralityPercentile("V0M");
 166     else
 167       cent = 99; // probably pp data
 168
 169     if (cent < 0) {
 170       AliWarning(Form("Centrality negative: %f, assuming 99", cent));
 171       cent = 99;
 172     }
 173   }
 174
 175   FindJets(tracks, clus, fAlgo, fRadius, cent);
 176 }
 177
 178 //________________________________________________________________________
 179 void AliEmcalJetTask::Terminate(Option_t *)
 180 {
 181   // Called once at the end of the analysis.
 182 }
 183
 184 //________________________________________________________________________
 185 void AliEmcalJetTask::FindJets(TObjArray *tracks, TObjArray *clus, Int_t algo, Double_t radius, Float_t /*cent*/)
 186 {
 187   // Find jets.
 188
 189   if (!tracks && !clus)
 190     return;
 191
 192   TString name("kt");
 193   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
 194   if (algo>=1) {
 195     name  = "antikt";
 196     jalgo = fastjet::antikt_algorithm;
 197   }
 198
 199   AliFJWrapper fjw(name, name);
 200   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
 201   fjw.SetR(radius);
 202   fjw.SetAlgorithm(jalgo);
 203   fjw.SetMaxRap(0.9);
 204   fjw.Clear();
 205
 206   if (tracks) {
 207     const Int_t Ntracks = tracks->GetEntries();
 208     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 209       AliVParticle *t = dynamic_cast<AliVParticle*>(tracks->At(iTracks));
 210       if (!t)
 211         continue;
 212
 213       if (t->Pt() < fMinJetTrackPt)
 214         continue;
 215       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);  // offset of 100 for consistency with cluster ids
 216     }
 217   }
 218
 219   if (clus) {
 220     Double_t vertex[3] = {0, 0, 0};
 221     fEvent->GetPrimaryVertex()->GetXYZ(vertex);
 222     const Int_t Nclus = clus->GetEntries();
 223     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 224       AliVCluster *c = dynamic_cast<AliVCluster*>(clus->At(iClus));
 225       if (!c)
 226         continue;
 227       if (!c->IsEMCAL())
 228         continue;
 229       TLorentzVector nPart;
 230       c->GetMomentum(nPart, vertex);
 231       if (nPart.Pt() < fMinJetClusPt)
 232         continue;
 233       fjw.AddInputVector(nPart.Px(), nPart.Py(), nPart.Pz(), nPart.P(), -iClus - 100);  // offset of 100 to skip ghost particles uid = -1
 234     }
 235   }
 236
 237   // run jet finder
 238   fjw.Run();
 239
 240   // get geometry
 241   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 242   if (!geom) {
 243     AliFatal("Can not create geometry");
 244     return;
 245   }
 246
 247   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
 248   for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
 249     if (jets_incl[ij].perp()<fMinJetPt)
 250       continue;
 251     if (fjw.GetJetArea(ij)<fMinJetArea)
 252       continue;
 253     AliEmcalJet *jet = new ((*fJets)[jetCount])
 254       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
 255     Double_t vertex[3] = {0, 0, 0};
 256     fEvent->GetPrimaryVertex()->GetXYZ(vertex);
 257     vector<fastjet::PseudoJet> constituents = fjw.GetJetConstituents(ij);
 258     Int_t nt            = 0;
 259     Int_t nc            = 0;
 260     Double_t neutralE   = 0;
 261     Double_t maxCh      = 0;
 262     Double_t maxNe      = 0;
 263     Int_t gall          = 0;
 264     Int_t gemc          = 0;
 265     Int_t ncharged      = 0;
 266     Int_t nneutral      = 0;
 267     Double_t MCpt       = 0;
 268
 269     jet->SetNumberOfTracks(constituents.size());
 270     jet->SetNumberOfClusters(constituents.size());
 271     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
 272       Int_t uid = constituents[ic].user_index();
 273
 274       if ((uid == -1) && (constituents[ic].kt2() < 1e-25)) { //ghost particle
 275         ++gall;
 276         Double_t gphi = constituents[ic].phi() * TMath::RadToDeg();
 277         Double_t geta = constituents[ic].eta();
 278         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
 279             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
 280           ++gemc;
 281         continue;
 282       } else if (uid > 0) {
 283         Int_t tid = uid - 100;
 284         AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
 285         if (t) {
 286           if (t->Charge() == 0) {
 287             neutralE += t->P();
 288             ++nneutral;
 289             if (t->Pt() > maxNe)
 290               maxNe = t->Pt();
 291           } else {
 292             ++ncharged;
 293             if (t->Pt() > maxCh)
 294               maxCh = t->Pt();
 295           }
 296
 297           if (t->InheritsFrom("AliMCParticle") || t->GetLabel() == 100) // MC particle
 298             MCpt += t->Pt();
 299
 300           jet->AddTrackAt(tid, nt);
 301           ++nt;
 302         }
 303       } else {
 304         Int_t cid = -uid - 100;
 305         AliVCluster *c = static_cast<AliVCluster*>(clus->At(cid));
 306         if (c) {
 307           TLorentzVector nP;
 308           c->GetMomentum(nP, vertex);
 309           neutralE += nP.P();
 310           if (nP.Pt() > maxNe)
 311             maxNe = nP.Pt();
 312
 313           if (c->Chi2() == 100) // MC particle
 314             MCpt += nP.Pt();
 315
 316           jet->AddClusterAt(cid, nc);
 317           ++nc;
 318           ++nneutral;
 319         }
 320       }
 321     }
 322     jet->SetNumberOfTracks(nt);
 323     jet->SetNumberOfClusters(nc);
 324     jet->SortConstituents();
 325     jet->SetMaxNeutralPt(maxNe);
 326     jet->SetMaxChargedPt(maxCh);
 327     jet->SetNEF(neutralE / jet->E());
 328     jet->SetArea(fjw.GetJetArea(ij));
 329     jet->SetNumberOfCharged(ncharged);
 330     jet->SetNumberOfNeutrals(nneutral);
 331     jet->SetMCPt(MCpt);
 332     if (gall > 0)
 333       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
 334     else
 335       jet->SetAreaEmc(-1);
 336     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
 337         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
 338         (jet->Eta() > geom->GetArm1EtaMin()) &&
 339         (jet->Eta() < geom->GetArm1EtaMax()))
 340       jet->SetAxisInEmcal(kTRUE);
 341     jetCount++;
 342   }
 343 }