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