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