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