AliPhysics  05c4c93 (05c4c93)
AliEmcalJetTask.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2016, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 #include <vector>
17 
18 #include <TClonesArray.h>
19 #include <TMath.h>
20 #include <TRandom3.h>
21 
22 #include <AliVCluster.h>
23 #include <AliVEvent.h>
24 #include <AliVParticle.h>
25 #include <AliEMCALGeometry.h>
26 
27 #include <AliAnalysisManager.h>
28 #include <AliVEventHandler.h>
29 #include <AliMultiInputEventHandler.h>
30 
31 #include "AliTLorentzVector.h"
32 #include "AliEmcalJet.h"
33 #include "AliEmcalParticle.h"
34 #include "AliFJWrapper.h"
35 #include "AliEmcalJetUtility.h"
36 #include "AliParticleContainer.h"
37 #include "AliClusterContainer.h"
40 
41 #include "AliEmcalJetTask.h"
42 
43 using std::cout;
44 using std::endl;
45 using std::cerr;
46 
48 ClassImp(AliEmcalJetTask);
50 
52 
59  fJetsTag(),
60  fJetType(AliJetContainer::kFullJet),
61  fJetAlgo(AliJetContainer::antikt_algorithm),
62  fRecombScheme(AliJetContainer::pt_scheme),
63  fRadius(0.4),
64  fMinJetArea(0.001),
65  fMinJetPt(1.0),
66  fJetPhiMin(-10),
67  fJetPhiMax(+10),
68  fJetEtaMin(-1),
69  fJetEtaMax(+1),
70  fGhostArea(0.005),
71  fTrackEfficiency(1.),
72  fUtilities(0),
73  fTrackEfficiencyOnlyForEmbedding(kFALSE),
74  fLocked(0),
75  fFillConstituents(kTRUE),
76  fJetsName(),
77  fIsInit(0),
78  fIsPSelSet(0),
79  fIsEmcPart(0),
80  fLegacyMode(kFALSE),
81  fFillGhost(kFALSE),
82  fJets(0),
83  fFastJetWrapper("AliEmcalJetTask","AliEmcalJetTask"),
84  fClusterContainerIndexMap(),
85  fParticleContainerIndexMap()
86 {
87 }
88 
95  fJetsTag("Jets"),
96  fJetType(AliJetContainer::kFullJet),
97  fJetAlgo(AliJetContainer::antikt_algorithm),
98  fRecombScheme(AliJetContainer::pt_scheme),
99  fRadius(0.4),
100  fMinJetArea(0.001),
101  fMinJetPt(1.0),
102  fJetPhiMin(-10),
103  fJetPhiMax(+10),
104  fJetEtaMin(-1),
105  fJetEtaMax(+1),
106  fGhostArea(0.005),
107  fTrackEfficiency(1.),
108  fUtilities(0),
110  fLocked(0),
111  fFillConstituents(kTRUE),
112  fJetsName(),
113  fIsInit(0),
114  fIsPSelSet(0),
115  fIsEmcPart(0),
116  fLegacyMode(kFALSE),
117  fFillGhost(kFALSE),
118  fJets(0),
119  fFastJetWrapper(name,name),
122 {
123 }
124 
129 {
130 }
131 
137 {
138  if (!fUtilities) fUtilities = new TObjArray();
139  if (fUtilities->FindObject(utility)) {
140  Error("AddUtility", "Jet utility %s already connected.", utility->GetName());
141  return utility;
142  }
143  fUtilities->Add(utility);
144  utility->SetJetTask(this);
145 
146  return utility;
147 }
148 
154 {
155  TIter next(fUtilities);
156  AliEmcalJetUtility *utility = 0;
157  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Init();
158 }
164 {
165  TIter next(fUtilities);
166  AliEmcalJetUtility *utility = 0;
167  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->InitEvent(fFastJetWrapper);
168 }
169 
176 {
177  TIter next(fUtilities);
178  AliEmcalJetUtility *utility = 0;
179  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Prepare(fFastJetWrapper);
180 }
181 
187 {
188  TIter next(fUtilities);
189  AliEmcalJetUtility *utility = 0;
190  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->ProcessJet(jet, ij, fFastJetWrapper);
191 }
192 
198 {
199  TIter next(fUtilities);
200  AliEmcalJetUtility *utility = 0;
201  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Terminate(fFastJetWrapper);
202 }
203 
209 {
210  InitEvent();
211  // clear the jet array (normally a null operation)
212  fJets->Delete();
213  Int_t n = FindJets();
214 
215  if (n == 0) return kFALSE;
216 
217  FillJetBranch();
218 
219  return kTRUE;
220 }
221 
230 {
231  if (fParticleCollArray.GetEntriesFast() == 0 && fClusterCollArray.GetEntriesFast() == 0){
232  AliError("No tracks or clusters, returning.");
233  return 0;
234  }
235 
237 
238  AliDebug(2,Form("Jet type = %d", fJetType));
239 
240  Int_t iColl = 1;
241  TIter nextPartColl(&fParticleCollArray);
242  AliParticleContainer* tracks = 0;
243  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
244  AliDebug(2,Form("Tracks from collection %d: '%s'. Embedded: %i, nTracks: %i", iColl-1, tracks->GetName(), tracks->GetIsEmbedding(), tracks->GetNParticles()));
246  for (AliParticleIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
247  // artificial inefficiency
248  if (fTrackEfficiency < 1.) {
249  if (fTrackEfficiencyOnlyForEmbedding == kFALSE || (fTrackEfficiencyOnlyForEmbedding == kTRUE && tracks->GetIsEmbedding())) {
250  Double_t rnd = gRandom->Rndm();
251  if (fTrackEfficiency < rnd) {
252  AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", it.current_index()));
253  continue;
254  }
255  }
256  }
257 
258  AliDebug(2,Form("Track %d accepted (label = %d, pt = %f, eta = %f, phi = %f, E = %f, m = %f, px = %f, py = %f, pz = %f)", it.current_index(), it->second->GetLabel(), it->first.Pt(), it->first.Eta(), it->first.Phi(), it->first.E(), it->first.M(), it->first.Px(), it->first.Py(), it->first.Pz()));
259  Int_t uid = it.current_index() + fgkConstIndexShift * iColl;
260  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
261  }
262  iColl++;
263  }
264 
265  iColl = 1;
266  TIter nextClusColl(&fClusterCollArray);
267  AliClusterContainer* clusters = 0;
268  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
269  AliDebug(2,Form("Clusters from collection %d: '%s'. Embedded: %i, nClusters: %i", iColl-1, clusters->GetName(), clusters->GetIsEmbedding(), clusters->GetNClusters()));
271  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
272  AliDebug(2,Form("Cluster %d accepted (label = %d, energy = %.3f)", it.current_index(), it->second->GetLabel(), it->first.E()));
273  Int_t uid = -it.current_index() - fgkConstIndexShift * iColl;
274  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
275  }
276  iColl++;
277  }
278 
279  if (fFastJetWrapper.GetInputVectors().size() == 0) return 0;
280 
281  // run jet finder
283 
284  return fFastJetWrapper.GetInclusiveJets().size();
285 }
286 
293 {
295 
296  // loop over fastjet jets
297  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper.GetInclusiveJets();
298  // sort jets according to jet pt
299  static Int_t indexes[9999] = {-1};
300  GetSortedArray(indexes, jets_incl);
301 
302  AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
303  for (UInt_t ijet = 0, jetCount = 0; ijet < jets_incl.size(); ++ijet) {
304  Int_t ij = indexes[ijet];
305  AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fFastJetWrapper.GetJetArea(ij)));
306 
307  if (jets_incl[ij].perp() < fMinJetPt) continue;
308  if (fFastJetWrapper.GetJetArea(ij) < fMinJetArea) continue;
309  if ((jets_incl[ij].eta() < fJetEtaMin) || (jets_incl[ij].eta() > fJetEtaMax) ||
310  (jets_incl[ij].phi() < fJetPhiMin) || (jets_incl[ij].phi() > fJetPhiMax))
311  continue;
312 
313  AliEmcalJet *jet = new ((*fJets)[jetCount])
314  AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
315  jet->SetLabel(ij);
316 
317  fastjet::PseudoJet area(fFastJetWrapper.GetJetAreaVector(ij));
318  jet->SetArea(area.perp());
319  jet->SetAreaEta(area.eta());
320  jet->SetAreaPhi(area.phi());
321  jet->SetAreaE(area.E());
323 
324  // Fill constituent info
325  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper.GetJetConstituents(ij));
326  FillJetConstituents(jet, constituents, constituents);
327 
328  if (fGeom) {
329  if ((jet->Phi() > fGeom->GetArm1PhiMin() * TMath::DegToRad()) &&
330  (jet->Phi() < fGeom->GetArm1PhiMax() * TMath::DegToRad()) &&
331  (jet->Eta() > fGeom->GetArm1EtaMin()) &&
332  (jet->Eta() < fGeom->GetArm1EtaMax()))
333  jet->SetAxisInEmcal(kTRUE);
334  }
335 
336  ExecuteUtilities(jet, ij);
337 
338  AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), jet->GetNumberOfConstituents()));
339  jetCount++;
340  }
341 
343 }
344 
351 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
352 {
353  static Float_t pt[9999] = {0};
354 
355  const Int_t n = (Int_t)array.size();
356 
357  if (n < 1)
358  return kFALSE;
359 
360  for (Int_t i = 0; i < n; i++)
361  pt[i] = array[i].perp();
362 
363  TMath::Sort(n, pt, indexes);
364 
365  return kTRUE;
366 }
367 
374 {
375  if (fTrackEfficiency < 1.) {
376  if (gRandom) delete gRandom;
377  gRandom = new TRandom3(0);
378  }
379 
381  std::cout << GetName() << ": Name of the jet container: " << fJetsName << std::endl;
382  std::cout << "Use this name in order to connect jet containers in your task to connect to the collection of jets found by this jet finder" << std::endl;
383  if(auto partcont = GetParticleContainer(0)) {
384  std::cout << "Found particle container with name " << partcont->GetName() << std::endl;
385  } else {
386  std::cout << "Not particle container found for task" << std::endl;
387  }
388  if(auto clustcont = GetClusterContainer(0)){
389  std::cout << "Found cluster container with name " << clustcont->GetName() << std::endl;
390  } else {
391  std::cout << "Not cluster container found for task" << std::endl;
392  }
393 
394  // add jets to event if not yet there
395  if (!(InputEvent()->FindListObject(fJetsName))) {
396  fJets = new TClonesArray("AliEmcalJet");
397  fJets->SetName(fJetsName);
398  ::Info("AliEmcalJetTask::ExecOnce", "Jet collection with name '%s' has been added to the event.", fJetsName.Data());
399  InputEvent()->AddObject(fJets);
400  }
401  else {
402  AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
403  return;
404  }
405 
406  // setup fj wrapper
407  fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts);
413 
414 
415  // setting legacy mode
416  if (fLegacyMode) {
418  }
419 
420  InitUtilities();
421 
423 
424  // Setup container utils. Must be called after AliAnalysisTaskEmcal::ExecOnce() so that the
425  // containers' arrays are setup.
428 }
429 
439 void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents,
440  std::vector<fastjet::PseudoJet>& constituents_unsub, Int_t flag, TString particlesSubName)
441 {
442  Int_t nt = 0;
443  Int_t nc = 0;
444  Double_t neutralE = 0.;
445  Double_t maxCh = 0.;
446  Double_t maxNe = 0.;
447  Int_t gall = 0;
448  Int_t gemc = 0;
449  Int_t cemc = 0;
450  Int_t ncharged = 0;
451  Int_t nneutral = 0;
452  Double_t mcpt = 0.;
453  Double_t emcpt = 0.;
454  TClonesArray * particles_sub = 0;
455 
456  Int_t uid = -1;
457 
458  jet->SetNumberOfTracks(constituents.size());
459  jet->SetNumberOfClusters(constituents.size());
460 
461  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
462 
463  if (flag == 0) {
464  uid = constituents[ic].user_index();
465  }
466  else {
467  if (constituents[ic].perp()<1.e-10) {
468  uid=-1;
469  }
470  else {
471  uid = constituents[ic].user_index();
472  }
473  if (uid==0) {
474  AliError("correspondence between un/subtracted constituent not found");
475  continue;
476  }
477  }
478 
479  AliDebug(3,Form("Processing constituent %d", uid));
480  if (uid == -1) { //ghost particle
481  ++gall;
482  if (fGeom) {
483  Double_t gphi = constituents[ic].phi();
484  if (gphi < 0) gphi += TMath::TwoPi();
485  gphi *= TMath::RadToDeg();
486  Double_t geta = constituents[ic].eta();
487  if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) &&
488  (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax()))
489  ++gemc;
490  }
491 
492  if (fFillGhost) jet->AddGhost(constituents[ic].px(),
493  constituents[ic].py(),
494  constituents[ic].pz(),
495  constituents[ic].e());
496  }
497  else if (uid >= fgkConstIndexShift) { // track constituent
498  Int_t iColl = uid / fgkConstIndexShift;
499  Int_t tid = uid - iColl * fgkConstIndexShift;
500  iColl--;
501  AliDebug(3,Form("Constituent %d is a track from collection %d and with ID %d", uid, iColl, tid));
502  AliParticleContainer* partCont = GetParticleContainer(iColl);
503  if (!partCont) {
504  AliError(Form("Could not find particle container %d",iColl));
505  continue;
506  }
507  AliVParticle *t = partCont->GetParticle(tid);
508  if (!t) {
509  AliError(Form("Could not find track %d",tid));
510  continue;
511  }
512  if(fFillConstituents){
513  jet->AddParticleConstituent(t, partCont->GetIsEmbedding(), fParticleContainerIndexMap.GlobalIndexFromLocalIndex(partCont, tid));
514  }
515 
516  Double_t cEta = t->Eta();
517  Double_t cPhi = t->Phi();
518  Double_t cPt = t->Pt();
519  Double_t cP = t->P();
520  if (t->Charge() == 0) {
521  neutralE += cP;
522  ++nneutral;
523  if (cPt > maxNe) maxNe = cPt;
524  } else {
525  ++ncharged;
526  if (cPt > maxCh) maxCh = cPt;
527  }
528 
529  // check if MC particle
530  if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt;
531 
532  if (fGeom) {
533  if (cPhi < 0) cPhi += TMath::TwoPi();
534  cPhi *= TMath::RadToDeg();
535  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
536  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
537  emcpt += cPt;
538  ++cemc;
539  }
540  }
541 
542  if (flag == 0 || particlesSubName == "") {
544  }
545  else {
546  // Get the particle container and array corresponding to the subtracted particles
547  partCont = GetParticleContainer(particlesSubName);
548  particles_sub = partCont->GetArray();
549  // Create the new particle in the particles_sub array and add it to the jet
550  Int_t part_sub_id = particles_sub->GetEntriesFast();
551  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(dynamic_cast<AliVTrack*>(t)); // SA: probably need to be fixed!!
552  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
553  jet->AddTrackAt(fParticleContainerIndexMap.GlobalIndexFromLocalIndex(partCont, part_sub_id), nt);
554  }
555 
556  ++nt;
557  }
558  else if (uid <= -fgkConstIndexShift) { // cluster constituent
559  Int_t iColl = -uid / fgkConstIndexShift;
560  Int_t cid = -uid - iColl * fgkConstIndexShift;
561  iColl--;
562  AliDebug(3,Form("Constituent %d is a cluster from collection %d and with ID %d", uid, iColl, cid));
563  AliClusterContainer* clusCont = GetClusterContainer(iColl);
564  AliVCluster *c = clusCont->GetCluster(cid);
565 
566  if (!c) continue;
567 
569  clusCont->GetMomentum(nP, cid);
570 
571  Double_t cEta = nP.Eta();
572  Double_t cPhi = nP.Phi_0_2pi();
573  Double_t cPt = nP.Pt();
574  Double_t cP = nP.P();
575  Double_t pvec[3] = {nP.Px(), nP.Py(), nP.Pz()};
576  if(fFillConstituents) jet->AddClusterConstituent(c, (AliVCluster::VCluUserDefEnergy_t)clusCont->GetDefaultClusterEnergy(), pvec, clusCont->GetIsEmbedding(), fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, cid));
577 
578  neutralE += cP;
579  if (cPt > maxNe) maxNe = cPt;
580 
581  // MC particle
582  if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
583 
584  if (fGeom) {
585  if (cPhi<0) cPhi += TMath::TwoPi();
586  cPhi *= TMath::RadToDeg();
587  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
588  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
589  emcpt += cPt;
590  ++cemc;
591  }
592  }
593 
594  if (flag == 0 || particlesSubName == "") {
596  }
597  else {
598  // Get the cluster container and array corresponding to the subtracted particles
599  clusCont = GetClusterContainer(particlesSubName);
600  particles_sub = clusCont->GetArray();
601  // Create the new particle in the particles_sub array and add it to the jet
602  Int_t part_sub_id = particles_sub->GetEntriesFast();
603  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c);
604  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
605  jet->AddClusterAt(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, part_sub_id), nc);
606  }
607 
608  ++nc;
609  ++nneutral;
610  }
611  else {
612  AliError(Form("%s: No logical way to end up here (uid = %d).", GetName(), uid));
613  continue;
614  }
615  }
616 
617  jet->SetNumberOfTracks(nt);
618  jet->SetNumberOfClusters(nc);
619  jet->SetNEF(neutralE / jet->E());
620  jet->SetMaxChargedPt(maxCh);
621  jet->SetMaxNeutralPt(maxNe);
622  if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall);
623  else jet->SetAreaEmc(-1);
624  jet->SetNEmc(cemc);
625  jet->SetNumberOfCharged(ncharged);
626  jet->SetNumberOfNeutrals(nneutral);
627  jet->SetMCPt(mcpt);
628  jet->SetPtEmc(emcpt);
629  jet->SortConstituents();
630 }
631 
639 {
640  if (fLocked) {
641  AliFatal("Jet finder task is locked! Changing properties is not allowed.");
642  return kTRUE;
643  }
644  else {
645  return kFALSE;
646  }
647 }
648 
656 {
657  if(!fIsPSelSet) {
658  fIsPSelSet = kTRUE;
659  fOfflineTriggerMask = offlineTriggerMask;
660  }
661  else {
662  AliWarning("Manually setting the event selection for jet finders is not allowed! Using trigger=old_trigger | your_trigger");
663  fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask;
664  }
665 }
666 
673 {
674  if (IsLocked()) return;
675 
676  TIter nextPartColl(&fParticleCollArray);
677  AliParticleContainer* tracks = 0;
678  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
679  tracks->SetParticleEtaLimits(emi, ema);
680  }
681 }
682 
688 {
689  if (IsLocked()) return;
690 
691  TIter nextClusColl(&fClusterCollArray);
692  AliClusterContainer* clusters = 0;
693  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
694  clusters->SetClusPtCut(min);
695  }
696 }
697 
703 {
704  if (IsLocked()) return;
705 
706  TIter nextClusColl(&fClusterCollArray);
707  AliClusterContainer* clusters = 0;
708  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
709  clusters->SetClusECut(min);
710  }
711 }
712 
718 {
719  if (IsLocked()) return;
720 
721  TIter nextPartColl(&fParticleCollArray);
722  AliParticleContainer* tracks = 0;
723  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
724  tracks->SetParticlePtCut(min);
725  }
726 }
727 
734 {
735  if (IsLocked()) return;
736 
737  TIter nextPartColl(&fParticleCollArray);
738  AliParticleContainer* tracks = 0;
739  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
740  tracks->SetParticlePhiLimits(pmi, pma);
741  }
742 }
743 
750 fastjet::JetAlgorithm AliEmcalJetTask::ConvertToFJAlgo(EJetAlgo_t algo)
751 {
752  switch(algo) {
754  return fastjet::kt_algorithm;
756  return fastjet::antikt_algorithm;
758  return fastjet::cambridge_algorithm;
760  return fastjet::genkt_algorithm;
762  return fastjet::cambridge_for_passive_algorithm;
764  return fastjet::genkt_for_passive_algorithm;
766  return fastjet::plugin_algorithm;
768  return fastjet::undefined_jet_algorithm;
769 
770  default:
771  ::Error("AliEmcalJetTask::ConvertToFJAlgo", "Jet algorithm %d not recognized!!!", algo);
772  return fastjet::undefined_jet_algorithm;
773  }
774 }
775 
782 fastjet::RecombinationScheme AliEmcalJetTask::ConvertToFJRecoScheme(ERecoScheme_t reco)
783 {
784  switch(reco) {
786  return fastjet::E_scheme;
787 
789  return fastjet::pt_scheme;
790 
792  return fastjet::pt2_scheme;
793 
795  return fastjet::Et_scheme;
796 
798  return fastjet::Et2_scheme;
799 
801  return fastjet::BIpt_scheme;
802 
804  return fastjet::BIpt2_scheme;
805 
807  return fastjet::external_scheme;
808 
809  default:
810  ::Error("AliEmcalJetTask::ConvertToFJRecoScheme", "Recombination scheme %d not recognized!!!", reco);
811  return fastjet::external_scheme;
812  }
813 }
814 
820 
821  //This method has to be called after the run number is known because it needs the EMCal geometry object.
822 
823  UInt_t jetAcceptanceType = AliEmcalJet::kUser; // all jets satify the "no acceptance cut" condition
824 
825  // Check if TPC
826  if( eta < 0.9 && eta > -0.9 ) {
827  jetAcceptanceType |= AliEmcalJet::kTPC;
828  // Check if TPCfid
829  if (eta < 0.9 - r && eta > -0.9 + r)
830  jetAcceptanceType |= AliEmcalJet::kTPCfid;
831  }
832 
833  // Check if EMCAL
834  if( IsJetInEmcal(eta, phi, 0) ) {
835  jetAcceptanceType |= AliEmcalJet::kEMCAL;
836  // Check if EMCALfid
837  if( IsJetInEmcal(eta, phi, r) )
838  jetAcceptanceType |= AliEmcalJet::kEMCALfid;
839  }
840 
841  // Check if DCAL (i.e. eta-phi rectangle spanning DCal, which includes most of PHOS)
842  if( IsJetInDcal(eta, phi, 0) ) {
843  jetAcceptanceType |= AliEmcalJet::kDCAL;
844  // Check if DCALfid
845  if( IsJetInDcal(eta, phi, r) )
846  jetAcceptanceType |= AliEmcalJet::kDCALfid;
847  }
848 
849  // Check if DCALonly (i.e. ONLY DCal, does not include any of PHOS region)
850  if( IsJetInDcalOnly(eta, phi, 0) ) {
851  jetAcceptanceType |= AliEmcalJet::kDCALonly;
852  // Check if DCALonlyfid
853  if( IsJetInDcalOnly(eta, phi, r) )
854  jetAcceptanceType |= AliEmcalJet::kDCALonlyfid;
855  }
856 
857  // Check if PHOS
858  if( IsJetInPhos(eta, phi, 0) ) {
859  jetAcceptanceType |= AliEmcalJet::kPHOS;
860  // Check if PHOSfid
861  if( IsJetInPhos(eta, phi, r) )
862  jetAcceptanceType |= AliEmcalJet::kPHOSfid;
863  }
864 
865  return jetAcceptanceType;
866 }
867 
872 {
873  if (!fGeom) return kFALSE;
874 
875  if ( eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
876  if(fRunNumber >= 177295 && fRunNumber <= 197470) {//small SM masked in 2012 and 2013
877  if ( phi < 3.135 - r && phi > 1.405 + r )
878  return kTRUE;
879  }
880  else {
881  if ( phi < fGeom->GetEMCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetArm1PhiMin() * TMath::DegToRad() + r)
882  return kTRUE;
883  }
884  }
885 
886  return kFALSE;
887 }
888 
893 {
894  if (!fGeom) return kFALSE;
895  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
896  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r)
897  return kTRUE;
898  }
899  return kFALSE;
900 }
901 
909 {
910  if (!fGeom) return kFALSE;
911 
912  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r) {
913  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r) {
914 
915  if (TMath::Abs(eta) > fGeom->GetDCALInnerExtandedEta() + r) {
916  return kTRUE;
917  }
918  if (r < 1e-6) {
919  if (phi > fGeom->GetEMCGeometry()->GetDCALStandardPhiMax() * TMath::DegToRad())
920  return kTRUE;
921  }
922 
923  }
924  }
925 
926  return kFALSE;
927 }
928 
933 {
934  Double_t etaMax = 0.130;
935  Double_t etaMin = -0.130;
936  Double_t phiMax = 320;
937  Double_t phiMin = 260; // Run 1
938  if (fRunNumber > 209121)
939  phiMin = 250; // Run 2
940 
941  if (eta < etaMax - r && eta > etaMin + r ) {
942  if (phi < phiMax * TMath::DegToRad() - r && phi > phiMin * TMath::DegToRad() + r)
943  return kTRUE;
944  }
945  return kFALSE;
946 }
947 
966  const TString nTracks, const TString nClusters,
967  const AliJetContainer::EJetAlgo_t jetAlgo, const Double_t radius, const AliJetContainer::EJetType_t jetType,
968  const Double_t minTrPt, const Double_t minClPt,
969  const Double_t ghostArea, const AliJetContainer::ERecoScheme_t reco,
970  const TString tag, const Double_t minJetPt,
971  const Bool_t lockTask, const Bool_t bFillGhosts
972 )
973 {
974  // Get the pointer to the existing analysis manager via the static access method
975  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
976  if (!mgr) {
977  ::Error("AddTaskEmcalJet", "No analysis manager to connect to.");
978  return 0;
979  }
980 
981  // Check the analysis type using the event handlers connected to the analysis manager
982  AliVEventHandler* handler = mgr->GetInputEventHandler();
983  if (!handler) {
984  ::Error("AddTaskEmcalJet", "This task requires an input event handler");
985  return 0;
986  }
987 
988  EDataType_t dataType = kUnknownDataType;
989 
990  if (handler->InheritsFrom("AliESDInputHandler")) {
991  dataType = kESD;
992  }
993  else if (handler->InheritsFrom("AliAODInputHandler")) {
994  dataType = kAOD;
995  }
996 
997  //-------------------------------------------------------
998  // Init the task and do settings
999  //-------------------------------------------------------
1000 
1001  TString trackName(nTracks);
1002  TString clusName(nClusters);
1003 
1004  if (trackName == "usedefault") {
1005  if (dataType == kESD) {
1006  trackName = "Tracks";
1007  }
1008  else if (dataType == kAOD) {
1009  trackName = "tracks";
1010  }
1011  else {
1012  trackName = "";
1013  }
1014  }
1015 
1016  if (clusName == "usedefault") {
1017  if (dataType == kESD) {
1018  clusName = "CaloClusters";
1019  }
1020  else if (dataType == kAOD) {
1021  clusName = "caloClusters";
1022  }
1023  else {
1024  clusName = "";
1025  }
1026  }
1027 
1028  AliParticleContainer* partCont = 0;
1029  if (trackName == "mcparticles") {
1030  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
1031  partCont = mcpartCont;
1032  }
1033  else if (trackName == "tracks" || trackName == "Tracks") {
1034  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
1035  partCont = trackCont;
1036  }
1037  else if (!trackName.IsNull()) {
1038  partCont = new AliParticleContainer(trackName);
1039  }
1040  if (partCont) partCont->SetParticlePtCut(minTrPt);
1041 
1042  AliClusterContainer* clusCont = 0;
1043  if (!clusName.IsNull()) {
1044  clusCont = new AliClusterContainer(clusName);
1045  clusCont->SetClusECut(0.);
1046  clusCont->SetClusPtCut(0.);
1047  clusCont->SetClusHadCorrEnergyCut(minClPt);
1048  clusCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
1049  }
1050 
1051  switch (jetType) {
1053  if (partCont) partCont->SetCharge(AliParticleContainer::kCharged);
1054  break;
1056  if (partCont) partCont->SetCharge(AliParticleContainer::kNeutral);
1057  break;
1058  default:
1059  break;
1060  }
1061 
1062  TString name = AliJetContainer::GenerateJetName(jetType, jetAlgo, reco, radius, partCont, clusCont, tag);
1063 
1064  Printf("Jet task name: %s", name.Data());
1065 
1066  AliEmcalJetTask* mgrTask = static_cast<AliEmcalJetTask *>(mgr->GetTask(name.Data()));
1067  if (mgrTask) return mgrTask;
1068 
1069  AliEmcalJetTask* jetTask = new AliEmcalJetTask(name);
1070  jetTask->SetJetType(jetType);
1071  jetTask->SetJetAlgo(jetAlgo);
1072  jetTask->SetRecombScheme(reco);
1073  jetTask->SetRadius(radius);
1074  if (partCont) jetTask->AdoptParticleContainer(partCont);
1075  if (clusCont) jetTask->AdoptClusterContainer(clusCont);
1076  jetTask->SetJetsName(tag);
1077  jetTask->SetMinJetPt(minJetPt);
1078  jetTask->SetGhostArea(ghostArea);
1079 
1080  if (bFillGhosts) jetTask->SetFillGhost();
1081  if (lockTask) jetTask->SetLocked();
1082 
1083  // Final settings, pass to manager and set the containers
1084 
1085  mgr->AddTask(jetTask);
1086 
1087  // Create containers for input/output
1088  AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer();
1089  mgr->ConnectInput(jetTask, 0, cinput);
1090 
1091  return jetTask;
1092 }
Bool_t fTrackEfficiencyOnlyForEmbedding
tituent Apply aritificial tracking inefficiency only for embedded tracks
void SetMaxNeutralPt(Double32_t t)
Definition: AliEmcalJet.h:262
void SetJetsName(const char *n)
void SetRecombScheme(ERecoScheme_t scheme)
fastjet::PseudoJet GetJetAreaVector(UInt_t idx) const
Double_t Area() const
Definition: AliEmcalJet.h:130
TObjArray fClusterCollArray
cluster collection array
void AddTrackAt(Int_t track, Int_t idx)
Definition: AliEmcalJet.h:275
void SetParticlePtCut(Double_t cut)
TClonesArray * fJets
!jet collection
virtual Int_t Run()
double Double_t
Definition: External.C:58
Bool_t fIsEmcPart
!=true if emcal particles are given as input (for clusters)
void SetEtaRange(Double_t emi, Double_t ema)
Bool_t IsLocked() const
Bool_t fIsInit
!=true if already initialized
PHOS acceptance.
Definition: AliEmcalJet.h:75
DCal acceptance – spans entire rectangular region in eta-phi (including most of PHOS) ...
Definition: AliEmcalJet.h:71
void AdoptParticleContainer(AliParticleContainer *cont)
Double_t Eta() const
Definition: AliEmcalJet.h:121
const AliClusterIterableMomentumContainer accepted_momentum() const
Base task in the EMCAL framework.
bidirectional stl iterator over the EMCAL iterable container
Double_t fJetEtaMin
minimum eta to keep jet in output
EJetType_t fJetType
jet type (full, charged, neutral)
Bool_t IsJetInEmcal(Double_t eta, Double_t phi, Double_t r)
Double_t Phi() const
Definition: AliEmcalJet.h:117
Int_t GetNParticles() const
Container with name, TClonesArray and cuts for particles.
void AddClusterAt(Int_t clus, Int_t idx)
Definition: AliEmcalJet.h:274
AliEmcalJetUtility * AddUtility(AliEmcalJetUtility *utility)
Bool_t fFillGhost
=true ghost particles will be filled in AliEmcalJet obj
Declaration of class AliTLorentzVector.
virtual void ProcessJet(AliEmcalJet *jet, Int_t ij, AliFJWrapper &fjw)=0
Double_t phiMin
void SetArea(Double_t a)
Definition: AliEmcalJet.h:256
void SetJetType(EJetType_t t)
void ExecuteUtilities(AliEmcalJet *jet, Int_t ij)
Full acceptance, i.e. no acceptance cut applied – left to user.
Definition: AliEmcalJet.h:77
TCanvas * c
Definition: TestFitELoss.C:172
void SetMaxRap(Double_t maxrap)
Definition: AliFJWrapper.h:126
virtual void Terminate(AliFJWrapper &fjw)=0
void SetMinJetPt(Double_t j)
void SetAreaE(Double_t a)
Definition: AliEmcalJet.h:259
AliEmcalContainerIndexMap< AliParticleContainer, AliVParticle > fParticleContainerIndexMap
! Mapping between index and particle containers
ERecoScheme_t fRecombScheme
recombination scheme used by fastjet
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
Bool_t fFillConstituents
If true jet consituents will be filled to the AliEmcalJet.
Double_t fTrackEfficiency
artificial tracking inefficiency (0...1)
void SetMCPt(Double_t p)
Definition: AliEmcalJet.h:269
Double_t E() const
Definition: AliEmcalJet.h:119
void SetJetAlgo(EJetAlgo_t a)
virtual ~AliEmcalJetTask()
TRandom * gRandom
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
int GlobalIndexFromLocalIndex(const U *inputObject, const int localIndex) const
void SetLabel(Int_t l)
Definition: AliEmcalJet.h:255
General jet finder task implementing a wrapper for FastJet.
static TString GenerateJetName(EJetType_t jetType, EJetAlgo_t jetAlgo, ERecoScheme_t recoScheme, Double_t radius, AliParticleContainer *partCont, AliClusterContainer *clusCont, TString tag)
UShort_t GetNumberOfConstituents() const
Definition: AliEmcalJet.h:140
Container for particles within the EMCAL framework.
Bool_t fIsPSelSet
!=true if physics selection was set
Int_t GetDefaultClusterEnergy() const
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
void AdoptClusterContainer(AliClusterContainer *cont)
TObjArray fParticleCollArray
particle/track collection array
void AddParticleConstituent(const AliVParticle *const part, Bool_t isFromEmbeddedEvent, UInt_t globalIndex)
Add new particle (track / mc particle) constituent to the given jet Note: this will append the consti...
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TPC fiducial acceptance (each eta edge narrowed by jet R)
Definition: AliEmcalJet.h:68
Bool_t fLocked
true if lock is set
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
void SetJetAcceptanceType(UInt_t type)
Definition: AliEmcalJet.h:366
Double_t fJetPhiMin
minimum phi to keep jet in output
UInt_t FindJetAcceptanceType(Double_t eta, Double_t phi, Double_t r)
int Int_t
Definition: External.C:63
Bool_t IsJetInDcal(Double_t eta, Double_t phi, Double_t r)
TPC acceptance.
Definition: AliEmcalJet.h:67
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Double_t GetJetArea(UInt_t idx) const
AliEMCALGeometry * fGeom
!emcal geometry
PHOS fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:76
void AddClusterConstituent(const AliVCluster *const clust, AliVCluster::VCluUserDefEnergy_t endef, Double_t *pvec, Bool_t isFromEmbeddedEvent, UInt_t globalIndex)
Add new cluster constituent to the given jet Note: this will append the constituent. No sorting according to particle is done.
Double_t Phi_0_2pi() const
EMCal acceptance.
Definition: AliEmcalJet.h:69
TString fJetsTag
tag of jet collection (usually = "Jets")
void SetLegacyMode(Bool_t mode)
Definition: AliFJWrapper.h:137
virtual AliVParticle * GetParticle(Int_t i=-1) const
TObjArray * fUtilities
jet utilities (gen subtractor, constituent subtractor etc.)
Double_t phiMax
void SetPtEtaPhiM(Double_t pt, Double_t eta, Double_t phi, Double_t m)
void SetAlgorithm(const fastjet::JetAlgorithm &algor)
Definition: AliFJWrapper.h:121
static const AliEmcalContainerIndexMap< TClonesArray, AliVCluster > & GetEmcalContainerIndexMap()
Get the EMCal container utils associated with particle containers.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
AliEmcalContainerIndexMap< AliClusterContainer, AliVCluster > fClusterContainerIndexMap
! Mapping between index and cluster containers
Int_t GetNClusters() const
Double_t Phi_0_2pi() const
Definition: AliEmcalJet.h:129
Double_t fJetEtaMax
maximum eta to keep jet in output
Double_t fGhostArea
ghost area
void SetMaxChargedPt(Double32_t t)
Definition: AliEmcalJet.h:263
virtual void Clear(const Option_t *="")
void SetNEF(Double_t nef)
Definition: AliEmcalJet.h:264
Bool_t GetSortedArray(Int_t indexes[], std::vector< fastjet::PseudoJet > array) const
void AddGhost(const Double_t dPx, const Double_t dPy, const Double_t dPz, const Double_t dE)
void SetNEmc(Int_t n)
Definition: AliEmcalJet.h:270
AliVCluster * GetCluster(Int_t i) const
Bool_t IsJetInDcalOnly(Double_t eta, Double_t phi, Double_t r)
virtual void Prepare(AliFJWrapper &fjw)=0
Int_t fMinMCLabel
minimum MC label value for the tracks/clusters being considered MC particles
void SetJetTask(AliEmcalJetTask *jetTask)
void SetMinJetClusE(Double_t min)
Double_t Pt() const
Definition: AliEmcalJet.h:109
void SelectCollisionCandidates(UInt_t offlineTriggerMask=AliVEvent::kMB)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
void SetGhostArea(Double_t gharea)
EJetAlgo_t fJetAlgo
jet algorithm (kt, akt, etc)
void SetNumberOfCharged(Int_t n)
Definition: AliEmcalJet.h:267
Double_t fJetPhiMax
maximum phi to keep jet in output
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
void SetAreaEta(Double_t a)
Definition: AliEmcalJet.h:257
static FJJetAlgo ConvertToFJAlgo(EJetAlgo_t algo)
void SetParticleEtaLimits(Double_t min, Double_t max)
const std::vector< fastjet::PseudoJet > & GetInputVectors() const
Definition: AliFJWrapper.h:31
void SetMinJetTrackPt(Double_t min)
Double_t fMinJetArea
min area to keep jet in output
TString fJetsName
!name of jet collection
void SetFillGhost(Bool_t b=kTRUE)
Double_t fMinJetPt
min jet pt to keep jet in output
const AliParticleIterableMomentumContainer accepted_momentum() const
void CopyMappingFrom(const AliEmcalContainerIndexMap< U2, V > &map, U *cont)
void FillJetConstituents(AliEmcalJet *jet, std::vector< fastjet::PseudoJet > &constituents, std::vector< fastjet::PseudoJet > &constituents_sub, Int_t flag=0, TString particlesSubName="")
Bool_t GetMomentum(TLorentzVector &mom, const AliVCluster *vc, Double_t mass) const
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:266
virtual void ExecOnce()
Perform steps needed to initialize the analysis.
Double_t fRadius
jet radius
Bool_t fLegacyMode
!=true to enable FJ 2.x behavior
void SetClusPtCut(Double_t cut)
Int_t fRunNumber
!run number (triggering RunChanged()
virtual void AddInputVector(Double_t px, Double_t py, Double_t pz, Double_t E, Int_t index=-99999)
void SortConstituents()
void SetMinJetClusPt(Double_t min)
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:72
bool Bool_t
Definition: External.C:53
void SetPhiRange(Double_t pmi, Double_t pma)
EDataType_t
Switch for the data type.
DCal acceptance – spans ONLY DCal (no PHOS or gap)
Definition: AliEmcalJet.h:73
void SetClusECut(Double_t cut)
void SetDefaultClusterEnergy(Int_t d)
void SetAreaPhi(Double_t a)
Definition: AliEmcalJet.h:258
virtual void InitEvent(AliFJWrapper &fjw)=0
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:123
void SetParticlePhiLimits(Double_t min, Double_t max)
void SetPtEmc(Double_t pt)
Definition: AliEmcalJet.h:271
Container structure for EMCAL clusters.
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:74
Container for MC-true particles within the EMCAL framework.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:70
static AliEmcalJetTask * AddTaskEmcalJet(const TString nTracks="usedefault", const TString nClusters="usedefault", const AliJetContainer::EJetAlgo_t jetAlgo=AliJetContainer::antikt_algorithm, const Double_t radius=0.4, const AliJetContainer::EJetType_t jetType=AliJetContainer::kFullJet, const Double_t minTrPt=0.15, const Double_t minClPt=0.30, const Double_t ghostArea=0.005, const AliJetContainer::ERecoScheme_t reco=AliJetContainer::pt_scheme, const TString tag="Jet", const Double_t minJetPt=0., const Bool_t lockTask=kTRUE, const Bool_t bFillGhosts=kFALSE)
void SetCharge(EChargeCut_t c)
Container for jet within the EMCAL jet framework.
void SetRadius(Double_t r)
AliFJWrapper fFastJetWrapper
!fastjet wrapper
void SetNumberOfClusters(Int_t n)
Definition: AliEmcalJet.h:265
virtual void Init()=0
Bool_t IsJetInPhos(Double_t eta, Double_t phi, Double_t r)
void SetAxisInEmcal(Bool_t b)
Definition: AliEmcalJet.h:261
static const AliEmcalContainerIndexMap< TClonesArray, AliVParticle > & GetEmcalContainerIndexMap()
Get the EMCal container utils associated with particle containers.
static const Int_t fgkConstIndexShift
!contituent index shift
void SetAreaEmc(Double_t a)
Definition: AliEmcalJet.h:260
void SetClusHadCorrEnergyCut(Double_t cut)
void SetNumberOfNeutrals(Int_t n)
Definition: AliEmcalJet.h:268