AliPhysics  ebc57ae (ebc57ae)
 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"
38 
39 #include "AliEmcalJetTask.h"
40 
41 using std::cout;
42 using std::endl;
43 using std::cerr;
44 
48 
50 
57  fJetsTag(),
58  fJetAlgo(AliJetContainer::antikt_algorithm),
59  fJetType(AliJetContainer::kFullJet),
60  fRecombScheme(AliJetContainer::pt_scheme),
61  fRadius(0.4),
62  fMinJetArea(0.001),
63  fMinJetPt(1.0),
64  fJetPhiMin(-10),
65  fJetPhiMax(+10),
66  fJetEtaMin(-1),
67  fJetEtaMax(+1),
68  fGhostArea(0.005),
69  fTrackEfficiency(1.),
70  fTrackEfficiencyOnlyForEmbedding(kFALSE),
71  fUtilities(0),
72  fLocked(0),
73  fJetsName(),
74  fIsInit(0),
75  fIsPSelSet(0),
76  fIsEmcPart(0),
77  fLegacyMode(kFALSE),
78  fFillGhost(kFALSE),
79  fJets(0),
80  fClusterContainerIndexMap(),
81  fParticleContainerIndexMap(),
82  fFastJetWrapper("AliEmcalJetTask","AliEmcalJetTask")
83 {
84 }
85 
92  fJetsTag("Jets"),
93  fJetAlgo(AliJetContainer::antikt_algorithm),
94  fJetType(AliJetContainer::kFullJet),
95  fRecombScheme(AliJetContainer::pt_scheme),
96  fRadius(0.4),
97  fMinJetArea(0.001),
98  fMinJetPt(1.0),
99  fJetPhiMin(-10),
100  fJetPhiMax(+10),
101  fJetEtaMin(-1),
102  fJetEtaMax(+1),
103  fGhostArea(0.005),
104  fTrackEfficiency(1.),
105  fTrackEfficiencyOnlyForEmbedding(kFALSE),
106  fUtilities(0),
107  fLocked(0),
108  fJetsName(),
109  fIsInit(0),
110  fIsPSelSet(0),
111  fIsEmcPart(0),
112  fLegacyMode(kFALSE),
113  fFillGhost(kFALSE),
114  fJets(0),
115  fClusterContainerIndexMap(),
116  fParticleContainerIndexMap(),
117  fFastJetWrapper(name,name)
118 {
119 }
120 
125 {
126 }
127 
133 {
134  if (!fUtilities) fUtilities = new TObjArray();
135  if (fUtilities->FindObject(utility)) {
136  Error("AddUtility", "Jet utility %s already connected.", utility->GetName());
137  return utility;
138  }
139  fUtilities->Add(utility);
140  utility->SetJetTask(this);
141 
142  return utility;
143 }
144 
150 {
151  TIter next(fUtilities);
152  AliEmcalJetUtility *utility = 0;
153  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Init();
154 }
160 {
161  TIter next(fUtilities);
162  AliEmcalJetUtility *utility = 0;
163  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->InitEvent(fFastJetWrapper);
164 }
165 
172 {
173  TIter next(fUtilities);
174  AliEmcalJetUtility *utility = 0;
175  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Prepare(fFastJetWrapper);
176 }
177 
183 {
184  TIter next(fUtilities);
185  AliEmcalJetUtility *utility = 0;
186  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->ProcessJet(jet, ij, fFastJetWrapper);
187 }
188 
194 {
195  TIter next(fUtilities);
196  AliEmcalJetUtility *utility = 0;
197  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Terminate(fFastJetWrapper);
198 }
199 
205 {
206  InitEvent();
207  // clear the jet array (normally a null operation)
208  fJets->Delete();
209  Int_t n = FindJets();
210 
211  if (n == 0) return kFALSE;
212 
213  FillJetBranch();
214 
215  return kTRUE;
216 }
217 
226 {
227  if (fParticleCollArray.GetEntriesFast() == 0 && fClusterCollArray.GetEntriesFast() == 0){
228  AliError("No tracks or clusters, returning.");
229  return 0;
230  }
231 
233 
234  AliDebug(2,Form("Jet type = %d", fJetType));
235 
236  Int_t iColl = 1;
237  TIter nextPartColl(&fParticleCollArray);
238  AliParticleContainer* tracks = 0;
239  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
240  AliDebug(2,Form("Tracks from collection %d: '%s'. Embedded: %i, nTracks: %i", iColl-1, tracks->GetName(), tracks->GetIsEmbedding(), tracks->GetNParticles()));
242  for (AliParticleIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
243  // artificial inefficiency
244  if (fTrackEfficiency < 1.) {
245  if (fTrackEfficiencyOnlyForEmbedding == kFALSE || (fTrackEfficiencyOnlyForEmbedding == kTRUE && tracks->GetIsEmbedding())) {
246  Double_t rnd = gRandom->Rndm();
247  if (fTrackEfficiency < rnd) {
248  AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", it.current_index()));
249  continue;
250  }
251  }
252  }
253 
254  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()));
255  Int_t uid = it.current_index() + fgkConstIndexShift * iColl;
256  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
257  }
258  iColl++;
259  }
260 
261  iColl = 1;
262  TIter nextClusColl(&fClusterCollArray);
263  AliClusterContainer* clusters = 0;
264  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
265  AliDebug(2,Form("Clusters from collection %d: '%s'. Embedded: %i, nClusters: %i", iColl-1, clusters->GetName(), clusters->GetIsEmbedding(), clusters->GetNClusters()));
267  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
268  AliDebug(2,Form("Cluster %d accepted (label = %d, energy = %.3f)", it.current_index(), it->second->GetLabel(), it->first.E()));
269  Int_t uid = -it.current_index() - fgkConstIndexShift * iColl;
270  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
271  }
272  iColl++;
273  }
274 
275  if (fFastJetWrapper.GetInputVectors().size() == 0) return 0;
276 
277  // run jet finder
279 
280  return fFastJetWrapper.GetInclusiveJets().size();
281 }
282 
289 {
291 
292  // loop over fastjet jets
293  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper.GetInclusiveJets();
294  // sort jets according to jet pt
295  static Int_t indexes[9999] = {-1};
296  GetSortedArray(indexes, jets_incl);
297 
298  AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
299  for (UInt_t ijet = 0, jetCount = 0; ijet < jets_incl.size(); ++ijet) {
300  Int_t ij = indexes[ijet];
301  AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fFastJetWrapper.GetJetArea(ij)));
302 
303  if (jets_incl[ij].perp() < fMinJetPt) continue;
304  if (fFastJetWrapper.GetJetArea(ij) < fMinJetArea) continue;
305  if ((jets_incl[ij].eta() < fJetEtaMin) || (jets_incl[ij].eta() > fJetEtaMax) ||
306  (jets_incl[ij].phi() < fJetPhiMin) || (jets_incl[ij].phi() > fJetPhiMax))
307  continue;
308 
309  AliEmcalJet *jet = new ((*fJets)[jetCount])
310  AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
311  jet->SetLabel(ij);
312 
313  fastjet::PseudoJet area(fFastJetWrapper.GetJetAreaVector(ij));
314  jet->SetArea(area.perp());
315  jet->SetAreaEta(area.eta());
316  jet->SetAreaPhi(area.phi());
317  jet->SetAreaE(area.E());
319 
320  // Fill constituent info
321  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper.GetJetConstituents(ij));
322  FillJetConstituents(jet, constituents, constituents);
323 
324  if (fGeom) {
325  if ((jet->Phi() > fGeom->GetArm1PhiMin() * TMath::DegToRad()) &&
326  (jet->Phi() < fGeom->GetArm1PhiMax() * TMath::DegToRad()) &&
327  (jet->Eta() > fGeom->GetArm1EtaMin()) &&
328  (jet->Eta() < fGeom->GetArm1EtaMax()))
329  jet->SetAxisInEmcal(kTRUE);
330  }
331 
332  ExecuteUtilities(jet, ij);
333 
334  AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), jet->GetNumberOfConstituents()));
335  jetCount++;
336  }
337 
339 }
340 
347 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
348 {
349  static Float_t pt[9999] = {0};
350 
351  const Int_t n = (Int_t)array.size();
352 
353  if (n < 1)
354  return kFALSE;
355 
356  for (Int_t i = 0; i < n; i++)
357  pt[i] = array[i].perp();
358 
359  TMath::Sort(n, pt, indexes);
360 
361  return kTRUE;
362 }
363 
370 {
371  if (fTrackEfficiency < 1.) {
372  if (gRandom) delete gRandom;
373  gRandom = new TRandom3(0);
374  }
375 
377 
378  // add jets to event if not yet there
379  if (!(InputEvent()->FindListObject(fJetsName))) {
380  fJets = new TClonesArray("AliEmcalJet");
381  fJets->SetName(fJetsName);
382  ::Info("AliEmcalJetTask::ExecOnce", "Jet collection with name '%s' has been added to the event.", fJetsName.Data());
383  InputEvent()->AddObject(fJets);
384  }
385  else {
386  AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
387  return;
388  }
389 
390  // setup fj wrapper
391  fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts);
397 
398 
399  // setting legacy mode
400  if (fLegacyMode) {
402  }
403 
404  InitUtilities();
405 
407 
408  // Setup container utils. Must be called after AliAnalysisTaskEmcal::ExecOnce() so that the
409  // containers' arrays are setup.
412 }
413 
423 void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents,
424  std::vector<fastjet::PseudoJet>& constituents_unsub, Int_t flag, TString particlesSubName)
425 {
426  Int_t nt = 0;
427  Int_t nc = 0;
428  Double_t neutralE = 0.;
429  Double_t maxCh = 0.;
430  Double_t maxNe = 0.;
431  Int_t gall = 0;
432  Int_t gemc = 0;
433  Int_t cemc = 0;
434  Int_t ncharged = 0;
435  Int_t nneutral = 0;
436  Double_t mcpt = 0.;
437  Double_t emcpt = 0.;
438  TClonesArray * particles_sub = 0;
439 
440  Int_t uid = -1;
441 
442  jet->SetNumberOfTracks(constituents.size());
443  jet->SetNumberOfClusters(constituents.size());
444 
445  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
446 
447  if (flag == 0) {
448  uid = constituents[ic].user_index();
449  }
450  else {
451  if (constituents[ic].perp()<1.e-10) {
452  uid=-1;
453  }
454  else {
455  uid = constituents[ic].user_index();
456  }
457  if (uid==0) {
458  AliError("correspondence between un/subtracted constituent not found");
459  continue;
460  }
461  }
462 
463  AliDebug(3,Form("Processing constituent %d", uid));
464  if (uid == -1) { //ghost particle
465  ++gall;
466  if (fGeom) {
467  Double_t gphi = constituents[ic].phi();
468  if (gphi < 0) gphi += TMath::TwoPi();
469  gphi *= TMath::RadToDeg();
470  Double_t geta = constituents[ic].eta();
471  if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) &&
472  (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax()))
473  ++gemc;
474  }
475 
476  if (fFillGhost) jet->AddGhost(constituents[ic].px(),
477  constituents[ic].py(),
478  constituents[ic].pz(),
479  constituents[ic].e());
480  }
481  else if (uid >= fgkConstIndexShift) { // track constituent
482  Int_t iColl = uid / fgkConstIndexShift;
483  Int_t tid = uid - iColl * fgkConstIndexShift;
484  iColl--;
485  AliDebug(3,Form("Constituent %d is a track from collection %d and with ID %d", uid, iColl, tid));
486  AliParticleContainer* partCont = GetParticleContainer(iColl);
487  if (!partCont) {
488  AliError(Form("Could not find particle container %d",iColl));
489  continue;
490  }
491  AliVParticle *t = partCont->GetParticle(tid);
492  if (!t) {
493  AliError(Form("Could not find track %d",tid));
494  continue;
495  }
496 
497  Double_t cEta = t->Eta();
498  Double_t cPhi = t->Phi();
499  Double_t cPt = t->Pt();
500  Double_t cP = t->P();
501  if (t->Charge() == 0) {
502  neutralE += cP;
503  ++nneutral;
504  if (cPt > maxNe) maxNe = cPt;
505  } else {
506  ++ncharged;
507  if (cPt > maxCh) maxCh = cPt;
508  }
509 
510  // check if MC particle
511  if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt;
512 
513  if (fGeom) {
514  if (cPhi < 0) cPhi += TMath::TwoPi();
515  cPhi *= TMath::RadToDeg();
516  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
517  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
518  emcpt += cPt;
519  ++cemc;
520  }
521  }
522 
523  if (flag == 0 || particlesSubName == "") {
525  }
526  else {
527  // Get the particle container and array corresponding to the subtracted particles
528  partCont = GetParticleContainer(particlesSubName);
529  particles_sub = partCont->GetArray();
530  // Create the new particle in the particles_sub array and add it to the jet
531  Int_t part_sub_id = particles_sub->GetEntriesFast();
532  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(dynamic_cast<AliVTrack*>(t)); // SA: probably need to be fixed!!
533  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
534  jet->AddTrackAt(fParticleContainerIndexMap.GlobalIndexFromLocalIndex(partCont, part_sub_id), nt);
535  }
536 
537  ++nt;
538  }
539  else if (uid <= -fgkConstIndexShift) { // cluster constituent
540  Int_t iColl = -uid / fgkConstIndexShift;
541  Int_t cid = -uid - iColl * fgkConstIndexShift;
542  iColl--;
543  AliDebug(3,Form("Constituent %d is a cluster from collection %d and with ID %d", uid, iColl, cid));
544  AliClusterContainer* clusCont = GetClusterContainer(iColl);
545  AliVCluster *c = clusCont->GetCluster(cid);
546 
547  if (!c) continue;
548 
550  clusCont->GetMomentum(nP, cid);
551 
552  Double_t cEta = nP.Eta();
553  Double_t cPhi = nP.Phi_0_2pi();
554  Double_t cPt = nP.Pt();
555  Double_t cP = nP.P();
556 
557  neutralE += cP;
558  if (cPt > maxNe) maxNe = cPt;
559 
560  // MC particle
561  if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
562 
563  if (fGeom) {
564  if (cPhi<0) cPhi += TMath::TwoPi();
565  cPhi *= TMath::RadToDeg();
566  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
567  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
568  emcpt += cPt;
569  ++cemc;
570  }
571  }
572 
573  if (flag == 0 || particlesSubName == "") {
575  }
576  else {
577  // Get the cluster container and array corresponding to the subtracted particles
578  clusCont = GetClusterContainer(particlesSubName);
579  particles_sub = clusCont->GetArray();
580  // Create the new particle in the particles_sub array and add it to the jet
581  Int_t part_sub_id = particles_sub->GetEntriesFast();
582  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c);
583  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
584  jet->AddClusterAt(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, part_sub_id), nc);
585  }
586 
587  ++nc;
588  ++nneutral;
589  }
590  else {
591  AliError(Form("%s: No logical way to end up here (uid = %d).", GetName(), uid));
592  continue;
593  }
594  }
595 
596  jet->SetNumberOfTracks(nt);
597  jet->SetNumberOfClusters(nc);
598  jet->SetNEF(neutralE / jet->E());
599  jet->SetMaxChargedPt(maxCh);
600  jet->SetMaxNeutralPt(maxNe);
601  if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall);
602  else jet->SetAreaEmc(-1);
603  jet->SetNEmc(cemc);
604  jet->SetNumberOfCharged(ncharged);
605  jet->SetNumberOfNeutrals(nneutral);
606  jet->SetMCPt(mcpt);
607  jet->SetPtEmc(emcpt);
608  jet->SortConstituents();
609 }
610 
618 {
619  if (fLocked) {
620  AliFatal("Jet finder task is locked! Changing properties is not allowed.");
621  return kTRUE;
622  }
623  else {
624  return kFALSE;
625  }
626 }
627 
635 {
636  if(!fIsPSelSet) {
637  fIsPSelSet = kTRUE;
638  fOfflineTriggerMask = offlineTriggerMask;
639  }
640  else {
641  AliWarning("Manually setting the event selection for jet finders is not allowed! Using trigger=old_trigger | your_trigger");
642  fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask;
643  }
644 }
645 
652 {
653  if (IsLocked()) return;
654 
655  TIter nextPartColl(&fParticleCollArray);
656  AliParticleContainer* tracks = 0;
657  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
658  tracks->SetParticleEtaLimits(emi, ema);
659  }
660 }
661 
667 {
668  if (IsLocked()) return;
669 
670  TIter nextClusColl(&fClusterCollArray);
671  AliClusterContainer* clusters = 0;
672  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
673  clusters->SetClusPtCut(min);
674  }
675 }
676 
682 {
683  if (IsLocked()) return;
684 
685  TIter nextClusColl(&fClusterCollArray);
686  AliClusterContainer* clusters = 0;
687  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
688  clusters->SetClusECut(min);
689  }
690 }
691 
697 {
698  if (IsLocked()) return;
699 
700  TIter nextPartColl(&fParticleCollArray);
701  AliParticleContainer* tracks = 0;
702  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
703  tracks->SetParticlePtCut(min);
704  }
705 }
706 
713 {
714  if (IsLocked()) return;
715 
716  TIter nextPartColl(&fParticleCollArray);
717  AliParticleContainer* tracks = 0;
718  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
719  tracks->SetParticlePhiLimits(pmi, pma);
720  }
721 }
722 
729 fastjet::JetAlgorithm AliEmcalJetTask::ConvertToFJAlgo(EJetAlgo_t algo)
730 {
731  switch(algo) {
733  return fastjet::kt_algorithm;
735  return fastjet::antikt_algorithm;
737  return fastjet::cambridge_algorithm;
739  return fastjet::genkt_algorithm;
741  return fastjet::cambridge_for_passive_algorithm;
743  return fastjet::genkt_for_passive_algorithm;
745  return fastjet::plugin_algorithm;
747  return fastjet::undefined_jet_algorithm;
748 
749  default:
750  ::Error("AliEmcalJetTask::ConvertToFJAlgo", "Jet algorithm %d not recognized!!!", algo);
751  return fastjet::undefined_jet_algorithm;
752  }
753 }
754 
761 fastjet::RecombinationScheme AliEmcalJetTask::ConvertToFJRecoScheme(ERecoScheme_t reco)
762 {
763  switch(reco) {
765  return fastjet::E_scheme;
766 
768  return fastjet::pt_scheme;
769 
771  return fastjet::pt2_scheme;
772 
774  return fastjet::Et_scheme;
775 
777  return fastjet::Et2_scheme;
778 
780  return fastjet::BIpt_scheme;
781 
783  return fastjet::BIpt2_scheme;
784 
786  return fastjet::external_scheme;
787 
788  default:
789  ::Error("AliEmcalJetTask::ConvertToFJRecoScheme", "Recombination scheme %d not recognized!!!", reco);
790  return fastjet::external_scheme;
791  }
792 }
793 
799 
800  //This method has to be called after the run number is known because it needs the EMCal geometry object.
801 
802  UInt_t jetAcceptanceType = AliEmcalJet::kUser; // all jets satify the "no acceptance cut" condition
803 
804  // Check if TPC
805  if( eta < 0.9 && eta > -0.9 ) {
806  jetAcceptanceType |= AliEmcalJet::kTPC;
807  // Check if TPCfid
808  if (eta < 0.9 - r && eta > -0.9 + r)
809  jetAcceptanceType |= AliEmcalJet::kTPCfid;
810  }
811 
812  // Check if EMCAL
813  if( IsJetInEmcal(eta, phi, 0) ) {
814  jetAcceptanceType |= AliEmcalJet::kEMCAL;
815  // Check if EMCALfid
816  if( IsJetInEmcal(eta, phi, r) )
817  jetAcceptanceType |= AliEmcalJet::kEMCALfid;
818  }
819 
820  // Check if DCAL (i.e. eta-phi rectangle spanning DCal, which includes most of PHOS)
821  if( IsJetInDcal(eta, phi, 0) ) {
822  jetAcceptanceType |= AliEmcalJet::kDCAL;
823  // Check if DCALfid
824  if( IsJetInDcal(eta, phi, r) )
825  jetAcceptanceType |= AliEmcalJet::kDCALfid;
826  }
827 
828  // Check if DCALonly (i.e. ONLY DCal, does not include any of PHOS region)
829  if( IsJetInDcalOnly(eta, phi, 0) ) {
830  jetAcceptanceType |= AliEmcalJet::kDCALonly;
831  // Check if DCALonlyfid
832  if( IsJetInDcalOnly(eta, phi, r) )
833  jetAcceptanceType |= AliEmcalJet::kDCALonlyfid;
834  }
835 
836  // Check if PHOS
837  if( IsJetInPhos(eta, phi, 0) ) {
838  jetAcceptanceType |= AliEmcalJet::kPHOS;
839  // Check if PHOSfid
840  if( IsJetInPhos(eta, phi, r) )
841  jetAcceptanceType |= AliEmcalJet::kPHOSfid;
842  }
843 
844  return jetAcceptanceType;
845 }
846 
851 {
852  if (!fGeom) return kFALSE;
853 
854  if ( eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
855  if(fRunNumber >= 177295 && fRunNumber <= 197470) {//small SM masked in 2012 and 2013
856  if ( phi < 3.135 - r && phi > 1.405 + r )
857  return kTRUE;
858  }
859  else {
860  if ( phi < fGeom->GetEMCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetArm1PhiMin() * TMath::DegToRad() + r)
861  return kTRUE;
862  }
863  }
864 
865  return kFALSE;
866 }
867 
872 {
873  if (!fGeom) return kFALSE;
874  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
875  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r)
876  return kTRUE;
877  }
878  return kFALSE;
879 }
880 
888 {
889  if (!fGeom) return kFALSE;
890 
891  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r) {
892  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r) {
893 
894  if (TMath::Abs(eta) > fGeom->GetDCALInnerExtandedEta() + r) {
895  return kTRUE;
896  }
897  if (r < 1e-6) {
898  if (phi > fGeom->GetEMCGeometry()->GetDCALStandardPhiMax() * TMath::DegToRad())
899  return kTRUE;
900  }
901 
902  }
903  }
904 
905  return kFALSE;
906 }
907 
912 {
913  Double_t etaMax = 0.130;
914  Double_t etaMin = -0.130;
915  Double_t phiMax = 320;
916  Double_t phiMin = 260; // Run 1
917  if (fRunNumber > 209121)
918  phiMin = 250; // Run 2
919 
920  if (eta < etaMax - r && eta > etaMin + r ) {
921  if (phi < phiMax * TMath::DegToRad() - r && phi > phiMin * TMath::DegToRad() + r)
922  return kTRUE;
923  }
924  return kFALSE;
925 }
926 
945  const TString nTracks, const TString nClusters,
946  const AliJetContainer::EJetAlgo_t jetAlgo, const Double_t radius, const AliJetContainer::EJetType_t jetType,
947  const Double_t minTrPt, const Double_t minClPt,
948  const Double_t ghostArea, const AliJetContainer::ERecoScheme_t reco,
949  const TString tag, const Double_t minJetPt,
950  const Bool_t lockTask, const Bool_t bFillGhosts
951 )
952 {
953  // Get the pointer to the existing analysis manager via the static access method
954  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
955  if (!mgr) {
956  ::Error("AddTaskEmcalJet", "No analysis manager to connect to.");
957  return 0;
958  }
959 
960  // Check the analysis type using the event handlers connected to the analysis manager
961  AliVEventHandler* handler = mgr->GetInputEventHandler();
962  if (!handler) {
963  ::Error("AddTaskEmcalJet", "This task requires an input event handler");
964  return 0;
965  }
966 
967  EDataType_t dataType = kUnknownDataType;
968 
969  if (handler->InheritsFrom("AliESDInputHandler")) {
970  dataType = kESD;
971  }
972  else if (handler->InheritsFrom("AliAODInputHandler")) {
973  dataType = kAOD;
974  }
975 
976  //-------------------------------------------------------
977  // Init the task and do settings
978  //-------------------------------------------------------
979 
980  TString trackName(nTracks);
981  TString clusName(nClusters);
982 
983  if (trackName == "usedefault") {
984  if (dataType == kESD) {
985  trackName = "Tracks";
986  }
987  else if (dataType == kAOD) {
988  trackName = "tracks";
989  }
990  else {
991  trackName = "";
992  }
993  }
994 
995  if (clusName == "usedefault") {
996  if (dataType == kESD) {
997  clusName = "CaloClusters";
998  }
999  else if (dataType == kAOD) {
1000  clusName = "caloClusters";
1001  }
1002  else {
1003  clusName = "";
1004  }
1005  }
1006 
1007  AliParticleContainer* partCont = 0;
1008  if (trackName == "mcparticles") {
1009  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
1010  partCont = mcpartCont;
1011  }
1012  else if (trackName == "tracks" || trackName == "Tracks") {
1013  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
1014  partCont = trackCont;
1015  }
1016  else if (!trackName.IsNull()) {
1017  partCont = new AliParticleContainer(trackName);
1018  }
1019  if (partCont) partCont->SetParticlePtCut(minTrPt);
1020 
1021  AliClusterContainer* clusCont = 0;
1022  if (!clusName.IsNull()) {
1023  clusCont = new AliClusterContainer(clusName);
1024  clusCont->SetClusECut(0.);
1025  clusCont->SetClusPtCut(0.);
1026  clusCont->SetClusHadCorrEnergyCut(minClPt);
1027  clusCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
1028  }
1029 
1030  switch (jetType) {
1032  if (partCont) partCont->SetCharge(AliParticleContainer::kCharged);
1033  break;
1035  if (partCont) partCont->SetCharge(AliParticleContainer::kNeutral);
1036  break;
1037  default:
1038  break;
1039  }
1040 
1041  TString name = AliJetContainer::GenerateJetName(jetType, jetAlgo, reco, radius, partCont, clusCont, tag);
1042 
1043  Printf("Jet task name: %s", name.Data());
1044 
1045  AliEmcalJetTask* mgrTask = static_cast<AliEmcalJetTask *>(mgr->GetTask(name.Data()));
1046  if (mgrTask) return mgrTask;
1047 
1048  AliEmcalJetTask* jetTask = new AliEmcalJetTask(name);
1049  jetTask->SetJetType(jetType);
1050  jetTask->SetJetAlgo(jetAlgo);
1051  jetTask->SetRecombScheme(reco);
1052  jetTask->SetRadius(radius);
1053  if (partCont) jetTask->AdoptParticleContainer(partCont);
1054  if (clusCont) jetTask->AdoptClusterContainer(clusCont);
1055  jetTask->SetJetsName(tag);
1056  jetTask->SetMinJetPt(minJetPt);
1057  jetTask->SetGhostArea(ghostArea);
1058 
1059  if (bFillGhosts) jetTask->SetFillGhost();
1060  if (lockTask) jetTask->SetLocked();
1061 
1062  // Final settings, pass to manager and set the containers
1063 
1064  mgr->AddTask(jetTask);
1065 
1066  // Create containers for input/output
1067  AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer();
1068  mgr->ConnectInput(jetTask, 0, cinput);
1069 
1070  TObjArray* cnt = mgr->GetContainers();
1071 
1072  return jetTask;
1073 }
Bool_t fTrackEfficiencyOnlyForEmbedding
void SetMaxNeutralPt(Double32_t t)
Definition: AliEmcalJet.h:191
void SetJetsName(const char *n)
void SetRecombScheme(ERecoScheme_t scheme)
fastjet::PseudoJet GetJetAreaVector(UInt_t idx) const
Double_t Area() const
Definition: AliEmcalJet.h:123
TObjArray fClusterCollArray
cluster collection array
void AddTrackAt(Int_t track, Int_t idx)
Definition: AliEmcalJet.h:204
void SetParticlePtCut(Double_t cut)
TClonesArray * fJets
=true ghost particles will be filled in AliEmcalJet obj
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:68
DCal acceptance – spans entire rectangular region in eta-phi (including most of PHOS) ...
Definition: AliEmcalJet.h:64
void AdoptParticleContainer(AliParticleContainer *cont)
Double_t Eta() const
Definition: AliEmcalJet.h:114
const AliClusterIterableMomentumContainer accepted_momentum() const
Base task in the EMCAL framework.
bidirectional stl iterator over the EMCAL iterable container
EJetType_t fJetType
Bool_t IsJetInEmcal(Double_t eta, Double_t phi, Double_t r)
Double_t Phi() const
Definition: AliEmcalJet.h:110
Int_t GetNParticles() const
Container with name, TClonesArray and cuts for particles.
void AddClusterAt(Int_t clus, Int_t idx)
Definition: AliEmcalJet.h:203
Bool_t fFillGhost
=true to enable FJ 2.x behavior
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:185
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:70
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:188
AliEmcalContainerIndexMap< AliParticleContainer, AliVParticle > fParticleContainerIndexMap
! Mapping between index and particle containers
ERecoScheme_t fRecombScheme
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:122
Double_t fTrackEfficiency
void SetMCPt(Double_t p)
Definition: AliEmcalJet.h:198
Double_t E() const
Definition: AliEmcalJet.h:112
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:184
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:133
Container for particles within the EMCAL framework.
Bool_t fIsPSelSet
=true if already initialized
void SetR(Double_t r)
Definition: AliFJWrapper.h:127
void AdoptClusterContainer(AliClusterContainer *cont)
TObjArray fParticleCollArray
particle/track collection array
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:61
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:34
void SetJetAcceptanceType(UInt_t type)
Definition: AliEmcalJet.h:261
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:60
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:69
Double_t Phi_0_2pi() const
EMCal acceptance.
Definition: AliEmcalJet.h:62
void SetLegacyMode(Bool_t mode)
Definition: AliFJWrapper.h:137
virtual AliVParticle * GetParticle(Int_t i=-1) const
TObjArray * fUtilities
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
contituent index shift
Int_t GetNClusters() const
Double_t Phi_0_2pi() const
Definition: AliEmcalJet.h:122
void SetMaxChargedPt(Double32_t t)
Definition: AliEmcalJet.h:192
virtual void Clear(const Option_t *="")
void SetNEF(Double_t nef)
Definition: AliEmcalJet.h:193
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:199
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:102
void SelectCollisionCandidates(UInt_t offlineTriggerMask=AliVEvent::kMB)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:125
ClassImp(AliAnalysisTaskDeltaPt) AliAnalysisTaskDeltaPt
void SetGhostArea(Double_t gharea)
EJetAlgo_t fJetAlgo
void SetNumberOfCharged(Int_t n)
Definition: AliEmcalJet.h:196
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
void SetAreaEta(Double_t a)
Definition: AliEmcalJet.h:186
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)
void SetFillGhost(Bool_t b=kTRUE)
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:44
void SetNumberOfTracks(Int_t n)
Definition: AliEmcalJet.h:195
virtual void ExecOnce()
Perform steps needed to initialize the analysis.
Bool_t fLegacyMode
=true if emcal particles are given as input (for clusters)
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:65
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:66
void SetClusECut(Double_t cut)
void SetDefaultClusterEnergy(Int_t d)
void SetAreaPhi(Double_t a)
Definition: AliEmcalJet.h:187
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:200
Container structure for EMCAL clusters.
DCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:67
Container for MC-true particles within the EMCAL framework.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:63
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
jet collection
void SetNumberOfClusters(Int_t n)
Definition: AliEmcalJet.h:194
virtual void Init()=0
Bool_t IsJetInPhos(Double_t eta, Double_t phi, Double_t r)
void SetAxisInEmcal(Bool_t b)
Definition: AliEmcalJet.h:190
static const AliEmcalContainerIndexMap< TClonesArray, AliVParticle > & GetEmcalContainerIndexMap()
Get the EMCal container utils associated with particle containers.
static const Int_t fgkConstIndexShift
fastjet wrapper
void SetAreaEmc(Double_t a)
Definition: AliEmcalJet.h:189
void SetClusHadCorrEnergyCut(Double_t cut)
void SetNumberOfNeutrals(Int_t n)
Definition: AliEmcalJet.h:197