PWGJE/EMCALJetTasks/AliEmcalJetTask.cxx

   1 // $Id$
   2 //
   3 // Emcal jet finder task.
   4 //
   5 // Authors: C.Loizides, S.Aiola
   6
   7 #include <vector>
   8 #include "AliEmcalJetTask.h"
   9
  10 #include <TChain.h>
  11 #include <TClonesArray.h>
  12 #include <TList.h>
  13 #include <TLorentzVector.h>
  14 #include <TMath.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
  28 ClassImp(AliEmcalJetTask)
  29
  30 //________________________________________________________________________
  31 AliEmcalJetTask::AliEmcalJetTask() :
  32   AliAnalysisTaskSE("AliEmcalJetTask"),
  33   fTracksName("Tracks"),
  34   fCaloName("CaloClusters"),
  35   fJetsName("Jets"),
  36   fJetType(kNone),
  37   fConstSel(0),
  38   fMCConstSel(0),
  39   fMarkConst(kFALSE),
  40   fRadius(0.4),
  41   fMinJetTrackPt(0.15),
  42   fMinJetClusPt(0.15),
  43   fPhiMin(0),
  44   fPhiMax(TMath::TwoPi()),
  45   fEtaMin(-0.9),
  46   fEtaMax(+0.9),
  47   fMinJetArea(0.005),
  48   fMinJetPt(1.0),
  49   fJetPhiMin(-10),
  50   fJetPhiMax(+10),
  51   fJetEtaMin(-1),
  52   fJetEtaMax(+1),
  53   fGhostArea(0.005),
  54   fMinMCLabel(0),
  55   fIsInit(0),
  56   fIsPSelSet(0),
  57   fIsMcPart(0),
  58   fIsEmcPart(0),
  59   fJets(0),
  60   fEvent(0),
  61   fTracks(0),
  62   fClus(0)
  63 {
  64   // Default constructor.
  65 }
  66
  67 //________________________________________________________________________
  68 AliEmcalJetTask::AliEmcalJetTask(const char *name) :
  69   AliAnalysisTaskSE(name),
  70   fTracksName("Tracks"),
  71   fCaloName("CaloClusters"),
  72   fJetsName("Jets"),
  73   fJetType(kAKT|kFullJet|kRX1Jet),
  74   fConstSel(0),
  75   fMCConstSel(0),
  76   fMarkConst(kFALSE),
  77   fRadius(0.4),
  78   fMinJetTrackPt(0.15),
  79   fMinJetClusPt(0.15),
  80   fPhiMin(0),
  81   fPhiMax(TMath::TwoPi()),
  82   fEtaMin(-0.9),
  83   fEtaMax(+0.9),
  84   fMinJetArea(0.001),
  85   fMinJetPt(1.0),
  86   fJetPhiMin(-10),
  87   fJetPhiMax(+10),
  88   fJetEtaMin(-1),
  89   fJetEtaMax(+1),
  90   fGhostArea(0.005),
  91   fMinMCLabel(0),
  92   fIsInit(0),
  93   fIsPSelSet(0),
  94   fIsMcPart(0),
  95   fIsEmcPart(0),
  96   fJets(0),
  97   fEvent(0),
  98   fTracks(0),
  99   fClus(0)
 100 {
 101   // Standard constructor.
 102
 103   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
 104 }
 105
 106 //________________________________________________________________________
 107 AliEmcalJetTask::~AliEmcalJetTask()
 108 {
 109   // Destructor
 110 }
 111
 112 //________________________________________________________________________
 113 void AliEmcalJetTask::UserCreateOutputObjects()
 114 {
 115   // Create user objects.
 116
 117   fJets = new TClonesArray("AliEmcalJet");
 118   fJets->SetName(fJetsName);
 119 }
 120
 121 //________________________________________________________________________
 122 void AliEmcalJetTask::UserExec(Option_t *)
 123 {
 124   // Main loop, called for each event.
 125
 126   if (!fIsInit) {
 127     if (!DoInit())
 128       return;
 129     fIsInit = kTRUE;
 130   }
 131
 132   FindJets();
 133 }
 134
 135 //________________________________________________________________________
 136 void AliEmcalJetTask::Terminate(Option_t *)
 137 {
 138   // Called once at the end of the analysis.
 139 }
 140
 141 //________________________________________________________________________
 142 void AliEmcalJetTask::FindJets()
 143 {
 144   // Find jets.
 145
 146   if (!fTracks && !fClus)
 147     return;
 148
 149   TString name("kt");
 150   fastjet::JetAlgorithm jalgo(fastjet::kt_algorithm);
 151   if ((fJetType & kAKT) != 0) {
 152     name  = "antikt";
 153     jalgo = fastjet::antikt_algorithm;
 154     AliDebug(1,"Using AKT algorithm");
 155   }
 156   else {
 157     AliDebug(1,"Using KT algorithm");
 158   }
 159
 160   if ((fJetType & kR020Jet) != 0)
 161     fRadius = 0.2;
 162   else if ((fJetType & kR030Jet) != 0)
 163     fRadius = 0.3;
 164   else if ((fJetType & kR040Jet) != 0)
 165     fRadius = 0.4;
 166
 167   // setup fj wrapper
 168   AliFJWrapper fjw(name, name);
 169   fjw.SetAreaType(fastjet::active_area_explicit_ghosts);
 170   fjw.SetGhostArea(fGhostArea);
 171   fjw.SetR(fRadius);
 172   fjw.SetAlgorithm(jalgo);
 173   fjw.SetMaxRap(fEtaMax);
 174   fjw.Clear();
 175
 176   // get primary vertex
 177   Double_t vertex[3] = {0, 0, 0};
 178   fEvent->GetPrimaryVertex()->GetXYZ(vertex);
 179
 180   AliDebug(2,Form("Jet type = %d", fJetType));
 181
 182   if ((fIsMcPart || ((fJetType & kFullJet) != 0) || ((fJetType & kChargedJet) != 0)) && fTracks) {
 183     const Int_t Ntracks = fTracks->GetEntries();
 184     for (Int_t iTracks = 0; iTracks < Ntracks; ++iTracks) {
 185       AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(iTracks));
 186       if (!t)
 187         continue;
 188       if (fIsMcPart) {
 189         if (((fJetType & kChargedJet) != 0) && (t->Charge() == 0)) {
 190           AliDebug(2,Form("Skipping track %d because it is neutral.", iTracks));
 191           continue;
 192         }
 193         if (((fJetType & kNeutralJet) != 0) && (t->Charge() != 0)) {
 194           AliDebug(2,Form("Skipping track %d because it is charged.", iTracks));
 195           continue;
 196         }
 197       }
 198       if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) {
 199         if (t->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
 200           AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
 201           continue;
 202         }
 203         else {
 204           AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fMCConstSel, t->TestBits(TObject::kBitMask)));
 205         }
 206       }
 207       else {
 208         if (t->TestBits(fConstSel) != (Int_t)fConstSel) {
 209           AliDebug(2,Form("Skipping track %d because it does not match the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
 210           continue;
 211         }
 212         else {
 213           AliDebug(2,Form("Track %d matches the bit mask (%d, %d)", iTracks, fConstSel, t->TestBits(TObject::kBitMask)));
 214         }
 215       }
 216       if (t->Pt() < fMinJetTrackPt)
 217         continue;
 218       Double_t eta = t->Eta();
 219       Double_t phi = t->Phi();
 220       if ((eta<fEtaMin) || (eta>fEtaMax) ||
 221           (phi<fPhiMin) || (phi>fPhiMax))
 222         continue;
 223
 224       // offset of 100 for consistency with cluster ids
 225       AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", iTracks, t->GetLabel(), t->Pt()));
 226       fjw.AddInputVector(t->Px(), t->Py(), t->Pz(), t->P(), iTracks + 100);
 227     }
 228   }
 229
 230   if ((((fJetType & kFullJet) != 0) || ((fJetType & kNeutralJet) != 0)) && fClus) {
 231     const Int_t Nclus = fClus->GetEntries();
 232     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 233       AliVCluster *c = 0;
 234       Double_t cEta=0,cPhi=0,cPt=0;
 235       Double_t cPx=0,cPy=0,cPz=0;
 236       if (fIsEmcPart) {
 237         AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(iClus));
 238         if (!ep)
 239           continue;
 240
 241         c = ep->GetCluster();
 242         if (!c)
 243           continue;
 244
 245         if (c->GetLabel() > fMinMCLabel) {
 246           if (ep->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
 247             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
 248             continue;
 249           }
 250           else {
 251             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, ep->TestBits(TObject::kBitMask)));
 252           }
 253         }
 254         else {
 255           if (ep->TestBits(fConstSel) != (Int_t)fConstSel) {
 256             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
 257             continue;
 258           }
 259           else {
 260             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, ep->TestBits(TObject::kBitMask)));
 261           }
 262         }
 263
 264         cEta = ep->Eta();
 265         cPhi = ep->Phi();
 266         cPt  = ep->Pt();
 267         cPx  = ep->Px();
 268         cPy  = ep->Py();
 269         cPz  = ep->Pz();
 270       } else {
 271         c = static_cast<AliVCluster*>(fClus->At(iClus));
 272         if (!c)
 273           continue;
 274
 275         if (c->GetLabel() > fMinMCLabel) {
 276           if (c->TestBits(fMCConstSel) != (Int_t)fMCConstSel) {
 277             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
 278             continue;
 279           }
 280           else {
 281             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fMCConstSel, c->TestBits(TObject::kBitMask)));
 282           }
 283         }
 284         else {
 285           if (c->TestBits(fConstSel) != (Int_t)fConstSel) {
 286             AliDebug(2,Form("Skipping cluster %d because it does not match the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
 287             continue;
 288           }
 289           else {
 290             AliDebug(2,Form("Cluster %d matches the bit mask (%d, %d)", iClus, fConstSel, c->TestBits(TObject::kBitMask)));
 291           }
 292         }
 293
 294         TLorentzVector nP;
 295         c->GetMomentum(nP, vertex);
 296         cEta = nP.Eta();
 297         cPhi = nP.Phi();
 298         cPt  = nP.Pt();
 299         cPx  = nP.Px();
 300         cPy  = nP.Py();
 301         cPz  = nP.Pz();
 302       }
 303       if (!c->IsEMCAL())
 304         continue;
 305       if (cPt < fMinJetClusPt)
 306         continue;
 307       if ((cEta<fEtaMin) || (cEta>fEtaMax) ||
 308           (cPhi<fPhiMin) || (cPhi>fPhiMax))
 309         continue;
 310       // offset of 100 to skip ghost particles uid = -1
 311       AliDebug(2,Form("Cluster %d accepted (label = %d)", iClus, c->GetLabel()));
 312       fjw.AddInputVector(cPx, cPy, cPz, TMath::Sqrt(cPx*cPx+cPy*cPy+cPz*cPz), -iClus - 100);
 313     }
 314   }
 315
 316   // run jet finder
 317   fjw.Run();
 318
 319   // get geometry
 320   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
 321   if (!geom) {
 322     AliFatal(Form("%s: Can not create geometry", GetName()));
 323     return;
 324   }
 325
 326   // loop over fastjet jets
 327   std::vector<fastjet::PseudoJet> jets_incl = fjw.GetInclusiveJets();
 328   // sort jets according to jet pt
 329   static Int_t indexes[9999] = {-1};
 330   GetSortedArray(indexes, jets_incl);
 331
 332   AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
 333   for (UInt_t ijet=0, jetCount=0; ijet<jets_incl.size(); ++ijet) {
 334     Int_t ij = indexes[ijet];
 335     AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fjw.GetJetArea(ij)));
 336
 337     if (jets_incl[ij].perp()<fMinJetPt)
 338       continue;
 339     if (fjw.GetJetArea(ij)<fMinJetArea)
 340       continue;
 341     if ((jets_incl[ij].eta()<fJetEtaMin) || (jets_incl[ij].eta()>fJetEtaMax) ||
 342         (jets_incl[ij].phi()<fJetPhiMin) || (jets_incl[ij].phi()>fJetPhiMax))
 343       continue;
 344
 345     AliEmcalJet *jet = new ((*fJets)[jetCount])
 346       AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
 347
 348     // loop over constituents
 349     std::vector<fastjet::PseudoJet> constituents(fjw.GetJetConstituents(ij));
 350     jet->SetNumberOfTracks(constituents.size());
 351     jet->SetNumberOfClusters(constituents.size());
 352
 353     Int_t nt            = 0;
 354     Int_t nc            = 0;
 355     Double_t neutralE   = 0;
 356     Double_t maxCh      = 0;
 357     Double_t maxNe      = 0;
 358     Int_t gall          = 0;
 359     Int_t gemc          = 0;
 360     Int_t cemc          = 0;
 361     Int_t ncharged      = 0;
 362     Int_t nneutral      = 0;
 363     Double_t mcpt       = 0;
 364     Double_t emcpt      = 0;
 365
 366     for(UInt_t ic = 0; ic < constituents.size(); ++ic) {
 367       Int_t uid = constituents[ic].user_index();
 368       AliDebug(3,Form("Processing constituent %d", uid));
 369       if ((uid == -1) /*&& (constituents[ic].kt2() < 1e-25)*/) { //ghost particle
 370         ++gall;
 371         Double_t gphi = constituents[ic].phi();
 372         if (gphi<0)
 373           gphi += TMath::TwoPi();
 374         gphi *= TMath::RadToDeg();
 375         Double_t geta = constituents[ic].eta();
 376         if ((gphi > geom->GetArm1PhiMin()) && (gphi < geom->GetArm1PhiMax()) &&
 377             (geta > geom->GetArm1EtaMin()) && (geta < geom->GetArm1EtaMax()))
 378           ++gemc;
 379       } else if ((uid > 0) && fTracks) { // track constituent
 380         Int_t tid = uid - 100;
 381         AliVParticle *t = static_cast<AliVParticle*>(fTracks->At(tid));
 382         if (!t) {
 383           AliError(Form("Could not find track %d",tid));
 384           continue;
 385         }
 386         if (jetCount < fMarkConst) {
 387           AliDebug(2,Form("Marking track %d with bit map %d", tid, fJetType));
 388           t->SetBit(fJetType);
 389         }
 390         Double_t cEta = t->Eta();
 391         Double_t cPhi = t->Phi();
 392         Double_t cPt  = t->Pt();
 393         Double_t cP   = t->P();
 394         if (t->Charge() == 0) {
 395           neutralE += cP;
 396           ++nneutral;
 397           if (cPt > maxNe)
 398             maxNe = cPt;
 399         } else {
 400           ++ncharged;
 401           if (cPt > maxCh)
 402             maxCh = cPt;
 403         }
 404         if (fIsMcPart || TMath::Abs(t->GetLabel()) > fMinMCLabel) // check if MC particle
 405           mcpt += cPt;
 406
 407         if (cPhi<0)
 408           cPhi += TMath::TwoPi();
 409         cPhi *= TMath::RadToDeg();
 410         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
 411             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
 412           emcpt += cPt;
 413           ++cemc;
 414         }
 415
 416         jet->AddTrackAt(tid, nt);
 417         ++nt;
 418       } else if (fClus) { // cluster constituent
 419         Int_t cid = -uid - 100;
 420         AliVCluster *c = 0;
 421         Double_t cEta=0,cPhi=0,cPt=0,cP=0;
 422         if (fIsEmcPart) {
 423           AliEmcalParticle *ep = static_cast<AliEmcalParticle*>(fClus->At(cid));
 424           if (!ep)
 425             continue;
 426           c = ep->GetCluster();
 427           if (!c)
 428             continue;
 429           if (jetCount < fMarkConst)
 430             ep->SetBit(fJetType);
 431           cEta = ep->Eta();
 432           cPhi = ep->Phi();
 433           cPt  = ep->Pt();
 434           cP   = ep->P();
 435         } else {
 436           c = static_cast<AliVCluster*>(fClus->At(cid));
 437           if (!c)
 438             continue;
 439           if (jetCount < fMarkConst)
 440             c->SetBit(fJetType);
 441           TLorentzVector nP;
 442           c->GetMomentum(nP, vertex);
 443           cEta = nP.Eta();
 444           cPhi = nP.Phi();
 445           cPt  = nP.Pt();
 446           cP   = nP.P();
 447         }
 448
 449         neutralE += cP;
 450         if (cPt > maxNe)
 451           maxNe = cPt;
 452
 453         if (c->GetLabel() > fMinMCLabel) // MC particle
 454           mcpt += cPt;
 455
 456         if (cPhi<0)
 457           cPhi += TMath::TwoPi();
 458         cPhi *= TMath::RadToDeg();
 459         if ((cPhi > geom->GetArm1PhiMin()) && (cPhi < geom->GetArm1PhiMax()) &&
 460             (cEta > geom->GetArm1EtaMin()) && (cEta < geom->GetArm1EtaMax())) {
 461           emcpt += cPt;
 462           ++cemc;
 463         }
 464
 465         jet->AddClusterAt(cid, nc);
 466         ++nc;
 467         ++nneutral;
 468       } else {
 469         AliError(Form("%s: No logical way to end up here.", GetName()));
 470         continue;
 471       }
 472     } /* end of constituent loop */
 473
 474     jet->SetNumberOfTracks(nt);
 475     jet->SetNumberOfClusters(nc);
 476     jet->SortConstituents();
 477     jet->SetMaxNeutralPt(maxNe);
 478     jet->SetMaxChargedPt(maxCh);
 479     jet->SetNEF(neutralE / jet->E());
 480     fastjet::PseudoJet area(fjw.GetJetAreaVector(ij));
 481     jet->SetArea(area.perp());
 482     jet->SetAreaEta(area.eta());
 483     jet->SetAreaPhi(area.phi());
 484     jet->SetNumberOfCharged(ncharged);
 485     jet->SetNumberOfNeutrals(nneutral);
 486     jet->SetMCPt(mcpt);
 487     jet->SetNEmc(cemc);
 488     jet->SetPtEmc(emcpt);
 489
 490     if (gall > 0)
 491       jet->SetAreaEmc(fjw.GetJetArea(ij) * gemc / gall);
 492     else
 493       jet->SetAreaEmc(-1);
 494     if ((jet->Phi() > geom->GetArm1PhiMin() * TMath::DegToRad()) &&
 495         (jet->Phi() < geom->GetArm1PhiMax() * TMath::DegToRad()) &&
 496         (jet->Eta() > geom->GetArm1EtaMin()) &&
 497         (jet->Eta() < geom->GetArm1EtaMax()))
 498       jet->SetAxisInEmcal(kTRUE);
 499
 500     AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), (Int_t)constituents.size()));
 501     jetCount++;
 502   }
 503   //fJets->Sort();
 504 }
 505
 506 //________________________________________________________________________
 507 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
 508 {
 509   // Get the leading jets.
 510
 511   static Float_t pt[9999] = {0};
 512
 513   const Int_t n = (Int_t)array.size();
 514
 515   if (n < 1)
 516     return kFALSE;
 517
 518   for (Int_t i = 0; i < n; i++)
 519     pt[i] = array[i].perp();
 520
 521   TMath::Sort(n, pt, indexes);
 522
 523   return kTRUE;
 524 }
 525
 526 //________________________________________________________________________
 527 Bool_t AliEmcalJetTask::DoInit()
 528 {
 529   // Init. Return true if successful.
 530
 531   // get input collections
 532   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 533
 534   // get the event
 535   fEvent = InputEvent();
 536   if (!fEvent) {
 537     AliError(Form("%s: Could not retrieve event! Returning", GetName()));
 538     return 0;
 539   }
 540
 541   // add jets to event if not yet there
 542   fJets->Delete();
 543   if (!(fEvent->FindListObject(fJetsName)))
 544     fEvent->AddObject(fJets);
 545   else {
 546     AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
 547     return 0;
 548   }
 549
 550   if (fTracksName == "Tracks")
 551     am->LoadBranch("Tracks");
 552   if (!fTracks && !fTracksName.IsNull()) {
 553     fTracks = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fTracksName));
 554     if (!fTracks) {
 555       AliError(Form("%s: Pointer to tracks %s == 0", GetName(), fTracksName.Data()));
 556       return 0;
 557     }
 558   }
 559   if (fTracks) {
 560     TClass cls(fTracks->GetClass()->GetName());
 561     if (cls.InheritsFrom("AliMCParticle") || cls.InheritsFrom("AliAODMCParticle"))
 562       fIsMcPart = 1;
 563   }
 564
 565   if (fCaloName == "CaloClusters")
 566     am->LoadBranch("CaloClusters");
 567   if (!fClus && !fCaloName.IsNull()) {
 568     fClus = dynamic_cast<TClonesArray*>(fEvent->FindListObject(fCaloName));
 569     if (!fClus) {
 570       AliError(Form("%s: Pointer to clus %s == 0", GetName(), fCaloName.Data()));
 571       return 0;
 572     }
 573   }
 574   if (fClus) {
 575     TClass cls(fClus->GetClass()->GetName());
 576     if (cls.InheritsFrom("AliEmcalParticle"))
 577       fIsEmcPart = 1;
 578   }
 579
 580   return 1;
 581 }