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