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 "AliRhoParameter.h"
  28 using std::cout;
  29 using std::endl;
  30 using std::cerr;
  31
  32 ClassImp(AliEmcalJetTask)
  33
  34 //________________________________________________________________________
  35 AliEmcalJetTask::AliEmcalJetTask() :
  36   AliAnalysisTaskSE("AliEmcalJetTask"),
  37   fTracksName("Tracks"),
  38   fCaloName("CaloClusters"),
  39   fJetsName("Jets"),
  40   fJetsSubName(""),
  41   fJetType(kNone),
  42   fConstSel(0),
  43   fMCConstSel(0),
  44   fMarkConst(kFALSE),
  45   fRadius(0.4),
  46   fMinJetTrackPt(0.15),
  47   fMinJetClusPt(0.15),
  48   fPhiMin(0),
  49   fPhiMax(TMath::TwoPi()),
  50   fEtaMin(-0.9),
  51   fEtaMax(+0.9),
  52   fMinJetArea(0.005),
  53   fMinJetPt(1.0),
  54   fJetPhiMin(-10),
  55   fJetPhiMax(+10),
  56   fJetEtaMin(-1),
  57   fJetEtaMax(+1),
  58   fGhostArea(0.005),
  59   fMinMCLabel(0),
  60   fRecombScheme(fastjet::pt_scheme),
  61   fTrackEfficiency(1.),
  62   fMCFlag(0),
  63   fGeneratorIndex(-1),
  64   fIsInit(0),
  65   fLocked(0),
  66   fIsPSelSet(0),
  67   fIsMcPart(0),
  68   fIsEmcPart(0),
  69   fLegacyMode(kFALSE),
  70   fCodeDebug(kFALSE),
  71   fDoGenericSubtractionJetMass(kFALSE),
  72   fDoGenericSubtractionGR(kFALSE),
  73   fDoConstituentSubtraction(kFALSE),
  74   fUseExternalBkg(kFALSE),
  75   fRhoName(""),
  76   fRhomName(""),
  77   fRho(0),
  78   fRhom(0),
  79   fRMax(0.4),
  80   fDRStep(0.04),
  81   fPtMinGR(40.),
  82   fJets(0),
  83   fJetsSub(0),
  84   fEvent(0),
  85   fTracks(0),
  86   fClus(0),
  87   fRhoParam(0),
  88   fRhomParam(0)
  89 {
  90   // Default constructor.
  91 }
  92
  93 //________________________________________________________________________
  94 AliEmcalJetTask::AliEmcalJetTask(const char *name) :
  95   AliAnalysisTaskSE(name),
  96   fTracksName("Tracks"),
  97   fCaloName("CaloClusters"),
  98   fJetsName("Jets"),
  99   fJetsSubName(""),
 100   fJetType(kAKT|kFullJet|kRX1Jet),
 101   fConstSel(0),
 102   fMCConstSel(0),
 103   fMarkConst(kFALSE),
 104   fRadius(0.4),
 105   fMinJetTrackPt(0.15),
 106   fMinJetClusPt(0.15),
 107   fPhiMin(0),
 108   fPhiMax(TMath::TwoPi()),
 109   fEtaMin(-0.9),
 110   fEtaMax(+0.9),
 111   fMinJetArea(0.001),
 112   fMinJetPt(1.0),
 113   fJetPhiMin(-10),
 114   fJetPhiMax(+10),
 115   fJetEtaMin(-1),
 116   fJetEtaMax(+1),
 117   fGhostArea(0.005),
 118   fMinMCLabel(0),
 119   fRecombScheme(fastjet::pt_scheme),
 120   fTrackEfficiency(1.),
 121   fMCFlag(0),
 122   fGeneratorIndex(-1),
 123   fIsInit(0),
 124   fLocked(0),
 125   fIsPSelSet(0),
 126   fIsMcPart(0),
 127   fIsEmcPart(0),
 128   fLegacyMode(kFALSE),
 129   fCodeDebug(kFALSE),
 130   fDoGenericSubtractionJetMass(kFALSE),
 131   fDoGenericSubtractionGR(kFALSE),
 132   fDoConstituentSubtraction(kFALSE),
 133   fUseExternalBkg(kFALSE),
 134   fRhoName(""),
 135   fRhomName(""),
 136   fRho(0),
 137   fRhom(0),
 138   fRMax(0.4),
 139   fDRStep(0.04),
 140   fPtMinGR(40.),
 141   fJets(0),
 142   fJetsSub(0),
 143   fEvent(0),
 144   fTracks(0),
 145   fClus(0),
 146   fRhoParam(0),
 147   fRhomParam(0)
 148 {
 149   // Standard constructor.
 150
 151   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
 152 }
 153
 154 //________________________________________________________________________
 155 AliEmcalJetTask::~AliEmcalJetTask()
 156 {
 157   // Destructor
 158 }
 159
 160 //________________________________________________________________________
 161 void AliEmcalJetTask::UserCreateOutputObjects()
 162 {
 163   // Create user objects.
 164
 165   fJets = new TClonesArray("AliEmcalJet");
 166   fJets->SetName(fJetsName);
 167
 168   if(!fJetsSubName.IsNull()) {
 169     fJetsSub = new TClonesArray("AliEmcalJet");
 170     fJetsSub->SetName(fJetsSubName);
 171   }
 172 }
 173
 174 //________________________________________________________________________
 175 void AliEmcalJetTask::UserExec(Option_t *)
 176 {
 177   // Main loop, called for each event.
 178   if (!fIsInit) {
 179     if (!DoInit())
 180       return;
 181     fIsInit = kTRUE;
 182   }
 183
 184   // clear the jet array (normally a null operation)
 185   fJets->Delete();
 186   if(fJetsSub) fJetsSub->Delete();
 187
 188   FindJets();
 189 }
 190
 191 //________________________________________________________________________
 192 void AliEmcalJetTask::Terminate(Option_t *)
 193 {
 194   // Called once at the end of the analysis.
 195 }
 196
 197 //________________________________________________________________________
 198 void AliEmcalJetTask::FindJets()
 199 {
 200   // Find jets.
 201   if (!fTracks && !fClus){
 202     cout << "WARNING NO TRACKS OR CLUSTERS:"  <<endl;
 203     return;
 204   }
 205
 206  if (fRhoParam)
 207     fRho = fRhoParam->GetVal();
 208   if (fRhomParam)
 209     fRhom = fRhomParam->GetVal();
 210
 211   TString name("kt");
 212   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
 213   if ((fJetType & kAKT) != 0) {
 214     name  = "antikt";
 215     jalgo = fastjet::antikt_algorithm;
 216     AliDebug(1,"Using AKT algorithm");
 217   }
 218   else {
 219     AliDebug(1,"Using KT algorithm");
 220   }
 221
 222   if ((fJetType & kR020Jet) != 0)
 223     fRadius = 0.2;
 224   else if ((fJetType & kR030Jet) != 0)
 225     fRadius = 0.3;
 226   else if ((fJetType & kR040Jet) != 0)
 227     fRadius = 0.4;
 228
 229   // setup fj wrapper
 230   AliFJWrapper fjw(name, name);
 231   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
 232   fjw.SetGhostArea(fGhostArea);
 233   fjw.SetR(fRadius);
 234   fjw.SetAlgorithm(jalgo);
 235   fjw.SetRecombScheme(static_cast<fastjet::RecombinationScheme>(fRecombScheme));
 236   fjw.SetMaxRap(fEtaMax);
 237   fjw.Clear();
 238
 239   // get primary vertex
 240   Double_t vertex[3] = {0, 0, 0};
 241   if(fEvent->GetPrimaryVertex()) fEvent->GetPrimaryVertex()->GetXYZ(vertex);
 242
 243   AliDebug(2,Form("Jet type = %d", fJetType));
 244
 245   if ((fIsMcPart || ((fJetType & kFullJet) != 0) || ((fJetType & kChargedJet) != 0)) && fTracks) {
 246     const Int_t Ntracks = fTracks->GetEntries();
 247     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 248       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
 249       if (!t)
 250         continue;
 251       if (fIsMcPart) {
 252         if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
 253           AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
 254           continue;
 255         }
 256         if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
 257           AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
 258           continue;
 259         }
 260       }
 261       if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) {
 262         if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
 263           AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
 264           continue;
 265         }
 266         else {
 267           AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
 268         }
 269       }
 270       else {
 271         if (t->TestBits(fConstSel) != (Int_t)fConstSel) {
 272           AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
 273           continue;
 274         }
 275         else {
 276           AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
 277         }
 278       }
 279       if ((t->GetFlag() & fMCFlag) != fMCFlag) {
 280         AliDebug(2,Form("Skipping track %d because it does not match the MC flags", iTracks));
 281         continue;
 282       }
 283       if (fGeneratorIndex >= 0 && t->GetGeneratorIndex() != fGeneratorIndex) {
 284         AliDebug(2,Form("Skipping track %d because it does not match the MC generator index", iTracks));
 285         continue;
 286       }
 287       if (t->Pt() < fMinJetTrackPt)
 288         continue;
 289       Double_t eta = t->Eta();
 290       Double_t phi = t->Phi();
 291       if ((eta<fEtaMin) || (eta>fEtaMax) ||
 292           (phi<fPhiMin) || (phi>fPhiMax))
 293         continue;
 294
 295       // artificial inefficiency
 296       if (fTrackEfficiency < 1.) {
 297         Double_t rnd = gRandom->Rndm();
 298         if (fTrackEfficiency < rnd) {
 299           AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", iTracks));
 300           continue;
 301         }
 302       }
 303
 304       // offset of 100 for consistency with cluster ids
 305       AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt()));
 306       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->E(), iTracks + 100);
 307     }
 308   }
 309
 310   if ((((fJetType & kFullJet) != 0) || ((fJetType & kNeutralJet) != 0)) && fClus) {
 311     const Int_t Nclus = fClus->GetEntries();
 312     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 313       AliVCluster *c = 0;
 314       Double_t cEta=0,cPhi=0,cPt=0;
 315       Double_t cPx=0,cPy=0,cPz=0;
 316       if (fIsEmcPart) {
 317         AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
 318         if (!ep)
 319           continue;
 320
 321         c = ep->GetCluster();
 322         if (!c)
 323           continue;
 324
 325         if (c->GetLabel() > fMinMCLabel) {
 326           if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
 327             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
 328             continue;
 329           }
 330           else {
 331             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
 332           }
 333         }
 334         else {
 335           if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
 336             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
 337             continue;
 338           }
 339           else {
 340             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
 341           }
 342         }
 343
 344         cEta = ep->Eta();
 345         cPhi = ep->Phi();
 346         cPt  = ep->Pt();
 347         cPx  = ep->Px();
 348         cPy  = ep->Py();
 349         cPz  = ep->Pz();
 350       } else {
 351         c = static_cast<AliVCluster*>(fClus->At(iClus));
 352         if (!c)
 353           continue;
 354
 355         if (c->GetLabel() > fMinMCLabel) {
 356           if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
 357             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
 358             continue;
 359           }
 360           else {
 361             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
 362           }
 363         }
 364         else {
 365           if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
 366             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
 367             continue;
 368           }
 369           else {
 370             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
 371           }
 372         }
 373
 374         TLorentzVector nP;
 375         c->GetMomentum(nP, vertex);
 376         cEta = nP.Eta();
 377         cPhi = nP.Phi();
 378         cPt  = nP.Pt();
 379         cPx  = nP.Px();
 380         cPy  = nP.Py();
 381         cPz  = nP.Pz();
 382       }
 383       if (!c->IsEMCAL())
 384         continue;
 385       if (cPt < fMinJetClusPt)
 386         continue;
 387       if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
 388           (cPhi<fPhiMin) || (cPhi>fPhiMax))
 389         continue;
 390       // offset of 100 to skip ghost particles uid = -1
 391       AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
 392       fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
 393     }
 394   }
 395
 396   // setting legacy mode
 397   if (fLegacyMode) {
 398     fjw.SetLegacyMode(kTRUE);
 399   }
 400
 401   // run jet finder
 402   fjw.Run();
 403
 404   //run generic subtractor
 405   if(fDoGenericSubtractionJetMass) {
 406     fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
 407     fjw.DoGenericSubtractionJetMass();
 408   }
 409
 410   //run constituent subtractor
 411   if(fDoConstituentSubtraction) {
 412     fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
 413     fjw.DoConstituentSubtraction();
 414   }
 415
 416   // get geometry
 417   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 418   if (!geom) {
 419     AliFatal(Form("%s: Can not create geometry", GetName()));
 420     return;
 421   }
 422
 423   // loop over fastjet jets
 424   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
 425   // sort jets according to jet pt
 426   static Int_t indexes[9999] = {-1};
 427   GetSortedArray(indexes, jets_incl);
 428
 429   AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
 430   for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
 431     Int_t ij = indexes[ijet];
 432     AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
 433
 434     if (jets_incl[ij].perp()<fMinJetPt)
 435       continue;
 436     if (fjw.GetJetArea(ij)<fMinJetArea)
 437       continue;
 438     if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
 439         (jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
 440       continue;
 441
 442     AliEmcalJet *jet = new ((*fJets)[jetCount])
 443       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
 444     jet->SetLabel(ij);
 445
 446     //do generic subtraction if requested
 447 #ifdef FASTJET_VERSION
 448     if(fDoGenericSubtractionJetMass) {
 449       std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
 450       Int_t n = (Int_t)jetMassInfo.size();
 451       if(n>ij && n>0) {
 452         jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
 453         jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
 454         jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
 455         jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
 456       }
 457     }
 458     //here do generic subtraction for angular structure function
 459     Double_t ptcorr = jets_incl[ij].perp()-fjw.GetJetArea(ij)*fRho;
 460     fRMax = fRadius+0.2;
 461     if(fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
 462       fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
 463       fjw.SetRMaxAndStep(fRMax,fDRStep);
 464       fjw.DoGenericSubtractionGR(ij);
 465       std::vector<double> num = fjw.GetGRNumerator();
 466       std::vector<double> den = fjw.GetGRDenominator();
 467       std::vector<double> nums = fjw.GetGRNumeratorSub();
 468       std::vector<double> dens = fjw.GetGRDenominatorSub();
 469       //pass this to AliEmcalJet
 470       jet->SetGRNumSize(num.size());
 471       jet->SetGRDenSize(den.size());
 472       jet->SetGRNumSubSize(nums.size());
 473       jet->SetGRDenSubSize(dens.size());
 474       Int_t nsize = (Int_t)num.size();
 475       for(Int_t g = 0; g<nsize; ++g) {
 476         jet->AddGRNumAt(num[g],g);
 477         jet->AddGRNumSubAt(nums[g],g);
 478       }
 479       Int_t dsize = (Int_t)den.size();
 480       for(Int_t g = 0; g<dsize; ++g) {
 481         jet->AddGRDenAt(den[g],g);
 482         jet->AddGRDenSubAt(dens[g],g);
 483       }
 484     }
 485 #endif
 486
 487     // loop over constituents
 488     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
 489     jet->SetNumberOfTracks(constituents.size());
 490     jet->SetNumberOfClusters(constituents.size());
 491
 492     Int_t nt            = 0;
 493     Int_t nc            = 0;
 494     Double_t neutralE   = 0;
 495     Double_t maxCh      = 0;
 496     Double_t maxNe      = 0;
 497     Int_t gall          = 0;
 498     Int_t gemc          = 0;
 499     Int_t cemc          = 0;
 500     Int_t ncharged      = 0;
 501     Int_t nneutral      = 0;
 502     Double_t mcpt       = 0;
 503     Double_t emcpt      = 0;
 504
 505     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
 506       Int_t uid = constituents[ic].user_index();
 507       AliDebug(3,Form("Processing constituent %d", uid));
 508       if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
 509         ++gall;
 510         Double_t gphi = constituents[ic].phi();
 511         if (gphi<0)
 512           gphi += TMath::TwoPi();
 513         gphi *= TMath::RadToDeg();
 514         Double_t geta = constituents[ic].eta();
 515         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
 516             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
 517           ++gemc;
 518       } else if ((uid > 0) && fTracks) { // track constituent
 519         Int_t tid = uid - 100;
 520         AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
 521         if (!t) {
 522           AliError(Form("Could not find track %d",tid));
 523           continue;
 524         }
 525         if (jetCount < fMarkConst) {
 526           AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
 527           t->SetBit(fJetType);
 528         }
 529         Double_t cEta = t->Eta();
 530         Double_t cPhi = t->Phi();
 531         Double_t cPt  = t->Pt();
 532         Double_t cP   = t->P();
 533         if (t->Charge() == 0) {
 534           neutralE += cP;
 535           ++nneutral;
 536           if (cPt > maxNe)
 537             maxNe = cPt;
 538         } else {
 539           ++ncharged;
 540           if (cPt > maxCh)
 541             maxCh = cPt;
 542         }
 543         if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
 544           mcpt += cPt;
 545
 546         if (cPhi<0)
 547           cPhi += TMath::TwoPi();
 548         cPhi *= TMath::RadToDeg();
 549         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
 550             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
 551           emcpt += cPt;
 552           ++cemc;
 553         }
 554
 555         jet->AddTrackAt(tid, nt);
 556         ++nt;
 557       } else if (fClus) { // cluster constituent
 558         Int_t cid = -uid - 100;
 559         AliVCluster *c = 0;
 560         Double_t cEta=0,cPhi=0,cPt=0,cP=0;
 561         if (fIsEmcPart) {
 562           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
 563           if (!ep)
 564             continue;
 565           c = ep->GetCluster();
 566           if (!c)
 567             continue;
 568           if (jetCount < fMarkConst)
 569             ep->SetBit(fJetType);
 570           cEta = ep->Eta();
 571           cPhi = ep->Phi();
 572           cPt  = ep->Pt();
 573           cP   = ep->P();
 574         } else {
 575           c = static_cast<AliVCluster*>(fClus->At(cid));
 576           if (!c)
 577             continue;
 578           if (jetCount < fMarkConst)
 579             c->SetBit(fJetType);
 580           TLorentzVector nP;
 581           c->GetMomentum(nP, vertex);
 582           cEta = nP.Eta();
 583           cPhi = nP.Phi();
 584           cPt  = nP.Pt();
 585           cP   = nP.P();
 586         }
 587
 588         neutralE += cP;
 589         if (cPt > maxNe)
 590           maxNe = cPt;
 591
 592         if (c->GetLabel() > fMinMCLabel) // MC particle
 593           mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
 594
 595         if (cPhi<0)
 596           cPhi += TMath::TwoPi();
 597         cPhi *= TMath::RadToDeg();
 598         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
 599             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
 600           emcpt += cPt;
 601           ++cemc;
 602         }
 603
 604         jet->AddClusterAt(cid, nc);
 605         ++nc;
 606         ++nneutral;
 607       } else {
 608         AliError(Form("%s: No logical way to end up here.", GetName()));
 609         continue;
 610       }
 611     } /* end of constituent loop */
 612
 613     jet->SetNumberOfTracks(nt);
 614     jet->SetNumberOfClusters(nc);
 615     jet->SortConstituents();
 616     jet->SetMaxNeutralPt(maxNe);
 617     jet->SetMaxChargedPt(maxCh);
 618     jet->SetNEF(neutralE / jet->E());
 619     fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
 620     jet->SetArea(area.perp());
 621     jet->SetAreaEta(area.eta());
 622     jet->SetAreaPhi(area.phi());
 623     jet->SetNumberOfCharged(ncharged);
 624     jet->SetNumberOfNeutrals(nneutral);
 625     jet->SetMCPt(mcpt);
 626     jet->SetNEmc(cemc);
 627     jet->SetPtEmc(emcpt);
 628
 629     if (gall > 0)
 630       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
 631     else
 632       jet->SetAreaEmc(-1);
 633     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
 634         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
 635         (jet->Eta() > geom->GetArm1EtaMin()) &&
 636         (jet->Eta() < geom->GetArm1EtaMax()))
 637       jet->SetAxisInEmcal(kTRUE);
 638
 639     AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
 640     jetCount++;
 641   }
 642   //fJets->Sort();
 643
 644   //run constituent subtractor if requested
 645   if(fDoConstituentSubtraction) {
 646 #ifdef FASTJET_VERSION
 647     if(!fJetsSub) AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s",fJetsSubName.Data()));
 648     else {
 649       std::vector<fastjet::PseudoJet> jets_sub;
 650       jets_sub = fjw.GetConstituentSubtrJets();
 651       AliDebug(1,Form("%d constituent subtracted jets found", (Int_t)jets_sub.size()));
 652       for (UInt_t ijet=0, jetCount=0; ijet<jets_sub.size(); ++ijet) {
 653         //Only storing 4-vector and jet area of unsubtracted jet
 654         AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount])
 655           AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
 656         jet_sub->SetLabel(ijet);
 657         fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
 658         jet_sub->SetArea(area.perp());
 659         jet_sub->SetAreaEta(area.eta());
 660         jet_sub->SetAreaPhi(area.phi());
 661         jet_sub->SetAreaEmc(area.perp());
 662         jetCount++;
 663       }
 664     }
 665 #endif
 666   } //constituent subtraction
 667 }
 668
 669 //________________________________________________________________________
 670 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
 671 {
 672   // Get the leading jets.
 673
 674   static Float_t pt[9999] = {0};
 675
 676   const Int_t n = (Int_t)array.size();
 677
 678   if (n < 1)
 679     return kFALSE;
 680
 681   for (Int_t i = 0; i < n; i++)
 682     pt[i] = array[i].perp();
 683
 684   TMath::Sort(n, pt, indexes);
 685
 686   return kTRUE;
 687 }
 688
 689 //________________________________________________________________________
 690 Bool_t AliEmcalJetTask::DoInit()
 691 {
 692   // Init. Return true if successful.
 693
 694   if (fTrackEfficiency < 1.) {
 695     if (gRandom) delete gRandom;
 696     gRandom = new TRandom3(0);
 697   }
 698
 699   // get input collections
 700   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 701
 702   // get the event
 703   fEvent = InputEvent();
 704   if (!fEvent) {
 705     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
 706     return 0;
 707   }
 708
 709   // add jets to event if not yet there
 710   if (!(fEvent->FindListObject(fJetsName)))
 711     fEvent->AddObject(fJets);
 712   else {
 713     AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
 714     return 0;
 715   }
 716
 717   if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub)
 718     fEvent->AddObject(fJetsSub);
 719
 720   if (fTracksName == "Tracks")
 721     am->LoadBranch("Tracks");
 722   if (!fTracks && !fTracksName.IsNull()) {
 723     fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
 724     if (!fTracks) {
 725       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
 726       return 0;
 727     }
 728   }
 729   if (fTracks) {
 730     TClass cls(fTracks->GetClass()->GetName());
 731     if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
 732       fIsMcPart = 1;
 733   }
 734
 735   if (fCaloName == "CaloClusters")
 736     am->LoadBranch("CaloClusters");
 737   if (!fClus && !fCaloName.IsNull()) {
 738     fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
 739     if (!fClus) {
 740       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
 741       return 0;
 742     }
 743   }
 744   if (fClus) {
 745     TClass cls(fClus->GetClass()->GetName());
 746     if (cls.InheritsFrom("AliEmcalParticle"))
 747       fIsEmcPart = 1;
 748   }
 749
 750   if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
 751     fRhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
 752     if (!fRhoParam) {
 753       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
 754       return 0;
 755     }
 756   }
 757   if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
 758     fRhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhomName));
 759     if (!fRhomParam) {
 760       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhomName.Data()));
 761       return 0;
 762     }
 763   }
 764
 765   return 1;
 766 }