PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx

   1 //
   2 // Emcal jet finder task.
   3 //
   4 // Authors: C.Loizides, S.Aiola, M. Verweij
   5
   6 #include <vector>
   7 #include "AliEmcalJetTask.h"
   8
   9 #include <TChain.h>
  10 #include <TClonesArray.h>
  11 #include <TList.h>
  12 #include <TLorentzVector.h>
  13 #include <TMath.h>
  14 #include <TRandom3.h>
  15
  16 #include "AliAnalysisManager.h"
  17 #include "AliCentrality.h"
  18 #include "AliEMCALGeometry.h"
  19 #include "AliESDEvent.h"
  20 #include "AliEmcalJet.h"
  21 #include "AliEmcalParticle.h"
  22 #include "AliFJWrapper.h"
  23 #include "AliMCEvent.h"
  24 #include "AliVCluster.h"
  25 #include "AliVEvent.h"
  26 #include "AliVParticle.h"
  27 #include "AliEmcalJetUtility.h"
  28
  29 using std::cout;
  30 using std::endl;
  31 using std::cerr;
  32
  33 ClassImp(AliEmcalJetTask)
  34
  35 //________________________________________________________________________
  36 AliEmcalJetTask::AliEmcalJetTask() :
  37   AliAnalysisTaskSE("AliEmcalJetTask"),
  38   fTracksName("Tracks"),
  39   fCaloName("CaloClusters"),
  40   fJetsName("Jets"),
  41   fJetType(kAKT|kFullJet|kRX1Jet),
  42   fMinLabelTracks(-kMaxInt),
  43   fMaxLabelTracks(kMaxInt),
  44   fMinLabelClusters(-kMaxInt),
  45   fMaxLabelClusters(kMaxInt),
  46   fMinMCLabel(0),
  47   fRadius(0.4),
  48   fMinJetTrackPt(0.15),
  49   fMinJetClusPt(0.15),
  50   fPhiMin(0),
  51   fPhiMax(TMath::TwoPi()),
  52   fEtaMin(-0.9),
  53   fEtaMax(+0.9),
  54   fMinJetArea(0.001),
  55   fMinJetPt(1.0),
  56   fJetPhiMin(-10),
  57   fJetPhiMax(+10),
  58   fJetEtaMin(-1),
  59   fJetEtaMax(+1),
  60   fGhostArea(0.005),
  61   fRecombScheme(fastjet::pt_scheme),
  62   fTrackEfficiency(1.),
  63   fMCFlag(0),
  64   fGeneratorIndex(-1),
  65   fUtilities(0),
  66   fLocked(0),
  67   fIsInit(0),
  68   fIsPSelSet(0),
  69   fIsEmcPart(0),
  70   fLegacyMode(kFALSE),
  71   fGeom(0),
  72   fJets(0),
  73   fEvent(0),
  74   fTracks(0),
  75   fClus(0),
  76   fFastJetWrapper("AliEmcalJetTask","AliEmcalJetTask")
  77 {
  78   // Default constructor.
  79
  80   fVertex[0] = 0.;
  81   fVertex[1] = 0.;
  82   fVertex[2] = 0.;
  83 }
  84
  85 //________________________________________________________________________
  86 AliEmcalJetTask::AliEmcalJetTask(const char *name) :
  87   AliAnalysisTaskSE(name),
  88   fTracksName("Tracks"),
  89   fCaloName("CaloClusters"),
  90   fJetsName("Jets"),
  91   fJetType(kAKT|kFullJet|kRX1Jet),
  92   fMinLabelTracks(-kMaxInt),
  93   fMaxLabelTracks(kMaxInt),
  94   fMinLabelClusters(-kMaxInt),
  95   fMaxLabelClusters(kMaxInt),
  96   fMinMCLabel(0),
  97   fRadius(0.4),
  98   fMinJetTrackPt(0.15),
  99   fMinJetClusPt(0.15),
 100   fPhiMin(0),
 101   fPhiMax(TMath::TwoPi()),
 102   fEtaMin(-0.9),
 103   fEtaMax(+0.9),
 104   fMinJetArea(0.001),
 105   fMinJetPt(1.0),
 106   fJetPhiMin(-10),
 107   fJetPhiMax(+10),
 108   fJetEtaMin(-1),
 109   fJetEtaMax(+1),
 110   fGhostArea(0.005),
 111   fRecombScheme(fastjet::pt_scheme),
 112   fTrackEfficiency(1.),
 113   fMCFlag(0),
 114   fGeneratorIndex(-1),
 115   fUtilities(0),
 116   fLocked(0),
 117   fIsInit(0),
 118   fIsPSelSet(0),
 119   fIsEmcPart(0),
 120   fLegacyMode(kFALSE),
 121   fGeom(0),
 122   fJets(0),
 123   fEvent(0),
 124   fTracks(0),
 125   fClus(0),
 126   fFastJetWrapper(name,name)
 127 {
 128   // Standard constructor.
 129
 130   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
 131
 132   fVertex[0] = 0.;
 133   fVertex[1] = 0.;
 134   fVertex[2] = 0.;
 135 }
 136
 137 //________________________________________________________________________
 138 AliEmcalJetTask::~AliEmcalJetTask()
 139 {
 140   // Destructor
 141 }
 142
 143 //________________________________________________________________________
 144 AliEmcalJetUtility* AliEmcalJetTask::AddUtility(AliEmcalJetUtility* utility)
 145 {
 146   // Addition of utilities.
 147
 148   if (!fUtilities) fUtilities = new TObjArray();
 149   if (fUtilities->FindObject(utility)) {
 150      Error("AddSupply", "Jet utility %s already connected.", utility->GetName());
 151      return utility;
 152   }
 153   fUtilities->Add(utility);
 154   utility->SetJetTask(this);
 155
 156   return utility;
 157 }
 158
 159 //________________________________________________________________________
 160 void AliEmcalJetTask::InitUtilities()
 161 {
 162   TIter next(fUtilities);
 163   AliEmcalJetUtility *utility = 0;
 164   while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Init();
 165 }
 166
 167 //________________________________________________________________________
 168 void AliEmcalJetTask::PrepareUtilities()
 169 {
 170   TIter next(fUtilities);
 171   AliEmcalJetUtility *utility = 0;
 172   while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Prepare(fFastJetWrapper);
 173 }
 174
 175 //________________________________________________________________________
 176 void AliEmcalJetTask::ExecuteUtilities(AliEmcalJet* jet, Int_t ij)
 177 {
 178   TIter next(fUtilities);
 179   AliEmcalJetUtility *utility = 0;
 180   while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->ProcessJet(jet, ij, fFastJetWrapper);
 181 }
 182
 183 //________________________________________________________________________
 184 void AliEmcalJetTask::TerminateUtilities()
 185 {
 186   TIter next(fUtilities);
 187   AliEmcalJetUtility *utility = 0;
 188   while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Terminate(fFastJetWrapper);
 189 }
 190
 191 //________________________________________________________________________
 192 void AliEmcalJetTask::UserCreateOutputObjects()
 193 {
 194   // Create user objects.
 195 }
 196
 197 //________________________________________________________________________
 198 void AliEmcalJetTask::UserExec(Option_t *)
 199 {
 200   // Main loop, called for each event.
 201   if (!fIsInit) {
 202     if (!DoInit())
 203       return;
 204     fIsInit = kTRUE;
 205   }
 206
 207   // clear the jet array (normally a null operation)
 208   fJets->Delete();
 209
 210   // get primary vertex
 211   if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(fVertex);
 212
 213   FindJets();
 214
 215   FillJetBranch();
 216 }
 217
 218 //________________________________________________________________________
 219 void AliEmcalJetTask::FindJets()
 220 {
 221   // Find jets.
 222
 223   if (!fTracks && !fClus){
 224     AliError("No tracks or clusters, returning.");
 225     return;
 226   }
 227
 228   fFastJetWrapper.Clear();
 229
 230   AliDebug(2,Form("Jet type = %d", fJetType));
 231
 232   if (fTracks) {
 233     const Int_t Ntracks = fTracks->GetEntries();
 234     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 235       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
 236
 237       if (!t) continue;
 238
 239       if (t->Pt() < fMinJetTrackPt) continue;
 240       Double_t eta = t->Eta();
 241       Double_t phi = t->Phi();
 242       if ((eta<fEtaMin) || (eta>fEtaMax) ||
 243           (phi<fPhiMin) || (phi>fPhiMax))
 244         continue;
 245
 246       if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
 247         AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
 248         continue;
 249       }
 250
 251       if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
 252         AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
 253         continue;
 254       }
 255
 256       Int_t lab = TMath::Abs(t->GetLabel());
 257       if (lab < fMinLabelTracks || lab > fMaxLabelTracks) {
 258         AliDebug(2,Form("Skipping track %d because label %d is not in range (%d, %d)", iTracks, lab, fMinLabelTracks, fMaxLabelTracks));
 259         continue;
 260       }
 261
 262       if ((t->GetFlag() & fMCFlag) != fMCFlag) {
 263         AliDebug(2,Form("Skipping track %d because it does not match the MC flags", iTracks));
 264         continue;
 265       }
 266
 267       if (fGeneratorIndex >= 0 && t->GetGeneratorIndex() != fGeneratorIndex) {
 268         AliDebug(2,Form("Skipping track %d because it does not match the MC generator index", iTracks));
 269         continue;
 270       }
 271
 272       // artificial inefficiency
 273       if (fTrackEfficiency < 1.) {
 274         Double_t rnd = gRandom->Rndm();
 275         if (fTrackEfficiency < rnd) {
 276           AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", iTracks));
 277           continue;
 278         }
 279       }
 280
 281       // offset of 100 for consistency with cluster ids
 282       AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, lab, t->Pt()));
 283       fFastJetWrapper.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
 284     }
 285   }
 286
 287   if (fClus) {
 288     const Int_t Nclus = fClus->GetEntries();
 289     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 290       AliVCluster *c = 0;
 291       Double_t cEta=0,cPhi=0,cPt=0;
 292       Double_t cPx=0,cPy=0,cPz=0;
 293       if (fIsEmcPart) {
 294         AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
 295         if (!ep) continue;
 296         c = ep->GetCluster();
 297         if (!c) continue;
 298
 299         cEta = ep->Eta();
 300         cPhi = ep->Phi();
 301         cPt  = ep->Pt();
 302         cPx  = ep->Px();
 303         cPy  = ep->Py();
 304         cPz  = ep->Pz();
 305       }
 306       else {
 307         c = static_cast<AliVCluster*>(fClus->At(iClus));
 308         if (!c) continue;
 309
 310         TLorentzVector nP;
 311         c->GetMomentum(nP, fVertex);
 312         cEta = nP.Eta();
 313         cPhi = nP.Phi();
 314         cPt  = nP.Pt();
 315         cPx  = nP.Px();
 316         cPy  = nP.Py();
 317         cPz  = nP.Pz();
 318       }
 319
 320       if (!c->IsEMCAL()) continue;
 321       if (cPt < fMinJetClusPt) continue;
 322       if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
 323           (cPhi<fPhiMin) || (cPhi>fPhiMax))
 324         continue;
 325
 326       Int_t lab = TMath::Abs(c->GetLabel());
 327       if (lab < fMinLabelClusters || lab > fMaxLabelClusters) {
 328         AliDebug(2,Form("Skipping cluster %d because label %d is not in range (%d, %d)", iClus, lab, fMinLabelClusters, fMaxLabelClusters));
 329         continue;
 330       }
 331
 332       // offset of 100 to skip ghost particles uid = -1
 333       AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
 334       Double_t e = TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz);
 335       fFastJetWrapper.AddInputVector(cPx, cPy, cPz, e, -iClus - 100);
 336     }
 337   }
 338
 339   // run jet finder
 340   fFastJetWrapper.Run();
 341 }
 342
 343 //________________________________________________________________________
 344 void AliEmcalJetTask::FillJetBranch()
 345 {
 346   // Fill jet branch with fastjet jets.
 347
 348   PrepareUtilities();
 349
 350   // loop over fastjet jets
 351   std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper.GetInclusiveJets();
 352   // sort jets according to jet pt
 353   static Int_t indexes[9999] = {-1};
 354   GetSortedArray(indexes, jets_incl);
 355
 356   AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
 357   for (UInt_t ijet = 0, jetCount = 0; ijet < jets_incl.size(); ++ijet) {
 358     Int_t ij = indexes[ijet];
 359     AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fFastJetWrapper.GetJetArea(ij)));
 360
 361     if (jets_incl[ij].perp() < fMinJetPt) continue;
 362     if (fFastJetWrapper.GetJetArea(ij) < fMinJetArea) continue;
 363     if ((jets_incl[ij].eta() < fJetEtaMin) || (jets_incl[ij].eta() > fJetEtaMax) ||
 364         (jets_incl[ij].phi() < fJetPhiMin) || (jets_incl[ij].phi() > fJetPhiMax))
 365       continue;
 366
 367     AliEmcalJet *jet = new ((*fJets)[jetCount])
 368       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
 369     jet->SetLabel(ij);
 370
 371     fastjet::PseudoJet area(fFastJetWrapper.GetJetAreaVector(ij));
 372     jet->SetArea(area.perp());
 373     jet->SetAreaEta(area.eta());
 374     jet->SetAreaPhi(area.phi());
 375     jet->SetAreaE(area.E());
 376
 377     // Fill constituent info
 378     std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper.GetJetConstituents(ij));
 379     FillJetConstituents(jet, constituents, fTracks, fClus, constituents);
 380
 381     if (fGeom) {
 382       if ((jet->Phi() > fGeom->GetArm1PhiMin() * TMath::DegToRad()) &&
 383           (jet->Phi() < fGeom->GetArm1PhiMax() * TMath::DegToRad()) &&
 384           (jet->Eta() > fGeom->GetArm1EtaMin()) &&
 385           (jet->Eta() < fGeom->GetArm1EtaMax()))
 386         jet->SetAxisInEmcal(kTRUE);
 387     }
 388
 389     ExecuteUtilities(jet, ij);
 390
 391     AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), jet->GetNumberOfConstituents()));
 392     jetCount++;
 393   }
 394
 395   TerminateUtilities();
 396 }
 397
 398 //________________________________________________________________________
 399 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
 400 {
 401   // Get the leading jets.
 402
 403   static Float_t pt[9999] = {0};
 404
 405   const Int_t n = (Int_t)array.size();
 406
 407   if (n < 1)
 408     return kFALSE;
 409
 410   for (Int_t i = 0; i < n; i++)
 411     pt[i] = array[i].perp();
 412
 413   TMath::Sort(n, pt, indexes);
 414
 415   return kTRUE;
 416 }
 417
 418 //________________________________________________________________________
 419 Bool_t AliEmcalJetTask::DoInit()
 420 {
 421   // Init. Return true if successful.
 422
 423   if (fTrackEfficiency < 1.) {
 424     if (gRandom) delete gRandom;
 425     gRandom = new TRandom3(0);
 426   }
 427
 428   // get geometry
 429   fGeom = AliEMCALGeometry::GetInstance();
 430   if (!fGeom) {
 431     AliWarning(Form("%s: Can not create EMCal geometry, some features will not work...", GetName()));
 432   }
 433
 434   // get input collections
 435   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 436
 437   // get the event
 438   fEvent = InputEvent();
 439   if (!fEvent) {
 440     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
 441     return 0;
 442   }
 443
 444   if (fTracksName == "Tracks")
 445     am->LoadBranch("Tracks");
 446   if (!fTracks && !fTracksName.IsNull()) {
 447     fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
 448     if (!fTracks) {
 449       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
 450       return 0;
 451     }
 452   }
 453
 454   if (fCaloName == "CaloClusters")
 455     am->LoadBranch("CaloClusters");
 456   if (!fClus && !fCaloName.IsNull()) {
 457     fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
 458     if (!fClus) {
 459       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
 460       return 0;
 461     }
 462   }
 463   if (fClus) {
 464     TClass cls(fClus->GetClass()->GetName());
 465     if (cls.InheritsFrom("AliEmcalParticle"))
 466       fIsEmcPart = 1;
 467   }
 468
 469   // add jets to event if not yet there
 470   if (!(fEvent->FindListObject(fJetsName))) {
 471     fJets = new TClonesArray("AliEmcalJet");
 472     fJets->SetName(fJetsName);
 473     fEvent->AddObject(fJets);
 474   }
 475   else {
 476     AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
 477     return 0;
 478   }
 479
 480   TString name("kt");
 481   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
 482   if ((fJetType & kAKT) != 0) {
 483     name  = "antikt";
 484     jalgo = fastjet::antikt_algorithm;
 485     AliDebug(1,"Using AKT algorithm");
 486   }
 487   else {
 488     AliDebug(1,"Using KT algorithm");
 489   }
 490
 491   if ((fJetType & kR020Jet) != 0)
 492     fRadius = 0.2;
 493   else if ((fJetType & kR030Jet) != 0)
 494     fRadius = 0.3;
 495   else if ((fJetType & kR040Jet) != 0)
 496     fRadius = 0.4;
 497
 498   // setup fj wrapper
 499   fFastJetWrapper.SetName(name);
 500   fFastJetWrapper.SetTitle(name);
 501   fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts);
 502   fFastJetWrapper.SetGhostArea(fGhostArea);
 503   fFastJetWrapper.SetR(fRadius);
 504   fFastJetWrapper.SetAlgorithm(jalgo);
 505   fFastJetWrapper.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
 506   fFastJetWrapper.SetMaxRap(fEtaMax);
 507
 508   // setting legacy mode
 509   if (fLegacyMode) {
 510     fFastJetWrapper.SetLegacyMode(kTRUE);
 511   }
 512
 513   InitUtilities();
 514
 515   return kTRUE;
 516 }
 517
 518 //___________________________________________________________________________________________________________________
 519 void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents, TClonesArray *tracks, TClonesArray *clusters,
 520                                           std::vector<fastjet::PseudoJet>& constituents_unsub, Int_t flag, TClonesArray *particles_sub)
 521 {
 522   Int_t nt            = 0;
 523   Int_t nc            = 0;
 524   Double_t neutralE   = 0.;
 525   Double_t maxCh      = 0.;
 526   Double_t maxNe      = 0.;
 527   Int_t gall          = 0;
 528   Int_t gemc          = 0;
 529   Int_t cemc          = 0;
 530   Int_t ncharged      = 0;
 531   Int_t nneutral      = 0;
 532   Double_t mcpt       = 0.;
 533   Double_t emcpt      = 0.;
 534
 535   Int_t uid   = -1;
 536
 537   jet->SetNumberOfTracks(constituents.size());
 538   jet->SetNumberOfClusters(constituents.size());
 539
 540   for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
 541
 542     if (flag == 0) uid = constituents[ic].user_index();
 543     else uid = GetIndexSub(constituents[ic].phi(), constituents_unsub);
 544
 545     AliDebug(3,Form("Processing constituent %d", uid));
 546     if (uid == -1) { //ghost particle
 547       ++gall;
 548       if (fGeom) {
 549         Double_t gphi = constituents[ic].phi();
 550         if (gphi < 0) gphi += TMath::TwoPi();
 551         gphi *= TMath::RadToDeg();
 552         Double_t geta = constituents[ic].eta();
 553         if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) &&
 554             (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax()))
 555           ++gemc;
 556       }
 557     }
 558     else if ((uid >= 100) && tracks) { // track constituent
 559       Int_t tid = uid - 100;
 560       AliVParticle *t = static_cast<AliVParticle*>(tracks->At(tid));
 561       if (!t) {
 562         AliError(Form("Could not find track %d",tid));
 563         continue;
 564       }
 565
 566       Double_t cEta = t->Eta();
 567       Double_t cPhi = t->Phi();
 568       Double_t cPt  = t->Pt();
 569       Double_t cP   = t->P();
 570       if (t->Charge() == 0) {
 571         neutralE += cP;
 572         ++nneutral;
 573         if (cPt > maxNe) maxNe = cPt;
 574       } else {
 575         ++ncharged;
 576         if (cPt > maxCh) maxCh = cPt;
 577       }
 578
 579       // check if MC particle
 580       if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt;
 581
 582       if (fGeom) {
 583         if (cPhi < 0) cPhi += TMath::TwoPi();
 584         cPhi *= TMath::RadToDeg();
 585         if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
 586             (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
 587           emcpt += cPt;
 588           ++cemc;
 589         }
 590       }
 591
 592       if (flag == 0 || particles_sub == 0) {
 593         jet->AddTrackAt(tid, nt);
 594       }
 595       else {
 596         Int_t part_sub_id = particles_sub->GetEntriesFast();
 597         AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(t);
 598         part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
 599         jet->AddTrackAt(part_sub_id, nt);
 600       }
 601
 602       ++nt;
 603     }
 604     else if ((uid <= -100) && clusters) { // cluster constituent
 605       Int_t cid = -uid - 100;
 606       AliVCluster *c = 0;
 607       Double_t cEta=0,cPhi=0,cPt=0,cP=0;
 608       if (fIsEmcPart) {
 609         AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
 610         if (!ep) continue;
 611         c = ep->GetCluster();
 612         if (!c) continue;
 613
 614         cEta = ep->Eta();
 615         cPhi = ep->Phi();
 616         cPt  = ep->Pt();
 617         cP   = ep->P();
 618       }
 619       else {
 620         c = static_cast<AliVCluster*>(fClus->At(cid));
 621         if (!c) continue;
 622
 623         TLorentzVector nP;
 624         c->GetMomentum(nP, fVertex);
 625         cEta = nP.Eta();
 626         cPhi = nP.Phi();
 627         cPt  = nP.Pt();
 628         cP   = nP.P();
 629       }
 630
 631       neutralE += cP;
 632       if (cPt > maxNe) maxNe = cPt;
 633
 634       // MC particle
 635       if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
 636
 637       if (fGeom) {
 638         if (cPhi<0) cPhi += TMath::TwoPi();
 639         cPhi *= TMath::RadToDeg();
 640         if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
 641             (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
 642           emcpt += cPt;
 643           ++cemc;
 644         }
 645       }
 646
 647       if (flag == 0 || particles_sub == 0) {
 648         jet->AddClusterAt(cid, nc);
 649       }
 650       else {
 651         Int_t part_sub_id = particles_sub->GetEntriesFast();
 652         AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c);
 653         part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
 654         jet->AddTrackAt(part_sub_id, nt);
 655       }
 656
 657       ++nc;
 658       ++nneutral;
 659     }
 660     else {
 661       AliError(Form("%s: No logical way to end up here.", GetName()));
 662       continue;
 663     }
 664   }
 665
 666   jet->SetNumberOfTracks(nt);
 667   jet->SetNumberOfClusters(nc);
 668   jet->SetNEF(neutralE / jet->E());
 669   jet->SetMaxChargedPt(maxCh);
 670   jet->SetMaxNeutralPt(maxNe);
 671   if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall);
 672   else jet->SetAreaEmc(-1);
 673   jet->SetNEmc(cemc);
 674   jet->SetNumberOfCharged(ncharged);
 675   jet->SetNumberOfNeutrals(nneutral);
 676   jet->SetMCPt(mcpt);
 677   jet->SetPtEmc(emcpt);
 678   jet->SortConstituents();
 679 }
 680
 681 //______________________________________________________________________________________
 682 Int_t AliEmcalJetTask::GetIndexSub(Double_t phi_sub, std::vector<fastjet::PseudoJet>& constituents_unsub)
 683 {
 684   for (UInt_t ii = 0; ii < constituents_unsub.size(); ii++) {
 685     Double_t phi_unsub = constituents_unsub[ii].phi();
 686     if (phi_sub == phi_unsub) return constituents_unsub[ii].user_index();
 687   }
 688
 689   return -1;
 690 }