AliPhysics  d497afb (d497afb)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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),
109  fTrackEfficiencyOnlyForEmbedding(kFALSE),
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),
120  fClusterContainerIndexMap(),
121  fParticleContainerIndexMap()
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 
382  // add jets to event if not yet there
383  if (!(InputEvent()->FindListObject(fJetsName))) {
384  fJets = new TClonesArray("AliEmcalJet");
385  fJets->SetName(fJetsName);
386  ::Info("AliEmcalJetTask::ExecOnce", "Jet collection with name '%s' has been added to the event.", fJetsName.Data());
387  InputEvent()->AddObject(fJets);
388  }
389  else {
390  AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
391  return;
392  }
393 
394  // setup fj wrapper
395  fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts);
401 
402 
403  // setting legacy mode
404  if (fLegacyMode) {
406  }
407 
408  InitUtilities();
409 
411 
412  // Setup container utils. Must be called after AliAnalysisTaskEmcal::ExecOnce() so that the
413  // containers' arrays are setup.
416 }
417 
427 void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents,
428  std::vector<fastjet::PseudoJet>& constituents_unsub, Int_t flag, TString particlesSubName)
429 {
430  Int_t nt = 0;
431  Int_t nc = 0;
432  Double_t neutralE = 0.;
433  Double_t maxCh = 0.;
434  Double_t maxNe = 0.;
435  Int_t gall = 0;
436  Int_t gemc = 0;
437  Int_t cemc = 0;
438  Int_t ncharged = 0;
439  Int_t nneutral = 0;
440  Double_t mcpt = 0.;
441  Double_t emcpt = 0.;
442  TClonesArray * particles_sub = 0;
443 
444  Int_t uid = -1;
445 
446  jet->SetNumberOfTracks(constituents.size());
447  jet->SetNumberOfClusters(constituents.size());
448 
449  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
450 
451  if (flag == 0) {
452  uid = constituents[ic].user_index();
453  }
454  else {
455  if (constituents[ic].perp()<1.e-10) {
456  uid=-1;
457  }
458  else {
459  uid = constituents[ic].user_index();
460  }
461  if (uid==0) {
462  AliError("correspondence between un/subtracted constituent not found");
463  continue;
464  }
465  }
466 
467  AliDebug(3,Form("Processing constituent %d", uid));
468  if (uid == -1) { //ghost particle
469  ++gall;
470  if (fGeom) {
471  Double_t gphi = constituents[ic].phi();
472  if (gphi < 0) gphi += TMath::TwoPi();
473  gphi *= TMath::RadToDeg();
474  Double_t geta = constituents[ic].eta();
475  if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) &&
476  (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax()))
477  ++gemc;
478  }
479 
480  if (fFillGhost) jet->AddGhost(constituents[ic].px(),
481  constituents[ic].py(),
482  constituents[ic].pz(),
483  constituents[ic].e());
484  }
485  else if (uid >= fgkConstIndexShift) { // track constituent
486  Int_t iColl = uid / fgkConstIndexShift;
487  Int_t tid = uid - iColl * fgkConstIndexShift;
488  iColl--;
489  AliDebug(3,Form("Constituent %d is a track from collection %d and with ID %d", uid, iColl, tid));
490  AliParticleContainer* partCont = GetParticleContainer(iColl);
491  if (!partCont) {
492  AliError(Form("Could not find particle container %d",iColl));
493  continue;
494  }
495  AliVParticle *t = partCont->GetParticle(tid);
496  if (!t) {
497  AliError(Form("Could not find track %d",tid));
498  continue;
499  }
500  if(fFillConstituents){
501  jet->AddParticleConstituent(t, partCont->GetIsEmbedding(), fParticleContainerIndexMap.GlobalIndexFromLocalIndex(partCont, tid));
502  }
503 
504  Double_t cEta = t->Eta();
505  Double_t cPhi = t->Phi();
506  Double_t cPt = t->Pt();
507  Double_t cP = t->P();
508  if (t->Charge() == 0) {
509  neutralE += cP;
510  ++nneutral;
511  if (cPt > maxNe) maxNe = cPt;
512  } else {
513  ++ncharged;
514  if (cPt > maxCh) maxCh = cPt;
515  }
516 
517  // check if MC particle
518  if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt;
519 
520  if (fGeom) {
521  if (cPhi < 0) cPhi += TMath::TwoPi();
522  cPhi *= TMath::RadToDeg();
523  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
524  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
525  emcpt += cPt;
526  ++cemc;
527  }
528  }
529 
530  if (flag == 0 || particlesSubName == "") {
532  }
533  else {
534  // Get the particle container and array corresponding to the subtracted particles
535  partCont = GetParticleContainer(particlesSubName);
536  particles_sub = partCont->GetArray();
537  // Create the new particle in the particles_sub array and add it to the jet
538  Int_t part_sub_id = particles_sub->GetEntriesFast();
539  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(dynamic_cast<AliVTrack*>(t)); // SA: probably need to be fixed!!
540  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
541  jet->AddTrackAt(fParticleContainerIndexMap.GlobalIndexFromLocalIndex(partCont, part_sub_id), nt);
542  }
543 
544  ++nt;
545  }
546  else if (uid <= -fgkConstIndexShift) { // cluster constituent
547  Int_t iColl = -uid / fgkConstIndexShift;
548  Int_t cid = -uid - iColl * fgkConstIndexShift;
549  iColl--;
550  AliDebug(3,Form("Constituent %d is a cluster from collection %d and with ID %d", uid, iColl, cid));
551  AliClusterContainer* clusCont = GetClusterContainer(iColl);
552  AliVCluster *c = clusCont->GetCluster(cid);
553 
554  if (!c) continue;
555 
557  clusCont->GetMomentum(nP, cid);
558 
559  Double_t cEta = nP.Eta();
560  Double_t cPhi = nP.Phi_0_2pi();
561  Double_t cPt = nP.Pt();
562  Double_t cP = nP.P();
563  Double_t pvec[3] = {nP.Px(), nP.Py(), nP.Pz()};
564  if(fFillConstituents) jet->AddClusterConstituent(c, (AliVCluster::VCluUserDefEnergy_t)clusCont->GetDefaultClusterEnergy(), pvec, clusCont->GetIsEmbedding(), fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, cid));
565 
566  neutralE += cP;
567  if (cPt > maxNe) maxNe = cPt;
568 
569  // MC particle
570  if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
571 
572  if (fGeom) {
573  if (cPhi<0) cPhi += TMath::TwoPi();
574  cPhi *= TMath::RadToDeg();
575  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
576  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
577  emcpt += cPt;
578  ++cemc;
579  }
580  }
581 
582  if (flag == 0 || particlesSubName == "") {
584  }
585  else {
586  // Get the cluster container and array corresponding to the subtracted particles
587  clusCont = GetClusterContainer(particlesSubName);
588  particles_sub = clusCont->GetArray();
589  // Create the new particle in the particles_sub array and add it to the jet
590  Int_t part_sub_id = particles_sub->GetEntriesFast();
591  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c);
592  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
593  jet->AddClusterAt(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, part_sub_id), nc);
594  }
595 
596  ++nc;
597  ++nneutral;
598  }
599  else {
600  AliError(Form("%s: No logical way to end up here (uid = %d).", GetName(), uid));
601  continue;
602  }
603  }
604 
605  jet->SetNumberOfTracks(nt);
606  jet->SetNumberOfClusters(nc);
607  jet->SetNEF(neutralE / jet->E());
608  jet->SetMaxChargedPt(maxCh);
609  jet->SetMaxNeutralPt(maxNe);
610  if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall);
611  else jet->SetAreaEmc(-1);
612  jet->SetNEmc(cemc);
613  jet->SetNumberOfCharged(ncharged);
614  jet->SetNumberOfNeutrals(nneutral);
615  jet->SetMCPt(mcpt);
616  jet->SetPtEmc(emcpt);
617  jet->SortConstituents();
618 }
619 
627 {
628  if (fLocked) {
629  AliFatal("Jet finder task is locked! Changing properties is not allowed.");
630  return kTRUE;
631  }
632  else {
633  return kFALSE;
634  }
635 }
636 
644 {
645  if(!fIsPSelSet) {
646  fIsPSelSet = kTRUE;
647  fOfflineTriggerMask = offlineTriggerMask;
648  }
649  else {
650  AliWarning("Manually setting the event selection for jet finders is not allowed! Using trigger=old_trigger | your_trigger");
651  fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask;
652  }
653 }
654 
661 {
662  if (IsLocked()) return;
663 
664  TIter nextPartColl(&fParticleCollArray);
665  AliParticleContainer* tracks = 0;
666  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
667  tracks->SetParticleEtaLimits(emi, ema);
668  }
669 }
670 
676 {
677  if (IsLocked()) return;
678 
679  TIter nextClusColl(&fClusterCollArray);
680  AliClusterContainer* clusters = 0;
681  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
682  clusters->SetClusPtCut(min);
683  }
684 }
685 
691 {
692  if (IsLocked()) return;
693 
694  TIter nextClusColl(&fClusterCollArray);
695  AliClusterContainer* clusters = 0;
696  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
697  clusters->SetClusECut(min);
698  }
699 }
700 
706 {
707  if (IsLocked()) return;
708 
709  TIter nextPartColl(&fParticleCollArray);
710  AliParticleContainer* tracks = 0;
711  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
712  tracks->SetParticlePtCut(min);
713  }
714 }
715 
722 {
723  if (IsLocked()) return;
724 
725  TIter nextPartColl(&fParticleCollArray);
726  AliParticleContainer* tracks = 0;
727  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
728  tracks->SetParticlePhiLimits(pmi, pma);
729  }
730 }
731 
738 fastjet::JetAlgorithm AliEmcalJetTask::ConvertToFJAlgo(EJetAlgo_t algo)
739 {
740  switch(algo) {
742  return fastjet::kt_algorithm;
744  return fastjet::antikt_algorithm;
746  return fastjet::cambridge_algorithm;
748  return fastjet::genkt_algorithm;
750  return fastjet::cambridge_for_passive_algorithm;
752  return fastjet::genkt_for_passive_algorithm;
754  return fastjet::plugin_algorithm;
756  return fastjet::undefined_jet_algorithm;
757 
758  default:
759  ::Error("AliEmcalJetTask::ConvertToFJAlgo", "Jet algorithm %d not recognized!!!", algo);
760  return fastjet::undefined_jet_algorithm;
761  }
762 }
763 
770 fastjet::RecombinationScheme AliEmcalJetTask::ConvertToFJRecoScheme(ERecoScheme_t reco)
771 {
772  switch(reco) {
774  return fastjet::E_scheme;
775 
777  return fastjet::pt_scheme;
778 
780  return fastjet::pt2_scheme;
781 
783  return fastjet::Et_scheme;
784 
786  return fastjet::Et2_scheme;
787 
789  return fastjet::BIpt_scheme;
790 
792  return fastjet::BIpt2_scheme;
793 
795  return fastjet::external_scheme;
796 
797  default:
798  ::Error("AliEmcalJetTask::ConvertToFJRecoScheme", "Recombination scheme %d not recognized!!!", reco);
799  return fastjet::external_scheme;
800  }
801 }
802 
808 
809  //This method has to be called after the run number is known because it needs the EMCal geometry object.
810 
811  UInt_t jetAcceptanceType = AliEmcalJet::kUser; // all jets satify the "no acceptance cut" condition
812 
813  // Check if TPC
814  if( eta < 0.9 && eta > -0.9 ) {
815  jetAcceptanceType |= AliEmcalJet::kTPC;
816  // Check if TPCfid
817  if (eta < 0.9 - r && eta > -0.9 + r)
818  jetAcceptanceType |= AliEmcalJet::kTPCfid;
819  }
820 
821  // Check if EMCAL
822  if( IsJetInEmcal(eta, phi, 0) ) {
823  jetAcceptanceType |= AliEmcalJet::kEMCAL;
824  // Check if EMCALfid
825  if( IsJetInEmcal(eta, phi, r) )
826  jetAcceptanceType |= AliEmcalJet::kEMCALfid;
827  }
828 
829  // Check if DCAL (i.e. eta-phi rectangle spanning DCal, which includes most of PHOS)
830  if( IsJetInDcal(eta, phi, 0) ) {
831  jetAcceptanceType |= AliEmcalJet::kDCAL;
832  // Check if DCALfid
833  if( IsJetInDcal(eta, phi, r) )
834  jetAcceptanceType |= AliEmcalJet::kDCALfid;
835  }
836 
837  // Check if DCALonly (i.e. ONLY DCal, does not include any of PHOS region)
838  if( IsJetInDcalOnly(eta, phi, 0) ) {
839  jetAcceptanceType |= AliEmcalJet::kDCALonly;
840  // Check if DCALonlyfid
841  if( IsJetInDcalOnly(eta, phi, r) )
842  jetAcceptanceType |= AliEmcalJet::kDCALonlyfid;
843  }
844 
845  // Check if PHOS
846  if( IsJetInPhos(eta, phi, 0) ) {
847  jetAcceptanceType |= AliEmcalJet::kPHOS;
848  // Check if PHOSfid
849  if( IsJetInPhos(eta, phi, r) )
850  jetAcceptanceType |= AliEmcalJet::kPHOSfid;
851  }
852 
853  return jetAcceptanceType;
854 }
855 
860 {
861  if (!fGeom) return kFALSE;
862 
863  if ( eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
864  if(fRunNumber >= 177295 && fRunNumber <= 197470) {//small SM masked in 2012 and 2013
865  if ( phi < 3.135 - r && phi > 1.405 + r )
866  return kTRUE;
867  }
868  else {
869  if ( phi < fGeom->GetEMCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetArm1PhiMin() * TMath::DegToRad() + r)
870  return kTRUE;
871  }
872  }
873 
874  return kFALSE;
875 }
876 
881 {
882  if (!fGeom) return kFALSE;
883  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
884  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r)
885  return kTRUE;
886  }
887  return kFALSE;
888 }
889 
897 {
898  if (!fGeom) return kFALSE;
899 
900  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r) {
901  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r) {
902 
903  if (TMath::Abs(eta) > fGeom->GetDCALInnerExtandedEta() + r) {
904  return kTRUE;
905  }
906  if (r < 1e-6) {
907  if (phi > fGeom->GetEMCGeometry()->GetDCALStandardPhiMax() * TMath::DegToRad())
908  return kTRUE;
909  }
910 
911  }
912  }
913 
914  return kFALSE;
915 }
916 
921 {
922  Double_t etaMax = 0.130;
923  Double_t etaMin = -0.130;
924  Double_t phiMax = 320;
925  Double_t phiMin = 260; // Run 1
926  if (fRunNumber > 209121)
927  phiMin = 250; // Run 2
928 
929  if (eta < etaMax - r && eta > etaMin + r ) {
930  if (phi < phiMax * TMath::DegToRad() - r && phi > phiMin * TMath::DegToRad() + r)
931  return kTRUE;
932  }
933  return kFALSE;
934 }
935 
954  const TString nTracks, const TString nClusters,
955  const AliJetContainer::EJetAlgo_t jetAlgo, const Double_t radius, const AliJetContainer::EJetType_t jetType,
956  const Double_t minTrPt, const Double_t minClPt,
957  const Double_t ghostArea, const AliJetContainer::ERecoScheme_t reco,
958  const TString tag, const Double_t minJetPt,
959  const Bool_t lockTask, const Bool_t bFillGhosts
960 )
961 {
962  // Get the pointer to the existing analysis manager via the static access method
963  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
964  if (!mgr) {
965  ::Error("AddTaskEmcalJet", "No analysis manager to connect to.");
966  return 0;
967  }
968 
969  // Check the analysis type using the event handlers connected to the analysis manager
970  AliVEventHandler* handler = mgr->GetInputEventHandler();
971  if (!handler) {
972  ::Error("AddTaskEmcalJet", "This task requires an input event handler");
973  return 0;
974  }
975 
976  EDataType_t dataType = kUnknownDataType;
977 
978  if (handler->InheritsFrom("AliESDInputHandler")) {
979  dataType = kESD;
980  }
981  else if (handler->InheritsFrom("AliAODInputHandler")) {
982  dataType = kAOD;
983  }
984 
985  //-------------------------------------------------------
986  // Init the task and do settings
987  //-------------------------------------------------------
988 
989  TString trackName(nTracks);
990  TString clusName(nClusters);
991 
992  if (trackName == "usedefault") {
993  if (dataType == kESD) {
994  trackName = "Tracks";
995  }
996  else if (dataType == kAOD) {
997  trackName = "tracks";
998  }
999  else {
1000  trackName = "";
1001  }
1002  }
1003 
1004  if (clusName == "usedefault") {
1005  if (dataType == kESD) {
1006  clusName = "CaloClusters";
1007  }
1008  else if (dataType == kAOD) {
1009  clusName = "caloClusters";
1010  }
1011  else {
1012  clusName = "";
1013  }
1014  }
1015 
1016  AliParticleContainer* partCont = 0;
1017  if (trackName == "mcparticles") {
1018  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
1019  partCont = mcpartCont;
1020  }
1021  else if (trackName == "tracks" || trackName == "Tracks") {
1022  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
1023  partCont = trackCont;
1024  }
1025  else if (!trackName.IsNull()) {
1026  partCont = new AliParticleContainer(trackName);
1027  }
1028  if (partCont) partCont->SetParticlePtCut(minTrPt);
1029 
1030  AliClusterContainer* clusCont = 0;
1031  if (!clusName.IsNull()) {
1032  clusCont = new AliClusterContainer(clusName);
1033  clusCont->SetClusECut(0.);
1034  clusCont->SetClusPtCut(0.);
1035  clusCont->SetClusHadCorrEnergyCut(minClPt);
1036  clusCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
1037  }
1038 
1039  switch (jetType) {
1041  if (partCont) partCont->SetCharge(AliParticleContainer::kCharged);
1042  break;
1044  if (partCont) partCont->SetCharge(AliParticleContainer::kNeutral);
1045  break;
1046  default:
1047  break;
1048  }
1049 
1050  TString name = AliJetContainer::GenerateJetName(jetType, jetAlgo, reco, radius, partCont, clusCont, tag);
1051 
1052  Printf("Jet task name: %s", name.Data());
1053 
1054  AliEmcalJetTask* mgrTask = static_cast<AliEmcalJetTask *>(mgr->GetTask(name.Data()));
1055  if (mgrTask) return mgrTask;
1056 
1057  AliEmcalJetTask* jetTask = new AliEmcalJetTask(name);
1058  jetTask->SetJetType(jetType);
1059  jetTask->SetJetAlgo(jetAlgo);
1060  jetTask->SetRecombScheme(reco);
1061  jetTask->SetRadius(radius);
1062  if (partCont) jetTask->AdoptParticleContainer(partCont);
1063  if (clusCont) jetTask->AdoptClusterContainer(clusCont);
1064  jetTask->SetJetsName(tag);
1065  jetTask->SetMinJetPt(minJetPt);
1066  jetTask->SetGhostArea(ghostArea);
1067 
1068  if (bFillGhosts) jetTask->SetFillGhost();
1069  if (lockTask) jetTask->SetLocked();
1070 
1071  // Final settings, pass to manager and set the containers
1072 
1073  mgr->AddTask(jetTask);
1074 
1075  // Create containers for input/output
1076  AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer();
1077  mgr->ConnectInput(jetTask, 0, cinput);
1078 
1079  return jetTask;
1080 }
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
void SetEtaRange(Double_t emi, Double_t ema)
Bool_t IsLocked() const
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
Bool_t fFillGhost
=true ghost particles will be filled in AliEmcalJet obj
AliEmcalJetUtility * AddUtility(AliEmcalJetUtility *utility)
Declaration of class AliTLorentzVector.
virtual void ProcessJet(AliEmcalJet *jet, Int_t ij, AliFJWrapper &fjw)=0
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.)
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