PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx

   1 // $Id$
   2 //
   3 // Emcal jet finder task.
   4 //
   5 // Authors: C.Loizides, S.Aiola
   6
   7 #include <vector>
   8 #include "AliEmcalJetTask.h"
   9
  10 #include <TChain.h>
  11 #include <TClonesArray.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14
  15 #include "AliAnalysisManager.h"
  16 #include "AliCentrality.h"
  17 #include "AliEMCALGeometry.h"
  18 #include "AliESDEvent.h"
  19 #include "AliEmcalJet.h"
  20 #include "AliEmcalParticle.h"
  21 #include "AliFJWrapper.h"
  22 #include "AliMCEvent.h"
  23 #include "AliVCluster.h"
  24 #include "AliVEvent.h"
  25 #include "AliVParticle.h"
  26
  27 ClassImp(AliEmcalJetTask)
  28
  29 //________________________________________________________________________
  30 AliEmcalJetTask::AliEmcalJetTask() :
  31   AliAnalysisTaskSE("AliEmcalJetTask"),
  32   fTracksName("Tracks"),
  33   fCaloName("CaloClusters"),
  34   fJetsName("Jets"),
  35   fAlgo(1),
  36   fRadius(0.4),
  37   fType(0),
  38   fMinJetTrackPt(0.15),
  39   fMinJetClusPt(0.15),
  40   fPhiMin(0),
  41   fPhiMax(TMath::TwoPi()),
  42   fEtaMin(-1),
  43   fEtaMax(1),
  44   fMinJetArea(0.01),
  45   fMinJetPt(1.0),
  46   fIsInit(0),
  47   fIsMcPart(0),
  48   fIsEmcPart(0),
  49   fJets(0),
  50   fEvent(0),
  51   fTracks(0),
  52   fClus(0)
  53 {
  54   // Default constructor.
  55 }
  56
  57 //________________________________________________________________________
  58 AliEmcalJetTask::AliEmcalJetTask(const char *name) :
  59   AliAnalysisTaskSE(name),
  60   fTracksName("Tracks"),
  61   fCaloName("CaloClusters"),
  62   fJetsName("Jets"),
  63   fAlgo(1),
  64   fRadius(0.4),
  65   fType(0),
  66   fMinJetTrackPt(0.15),
  67   fMinJetClusPt(0.15),
  68   fPhiMin(0),
  69   fPhiMax(TMath::TwoPi()),
  70   fEtaMin(-1),
  71   fEtaMax(1),
  72   fMinJetArea(0.01),
  73   fMinJetPt(1.0),
  74   fIsInit(0),
  75   fIsMcPart(0),
  76   fIsEmcPart(0),
  77   fJets(0),
  78   fEvent(0),
  79   fTracks(0),
  80   fClus(0)
  81 {
  82   // Standard constructor.
  83
  84   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
  85 }
  86
  87 //________________________________________________________________________
  88 AliEmcalJetTask::~AliEmcalJetTask()
  89 {
  90   // Destructor
  91 }
  92
  93 //________________________________________________________________________
  94 void AliEmcalJetTask::UserCreateOutputObjects()
  95 {
  96   // Create user objects.
  97
  98   fJets = new TClonesArray("AliEmcalJet");
  99   fJets->SetName(fJetsName);
 100 }
 101
 102 //________________________________________________________________________
 103 void AliEmcalJetTask::UserExec(Option_t *)
 104 {
 105   // Main loop, called for each event.
 106
 107   if (!fIsInit) {
 108     if (!DoInit())
 109       return;
 110     fIsInit = kTRUE;
 111   }
 112
 113   FindJets();
 114 }
 115
 116 //________________________________________________________________________
 117 void AliEmcalJetTask::Terminate(Option_t *)
 118 {
 119   // Called once at the end of the analysis.
 120 }
 121
 122 //________________________________________________________________________
 123 void AliEmcalJetTask::FindJets()
 124 {
 125   // Find jets.
 126
 127   if (!fTracks && !fClus)
 128     return;
 129
 130   TString name("kt");
 131   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
 132   if (fAlgo>=1) {
 133     name  = "antikt";
 134     jalgo = fastjet::antikt_algorithm;
 135   }
 136
 137   // setup fj wrapper
 138   AliFJWrapper fjw(name, name);
 139   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
 140   fjw.SetR(fRadius);
 141   fjw.SetAlgorithm(jalgo);
 142   fjw.SetMaxRap(1);
 143   fjw.Clear();
 144
 145   // get primary vertex
 146   Double_t vertex[3] = {0, 0, 0};
 147   fEvent->GetPrimaryVertex()->GetXYZ(vertex);
 148
 149   if ((fIsMcPart || fType == 0 || fType == 1) && fTracks) {
 150     const Int_t Ntracks = fTracks->GetEntries();
 151     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 152       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
 153       if (!t)
 154         continue;
 155       if (fIsMcPart && fType == 1 && t->Charge() == 0)
 156         continue;
 157       if (fIsMcPart && fType == 2 && t->Charge() != 0)
 158         continue;
 159       if (t->Pt() < fMinJetTrackPt)
 160         continue;
 161       Double_t eta = t->Eta();
 162       Double_t phi = t->Phi();
 163       if ((eta<fEtaMin) && (eta>fEtaMax) &&
 164           (phi<fPhiMin) && (phi<fPhiMax))
 165         continue;
 166
 167       // offset of 100 for consistency with cluster ids
 168       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
 169     }
 170   }
 171
 172   if ((fType == 0 || fType == 2) && fClus) {
 173     const Int_t Nclus = fClus->GetEntries();
 174     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 175       AliVCluster *c = 0;
 176       Double_t cEta=0,cPhi=0,cPt=0;
 177       Double_t cPx=0,cPy=0,cPz=0;
 178       if (fIsEmcPart) {
 179         AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
 180         if (!ep)
 181           continue;
 182         c = ep->GetCluster();
 183         if (!c)
 184             continue;
 185         cEta = ep->Eta();
 186         cPhi = ep->Phi();
 187         cPt  = ep->Pt();
 188         cPx  = ep->Px();
 189         cPy  = ep->Py();
 190         cPz  = ep->Pz();
 191       } else {
 192         c = static_cast<AliVCluster*>(fClus->At(iClus));
 193         if (!c)
 194           continue;
 195         TLorentzVector nP;
 196         c->GetMomentum(nP, vertex);
 197         cEta = nP.Eta();
 198         cPhi = nP.Phi();
 199         cPt  = nP.Pt();
 200         cPx  = nP.Px();
 201         cPy  = nP.Py();
 202         cPz  = nP.Pz();
 203       }
 204       if (!c->IsEMCAL())
 205         continue;
 206       if (cPt < fMinJetClusPt)
 207         continue;
 208       if ((cEta<fEtaMin) && (cEta>fEtaMax) &&
 209           (cPhi<fPhiMin) && (cPhi<fPhiMax))
 210         continue;
 211       // offset of 100 to skip ghost particles uid = -1
 212       fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
 213     }
 214   }
 215
 216   // run jet finder
 217   fjw.Run();
 218
 219   // get geometry
 220   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 221   if (!geom) {
 222     AliFatal(Form("%s: Can not create geometry", GetName()));
 223     return;
 224   }
 225
 226   // loop over fastjet jets
 227   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
 228   for (UInt_t ij=0, jetCount=0; ij<jets_incl.size(); ++ij) {
 229     if (jets_incl[ij].perp()<fMinJetPt)
 230       continue;
 231     if (fjw.GetJetArea(ij)<fMinJetArea)
 232       continue;
 233
 234     AliEmcalJet *jet = new ((*fJets)[jetCount])
 235       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
 236
 237     // loop over constituents
 238     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
 239     jet->SetNumberOfTracks(constituents.size());
 240     jet->SetNumberOfClusters(constituents.size());
 241
 242     Int_t nt            = 0;
 243     Int_t nc            = 0;
 244     Double_t neutralE   = 0;
 245     Double_t maxCh      = 0;
 246     Double_t maxNe      = 0;
 247     Int_t gall          = 0;
 248     Int_t gemc          = 0;
 249     Int_t cemc          = 0;
 250     Int_t ncharged      = 0;
 251     Int_t nneutral      = 0;
 252     Double_t mcpt       = 0;
 253     Double_t emcpt      = 0;
 254
 255     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
 256       Int_t uid = constituents[ic].user_index();
 257
 258       if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
 259         ++gall;
 260         Double_t gphi = constituents[ic].phi();
 261         if (gphi<0)
 262           gphi += TMath::TwoPi();
 263         gphi *= TMath::RadToDeg();
 264         Double_t geta = constituents[ic].eta();
 265         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
 266             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
 267           ++gemc;
 268       } else if ((uid > 0) && fTracks) { // track constituent
 269         Int_t tid = uid - 100;
 270         AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
 271         if (!t)
 272           continue;
 273         Double_t cEta = t->Eta();
 274         Double_t cPhi = t->Phi();
 275         Double_t cPt  = t->Pt();
 276         Double_t cP   = t->P();
 277         if (t->Charge() == 0) {
 278           neutralE += cP;
 279           ++nneutral;
 280           if (cPt > maxNe)
 281             maxNe = cPt;
 282         } else {
 283           ++ncharged;
 284           if (cPt > maxCh)
 285             maxCh = cPt;
 286         }
 287
 288         if (fIsMcPart || t->GetLabel() == 100) // check if MC particle
 289           mcpt += cPt;
 290
 291         if (cPhi<0)
 292           cPhi += TMath::TwoPi();
 293         cPhi *= TMath::RadToDeg();
 294         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
 295             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
 296           emcpt += cPt;
 297           ++cemc;
 298         }
 299
 300         jet->AddTrackAt(tid, nt);
 301         ++nt;
 302       } else if (fClus) { // cluster constituent
 303         Int_t cid = -uid - 100;
 304         AliVCluster *c = 0;
 305         Double_t cEta=0,cPhi=0,cPt=0,cP=0;
 306         if (fIsEmcPart) {
 307           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
 308           if (!ep)
 309             continue;
 310           c = ep->GetCluster();
 311           if (!c)
 312             continue;
 313           cEta = ep->Eta();
 314           cPhi = ep->Phi();
 315           cPt  = ep->Pt();
 316           cP   = ep->P();
 317         } else {
 318           c = static_cast<AliVCluster*>(fClus->At(cid));
 319           if (!c)
 320             continue;
 321           TLorentzVector nP;
 322           c->GetMomentum(nP, vertex);
 323           cEta = nP.Eta();
 324           cPhi = nP.Phi();
 325           cPt  = nP.Pt();
 326           cP   = nP.P();
 327         }
 328
 329         neutralE += cP;
 330         if (cPt > maxNe)
 331           maxNe = cPt;
 332
 333         if (c->Chi2() == 100) // MC particle
 334           mcpt += cPt;
 335
 336         if (cPhi<0)
 337           cPhi += TMath::TwoPi();
 338         cPhi *= TMath::RadToDeg();
 339         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
 340             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
 341           emcpt += cPt;
 342           ++cemc;
 343         }
 344
 345         jet->AddClusterAt(cid, nc);
 346         ++nc;
 347         ++nneutral;
 348       } else {
 349         AliError(Form("%s: No logical way to end up here.", GetName()));
 350         continue;
 351       }
 352     } /* end of constituent loop */
 353
 354     jet->SetNumberOfTracks(nt);
 355     jet->SetNumberOfClusters(nc);
 356     jet->SortConstituents();
 357     jet->SetMaxNeutralPt(maxNe);
 358     jet->SetMaxChargedPt(maxCh);
 359     jet->SetNEF(neutralE / jet->E());
 360     fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
 361     jet->SetArea(area.perp());
 362     jet->SetAreaEta(area.eta());
 363     jet->SetAreaPhi(area.phi());
 364     jet->SetNumberOfCharged(ncharged);
 365     jet->SetNumberOfNeutrals(nneutral);
 366     jet->SetMCPt(mcpt);
 367     jet->SetNEmc(cemc);
 368     jet->SetPtEmc(emcpt);
 369
 370     if (gall > 0)
 371       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
 372     else
 373       jet->SetAreaEmc(-1);
 374     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
 375         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
 376         (jet->Eta() > geom->GetArm1EtaMin()) &&
 377         (jet->Eta() < geom->GetArm1EtaMax()))
 378       jet->SetAxisInEmcal(kTRUE);
 379     jetCount++;
 380   }
 381   fJets->Sort();
 382 }
 383
 384 //________________________________________________________________________
 385 Bool_t AliEmcalJetTask::DoInit()
 386 {
 387   // Init. Return true if successful.
 388
 389   // get input collections
 390   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 391
 392   // get the event
 393   fEvent = InputEvent();
 394   if (!fEvent) {
 395     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
 396     return 0;
 397   }
 398
 399   // add jets to event if not yet there
 400   if (!(fEvent->FindListObject(fJetsName)))
 401     fEvent->AddObject(fJets);
 402   else {
 403     AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
 404     return 0;
 405   }
 406
 407   if (fTracksName == "Tracks")
 408     am->LoadBranch("Tracks");
 409   if (!fTracks && !fTracksName.IsNull()) {
 410     fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
 411     if (!fTracks) {
 412       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
 413       return 0;
 414     }
 415   }
 416   if (fTracks) {
 417     TClass cls(fTracks->GetClass()->GetName());
 418     if (cls.InheritsFrom("AliMCParticle"))
 419       fIsMcPart = 1;
 420   }
 421
 422   if (fCaloName == "CaloClusters")
 423     am->LoadBranch("CaloClusters");
 424   if (!fClus && !fCaloName.IsNull()) {
 425     fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
 426     if (!fClus) {
 427       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
 428       return 0;
 429     }
 430   }
 431   if (fClus) {
 432     TClass cls(fClus->GetClass()->GetName());
 433     if (cls.InheritsFrom("AliEmcalParticle"))
 434       fIsEmcPart = 1;
 435   }
 436
 437   return 1;
 438 }