AliPhysics  58f3d52 (58f3d52)
 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  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
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