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.DoGenericSubtractionJetSigma2();
 423     fjw.DoGenericSubtractionJetConstituent();
 424     fjw.DoGenericSubtractionJetLeSub();
 425   }
 426
 427   //run constituent subtractor
 428   if(fDoConstituentSubtraction) {
 429     fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
 430     fjw.DoConstituentSubtraction();
 431   }
 432
 433   // get geometry
 434   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 435   if (!geom) {
 436     AliFatal(Form("%s: Can not create geometry", GetName()));
 437     return;
 438   }
 439
 440   // loop over fastjet jets
 441   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
 442   // sort jets according to jet pt
 443   static Int_t indexes[9999] = {-1};
 444   GetSortedArray(indexes, jets_incl);
 445
 446   AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
 447   for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
 448     Int_t ij = indexes[ijet];
 449     AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
 450
 451     if (jets_incl[ij].perp()<fMinJetPt)
 452       continue;
 453     if (fjw.GetJetArea(ij)<fMinJetArea)
 454       continue;
 455     if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
 456         (jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
 457       continue;
 458
 459     AliEmcalJet *jet = new ((*fJets)[jetCount])
 460       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
 461     jet->SetLabel(ij);
 462
 463     //do generic subtraction if requested
 464 #ifdef FASTJET_VERSION
 465     if(fDoGenericSubtractionJetMass) {
 466       std::vector<fastjet::contrib::GenericSubtractorInfo> jetMassInfo = fjw.GetGenSubtractorInfoJetMass();
 467       Int_t n = (Int_t)jetMassInfo.size();
 468       if(n>ij && n>0) {
 469         jet->SetFirstDerivative(jetMassInfo[ij].first_derivative());
 470         jet->SetSecondDerivative(jetMassInfo[ij].second_derivative());
 471         jet->SetFirstOrderSubtracted(jetMassInfo[ij].first_order_subtracted());
 472         jet->SetSecondOrderSubtracted(jetMassInfo[ij].second_order_subtracted());
 473       }
 474     }
 475
 476     //here do generic subtraction for angular structure function
 477     Double_t ptcorr = jets_incl[ij].perp()-fjw.GetJetArea(ij)*fRho;
 478     fRMax = fRadius+0.2;
 479     if(fDoGenericSubtractionGR && ptcorr>fPtMinGR) {
 480       fjw.SetUseExternalBkg(fUseExternalBkg,fRho,fRhom);
 481       fjw.SetRMaxAndStep(fRMax,fDRStep);
 482       fjw.DoGenericSubtractionGR(ij);
 483       std::vector<double> num = fjw.GetGRNumerator();
 484       std::vector<double> den = fjw.GetGRDenominator();
 485       std::vector<double> nums = fjw.GetGRNumeratorSub();
 486       std::vector<double> dens = fjw.GetGRDenominatorSub();
 487       //pass this to AliEmcalJet
 488       jet->SetGRNumSize(num.size());
 489       jet->SetGRDenSize(den.size());
 490       jet->SetGRNumSubSize(nums.size());
 491       jet->SetGRDenSubSize(dens.size());
 492       Int_t nsize = (Int_t)num.size();
 493       for(Int_t g = 0; g<nsize; ++g) {
 494         jet->AddGRNumAt(num[g],g);
 495         jet->AddGRNumSubAt(nums[g],g);
 496       }
 497       Int_t dsize = (Int_t)den.size();
 498       for(Int_t g = 0; g<dsize; ++g) {
 499         jet->AddGRDenAt(den[g],g);
 500         jet->AddGRDenSubAt(dens[g],g);
 501       }
 502     }
 503
 504    if(fDoGenericSubtractionExtraJetShapes){
 505       std::vector<fastjet::contrib::GenericSubtractorInfo> jetAngularityInfo = fjw.GetGenSubtractorInfoJetAngularity();
 506       Int_t na = (Int_t)jetAngularityInfo.size();
 507       if(na>ij && na>0) {
 508         jet->SetFirstDerivativeAngularity(jetAngularityInfo[ij].first_derivative());
 509         jet->SetSecondDerivativeAngularity(jetAngularityInfo[ij].second_derivative());
 510         jet->SetFirstOrderSubtractedAngularity(jetAngularityInfo[ij].first_order_subtracted());
 511         jet->SetSecondOrderSubtractedAngularity(jetAngularityInfo[ij].second_order_subtracted());
 512       }
 513
 514       std::vector<fastjet::contrib::GenericSubtractorInfo> jetpTDInfo = fjw.GetGenSubtractorInfoJetpTD();
 515       Int_t np = (Int_t)jetpTDInfo.size();
 516       if(np>ij && np>0) {
 517         jet->SetFirstDerivativepTD(jetpTDInfo[ij].first_derivative());
 518         jet->SetSecondDerivativepTD(jetpTDInfo[ij].second_derivative());
 519         jet->SetFirstOrderSubtractedpTD(jetpTDInfo[ij].first_order_subtracted());
 520         jet->SetSecondOrderSubtractedpTD(jetpTDInfo[ij].second_order_subtracted());
 521       }
 522
 523       std::vector<fastjet::contrib::GenericSubtractorInfo> jetCircularityInfo = fjw.GetGenSubtractorInfoJetCircularity();
 524       Int_t nc = (Int_t)jetCircularityInfo.size();
 525       if(nc>ij && nc>0) {
 526         jet->SetFirstDerivativeCircularity(jetCircularityInfo[ij].first_derivative());
 527         jet->SetSecondDerivativeCircularity(jetCircularityInfo[ij].second_derivative());
 528         jet->SetFirstOrderSubtractedCircularity(jetCircularityInfo[ij].first_order_subtracted());
 529         jet->SetSecondOrderSubtractedCircularity(jetCircularityInfo[ij].second_order_subtracted());
 530       }
 531
 532       std::vector<fastjet::contrib::GenericSubtractorInfo> jetSigma2Info = fjw.GetGenSubtractorInfoJetSigma2();
 533       Int_t ns = (Int_t)jetSigma2Info.size();
 534       if(ns>ij && ns>0) {
 535         jet->SetFirstDerivativeSigma2(jetSigma2Info[ij].first_derivative());
 536         jet->SetSecondDerivativeSigma2(jetSigma2Info[ij].second_derivative());
 537         jet->SetFirstOrderSubtractedSigma2(jetSigma2Info[ij].first_order_subtracted());
 538         jet->SetSecondOrderSubtractedSigma2(jetSigma2Info[ij].second_order_subtracted());
 539       }
 540
 541
 542       std::vector<fastjet::contrib::GenericSubtractorInfo> jetConstituentInfo = fjw.GetGenSubtractorInfoJetConstituent();
 543       Int_t nco= (Int_t)jetConstituentInfo.size();
 544       if(nco>ij && nco>0) {
 545         jet->SetFirstDerivativeConstituent(jetConstituentInfo[ij].first_derivative());
 546         jet->SetSecondDerivativeConstituent(jetConstituentInfo[ij].second_derivative());
 547         jet->SetFirstOrderSubtractedConstituent(jetConstituentInfo[ij].first_order_subtracted());
 548         jet->SetSecondOrderSubtractedConstituent(jetConstituentInfo[ij].second_order_subtracted());
 549       }
 550
 551       std::vector<fastjet::contrib::GenericSubtractorInfo> jetLeSubInfo = fjw.GetGenSubtractorInfoJetLeSub();
 552       Int_t nlsub= (Int_t)jetLeSubInfo.size();
 553       if(nlsub>ij && nlsub>0) {
 554         jet->SetFirstDerivativeLeSub(jetLeSubInfo[ij].first_derivative());
 555         jet->SetSecondDerivativeLeSub(jetLeSubInfo[ij].second_derivative());
 556         jet->SetFirstOrderSubtractedLeSub(jetLeSubInfo[ij].first_order_subtracted());
 557         jet->SetSecondOrderSubtractedLeSub(jetLeSubInfo[ij].second_order_subtracted());
 558       }
 559    }
 560 #endif
 561
 562     // Fill constituent info
 563     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
 564     jet->SetNumberOfTracks(constituents.size());
 565     jet->SetNumberOfClusters(constituents.size());
 566
 567     Int_t nt            = 0;
 568     Int_t nc            = 0;
 569     Double_t neutralE   = 0;
 570     Double_t maxCh      = 0;
 571     Double_t maxNe      = 0;
 572     Int_t gall          = 0;
 573     Int_t gemc          = 0;
 574     Int_t cemc          = 0;
 575     Int_t ncharged      = 0;
 576     Int_t nneutral      = 0;
 577     Double_t mcpt       = 0;
 578     Double_t emcpt      = 0;
 579
 580     FillJetConstituents(constituents,jet,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc,constituents,0);
 581     jet->SetNumberOfTracks(nt);
 582     jet->SetNumberOfClusters(nc);
 583     jet->SortConstituents();
 584     jet->SetMaxNeutralPt(maxNe);
 585     jet->SetMaxChargedPt(maxCh);
 586     jet->SetNEF(neutralE / jet->E());
 587     fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
 588     jet->SetArea(area.perp());
 589     jet->SetAreaEta(area.eta());
 590     jet->SetAreaPhi(area.phi());
 591     jet->SetNumberOfCharged(ncharged);
 592     jet->SetNumberOfNeutrals(nneutral);
 593     jet->SetMCPt(mcpt);
 594     jet->SetNEmc(cemc);
 595     jet->SetPtEmc(emcpt);
 596
 597     if (gall > 0)
 598       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
 599     else
 600       jet->SetAreaEmc(-1);
 601     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
 602         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
 603         (jet->Eta() > geom->GetArm1EtaMin()) &&
 604         (jet->Eta() < geom->GetArm1EtaMax()))
 605       jet->SetAxisInEmcal(kTRUE);
 606
 607     AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
 608     jetCount++;
 609   }
 610   //fJets->Sort();
 611
 612   //run constituent subtractor if requested
 613   if(fDoConstituentSubtraction) {
 614 #ifdef FASTJET_VERSION
 615     if(!fJetsSub) AliWarning(Form("No jet branch to write to for subtracted jets. fJetsSubName: %s",fJetsSubName.Data()));
 616     else {
 617       std::vector<fastjet::PseudoJet> jets_sub;
 618       jets_sub = fjw.GetConstituentSubtrJets();
 619       AliDebug(1,Form("%d constituent subtracted jets found", (Int_t)jets_sub.size()));
 620       for (UInt_t ijet=0, jetCount=0; ijet<jets_sub.size(); ++ijet) {
 621         //Only storing 4-vector and jet area of unsubtracted jet
 622         AliEmcalJet *jet_sub = new ((*fJetsSub)[jetCount])
 623           AliEmcalJet(jets_sub[ijet].perp(), jets_sub[ijet].eta(), jets_sub[ijet].phi(), jets_sub[ijet].m());
 624         jet_sub->SetLabel(ijet);
 625          // Fill constituent info
 626         std::vector<fastjet::PseudoJet> constituents_unsub(fjw.GetJetConstituents(ijet));
 627         std::vector<fastjet::PseudoJet> constituents_sub = jets_sub[ijet].constituents();
 628         jet_sub->SetNumberOfTracks(constituents_sub.size());
 629         jet_sub->SetNumberOfClusters(constituents_sub.size());
 630         Int_t nt            = 0;
 631         Int_t nc            = 0;
 632         Double_t neutralE   = 0;
 633         Double_t maxCh      = 0;
 634         Double_t maxNe      = 0;
 635         Int_t gall          = 0;
 636         Int_t gemc          = 0;
 637         Int_t cemc          = 0;
 638         Int_t ncharged      = 0;
 639         Int_t nneutral      = 0;
 640         Double_t mcpt       = 0;
 641         Double_t emcpt      = 0;
 642
 643         FillJetConstituents(constituents_sub,jet_sub,vertex,jetCount,nt,nc,maxCh,maxNe,ncharged,nneutral,neutralE,mcpt,cemc,emcpt,gall,gemc,constituents_unsub,1);
 644         jet_sub->SetNumberOfTracks(nt);
 645         jet_sub->SetNumberOfClusters(nc);
 646         jet_sub->SortConstituents();
 647
 648         fastjet::PseudoJet area(fjw.GetJetAreaVector(ijet));
 649         jet_sub->SetArea(area.perp());
 650         jet_sub->SetAreaEta(area.eta());
 651         jet_sub->SetAreaPhi(area.phi());
 652         jet_sub->SetAreaEmc(area.perp());
 653         jetCount++;
 654       }
 655     }
 656 #endif
 657   } //constituent subtraction
 658 }
 659
 660 //________________________________________________________________________
 661 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
 662 {
 663   // Get the leading jets.
 664
 665   static Float_t pt[9999] = {0};
 666
 667   const Int_t n = (Int_t)array.size();
 668
 669   if (n < 1)
 670     return kFALSE;
 671
 672   for (Int_t i = 0; i < n; i++)
 673     pt[i] = array[i].perp();
 674
 675   TMath::Sort(n, pt, indexes);
 676
 677   return kTRUE;
 678 }
 679
 680 //________________________________________________________________________
 681 Bool_t AliEmcalJetTask::DoInit()
 682 {
 683   // Init. Return true if successful.
 684
 685   if (fTrackEfficiency < 1.) {
 686     if (gRandom) delete gRandom;
 687     gRandom = new TRandom3(0);
 688   }
 689
 690   // get input collections
 691   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 692
 693   // get the event
 694   fEvent = InputEvent();
 695   if (!fEvent) {
 696     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
 697     return 0;
 698   }
 699
 700   if (!(fEvent->FindListObject(fJetsSubName)) && fJetsSub)
 701     fEvent->AddObject(fJetsSub);
 702
 703   if (fTracksName == "Tracks")
 704     am->LoadBranch("Tracks");
 705   if (!fTracks && !fTracksName.IsNull()) {
 706     fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
 707     if (!fTracks) {
 708       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
 709       return 0;
 710     }
 711   }
 712   if (fTracks) {
 713     TClass cls(fTracks->GetClass()->GetName());
 714     if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
 715       fIsMcPart = 1;
 716   }
 717
 718   if (fCaloName == "CaloClusters")
 719     am->LoadBranch("CaloClusters");
 720   if (!fClus && !fCaloName.IsNull()) {
 721     fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
 722     if (!fClus) {
 723       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
 724       return 0;
 725     }
 726   }
 727   if (fClus) {
 728     TClass cls(fClus->GetClass()->GetName());
 729     if (cls.InheritsFrom("AliEmcalParticle"))
 730       fIsEmcPart = 1;
 731   }
 732
 733   if (!fRhoName.IsNull() && !fRhoParam) { // get rho from the event
 734     fRhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhoName));
 735     if (!fRhoParam) {
 736       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhoName.Data()));
 737       return 0;
 738     }
 739   }
 740   if (!fRhomName.IsNull() && !fRhomParam) { // get rhom from the event
 741     fRhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(fRhomName));
 742     if (!fRhomParam) {
 743       AliError(Form("%s: Could not retrieve rho %s!", GetName(), fRhomName.Data()));
 744       return 0;
 745     }
 746   }
 747
 748   // add jets to event if not yet there
 749   if (!(fEvent->FindListObject(fJetsName)))
 750     fEvent->AddObject(fJets);
 751   else {
 752     AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
 753     return 0;
 754   }
 755
 756   return 1;
 757 }
 758
 759 //___________________________________________________________________________________________________________________
 760 void  AliEmcalJetTask::FillJetConstituents(std::vector<fastjet::PseudoJet>& constituents,AliEmcalJet *jet,Double_t vertex[3],UInt_t jetCount,Int_t& nt,Int_t& nc,Double_t& maxCh,Double_t& maxNe,Int_t& ncharged,Int_t& nneutral,Double_t& neutralE,Double_t& mcpt,Int_t& cemc,Double_t& emcpt,Int_t& gall,Int_t& gemc,std::vector<fastjet::PseudoJet>& constituentsunsub,Int_t flagsub){
 761     nt            = 0;
 762     nc            = 0;
 763     neutralE   = 0;
 764     maxCh      = 0;
 765     maxNe      = 0;
 766     gall          = 0;
 767     gemc       =0;
 768     cemc          = 0;
 769     ncharged      = 0;
 770     nneutral      = 0;
 771     mcpt       = 0;
 772     emcpt      = 0;
 773     Int_t uid   = -1;
 774    AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 775     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
 776       if(flagsub==0) uid = constituents[ic].user_index();
 777       if(flagsub!=0) uid=GetIndexSub(constituents[ic].perp(),constituentsunsub);
 778       AliDebug(3,Form("Processing constituent %d", uid));
 779       if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
 780         ++gall;
 781         Double_t gphi = constituents[ic].phi();
 782         if (gphi<0)
 783           gphi += TMath::TwoPi();
 784         gphi *= TMath::RadToDeg();
 785         Double_t geta = constituents[ic].eta();
 786         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
 787             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
 788           ++gemc;
 789       } else if ((uid > 0) && fTracks) { // track constituent
 790         Int_t tid = uid - 100;
 791         AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
 792         if (!t) {
 793           AliError(Form("Could not find track %d",tid));
 794           continue;
 795         }
 796         if (jetCount < fMarkConst) {
 797           AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
 798           t->SetBit(fJetType);
 799         }
 800         Double_t cEta = t->Eta();
 801         Double_t cPhi = t->Phi();
 802         Double_t cPt  = t->Pt();
 803         Double_t cP   = t->P();
 804         if (t->Charge() == 0) {
 805           neutralE += cP;
 806           ++nneutral;
 807           if (cPt > maxNe)
 808             maxNe = cPt;
 809         } else {
 810           ++ncharged;
 811           if (cPt > maxCh)
 812             maxCh = cPt;
 813         }
 814         if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
 815           mcpt += cPt;
 816
 817         if (cPhi<0)
 818           cPhi += TMath::TwoPi();
 819         cPhi *= TMath::RadToDeg();
 820         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
 821             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
 822           emcpt += cPt;
 823           ++cemc;
 824         }
 825         // cout<<"Adding the track"<<endl;
 826         jet->AddTrackAt(tid, nt);
 827         ++nt;
 828       } else if (fClus) { // cluster constituent
 829         Int_t cid = -uid - 100;
 830         AliVCluster *c = 0;
 831         Double_t cEta=0,cPhi=0,cPt=0,cP=0;
 832         if (fIsEmcPart) {
 833           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
 834           if (!ep)
 835             continue;
 836           c = ep->GetCluster();
 837           if (!c)
 838             continue;
 839           if (jetCount < fMarkConst)
 840             ep->SetBit(fJetType);
 841           cEta = ep->Eta();
 842           cPhi = ep->Phi();
 843           cPt  = ep->Pt();
 844           cP   = ep->P();
 845         } else {
 846           c = static_cast<AliVCluster*>(fClus->At(cid));
 847           if (!c)
 848             continue;
 849           if (jetCount < fMarkConst)
 850             c->SetBit(fJetType);
 851           TLorentzVector nP;
 852           c->GetMomentum(nP, vertex);
 853           cEta = nP.Eta();
 854           cPhi = nP.Phi();
 855           cPt  = nP.Pt();
 856           cP   = nP.P();
 857         }
 858
 859         neutralE += cP;
 860         if (cPt > maxNe)
 861           maxNe = cPt;
 862
 863         if (c->GetLabel() > fMinMCLabel) // MC particle
 864           mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
 865
 866         if (cPhi<0)
 867           cPhi += TMath::TwoPi();
 868         cPhi *= TMath::RadToDeg();
 869         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
 870             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
 871           emcpt += cPt;
 872           ++cemc;
 873         }
 874
 875         jet->AddClusterAt(cid, nc);
 876         ++nc;
 877         ++nneutral;
 878       } else {
 879         AliError(Form("%s: No logical way to end up here.", GetName()));
 880         continue;
 881       }
 882     }
 883 }
 884 //______________________________________________________________________________________
 885 Int_t AliEmcalJetTask::GetIndexSub(Double_t ptsub,std::vector<fastjet::PseudoJet>& constituentsunsub){
 886
 887     for(UInt_t ii=0;ii<constituentsunsub.size();ii++){
 888     Double_t ptunsub=constituentsunsub[ii].perp();
 889     //cout<<ptunsub<<" "<<ptsub<<endl;
 890     if(ptsub==ptunsub) return constituentsunsub[ii].user_index();
 891
 892     } return -1;}
 893 //______________________________________________________________________________________