AliPhysics  cdeda5a (cdeda5a)
 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 }
155 
162 {
163  TIter next(fUtilities);
164  AliEmcalJetUtility *utility = 0;
165  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Prepare(fFastJetWrapper);
166 }
167 
173 {
174  TIter next(fUtilities);
175  AliEmcalJetUtility *utility = 0;
176  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->ProcessJet(jet, ij, fFastJetWrapper);
177 }
178 
184 {
185  TIter next(fUtilities);
186  AliEmcalJetUtility *utility = 0;
187  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Terminate(fFastJetWrapper);
188 }
189 
195 {
196  // clear the jet array (normally a null operation)
197  fJets->Delete();
198 
199  Int_t n = FindJets();
200 
201  if (n == 0) return kFALSE;
202 
203  FillJetBranch();
204 
205  return kTRUE;
206 }
207 
216 {
217  if (fParticleCollArray.GetEntriesFast() == 0 && fClusterCollArray.GetEntriesFast() == 0){
218  AliError("No tracks or clusters, returning.");
219  return 0;
220  }
221 
223 
224  AliDebug(2,Form("Jet type = %d", fJetType));
225 
226  Int_t iColl = 1;
227  TIter nextPartColl(&fParticleCollArray);
228  AliParticleContainer* tracks = 0;
229  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
230  AliDebug(2,Form("Tracks from collection %d: '%s'. Embedded: %i, nTracks: %i", iColl-1, tracks->GetName(), tracks->GetIsEmbedding(), tracks->GetNParticles()));
232  for (AliParticleIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
233  // artificial inefficiency
234  if (fTrackEfficiency < 1.) {
235  if (fTrackEfficiencyOnlyForEmbedding == kFALSE || (fTrackEfficiencyOnlyForEmbedding == kTRUE && tracks->GetIsEmbedding())) {
236  Double_t rnd = gRandom->Rndm();
237  if (fTrackEfficiency < rnd) {
238  AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", it.current_index()));
239  continue;
240  }
241  }
242  }
243 
244  AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", it.current_index(), it->second->GetLabel(), it->first.Pt()));
245  Int_t uid = it.current_index() + fgkConstIndexShift * iColl;
246  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
247  }
248  iColl++;
249  }
250 
251  iColl = 1;
252  TIter nextClusColl(&fClusterCollArray);
253  AliClusterContainer* clusters = 0;
254  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
255  AliDebug(2,Form("Clusters from collection %d: '%s'. Embedded: %i, nClusters: %i", iColl-1, clusters->GetName(), clusters->GetIsEmbedding(), clusters->GetNClusters()));
257  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
258  AliDebug(2,Form("Cluster %d accepted (label = %d, energy = %.3f)", it.current_index(), it->second->GetLabel(), it->first.E()));
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  if (fFastJetWrapper.GetInputVectors().size() == 0) return 0;
266 
267  // run jet finder
269 
270  return fFastJetWrapper.GetInclusiveJets().size();
271 }
272 
279 {
281 
282  // loop over fastjet jets
283  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper.GetInclusiveJets();
284  // sort jets according to jet pt
285  static Int_t indexes[9999] = {-1};
286  GetSortedArray(indexes, jets_incl);
287 
288  AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
289  for (UInt_t ijet = 0, jetCount = 0; ijet < jets_incl.size(); ++ijet) {
290  Int_t ij = indexes[ijet];
291  AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fFastJetWrapper.GetJetArea(ij)));
292 
293  if (jets_incl[ij].perp() < fMinJetPt) continue;
294  if (fFastJetWrapper.GetJetArea(ij) < fMinJetArea) continue;
295  if ((jets_incl[ij].eta() < fJetEtaMin) || (jets_incl[ij].eta() > fJetEtaMax) ||
296  (jets_incl[ij].phi() < fJetPhiMin) || (jets_incl[ij].phi() > fJetPhiMax))
297  continue;
298 
299  AliEmcalJet *jet = new ((*fJets)[jetCount])
300  AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
301  jet->SetLabel(ij);
302 
303  fastjet::PseudoJet area(fFastJetWrapper.GetJetAreaVector(ij));
304  jet->SetArea(area.perp());
305  jet->SetAreaEta(area.eta());
306  jet->SetAreaPhi(area.phi());
307  jet->SetAreaE(area.E());
309 
310  // Fill constituent info
311  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper.GetJetConstituents(ij));
312  FillJetConstituents(jet, constituents, constituents);
313 
314  if (fGeom) {
315  if ((jet->Phi() > fGeom->GetArm1PhiMin() * TMath::DegToRad()) &&
316  (jet->Phi() < fGeom->GetArm1PhiMax() * TMath::DegToRad()) &&
317  (jet->Eta() > fGeom->GetArm1EtaMin()) &&
318  (jet->Eta() < fGeom->GetArm1EtaMax()))
319  jet->SetAxisInEmcal(kTRUE);
320  }
321 
322  ExecuteUtilities(jet, ij);
323 
324  AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), jet->GetNumberOfConstituents()));
325  jetCount++;
326  }
327 
329 }
330 
337 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
338 {
339  static Float_t pt[9999] = {0};
340 
341  const Int_t n = (Int_t)array.size();
342 
343  if (n < 1)
344  return kFALSE;
345 
346  for (Int_t i = 0; i < n; i++)
347  pt[i] = array[i].perp();
348 
349  TMath::Sort(n, pt, indexes);
350 
351  return kTRUE;
352 }
353 
360 {
361  if (fTrackEfficiency < 1.) {
362  if (gRandom) delete gRandom;
363  gRandom = new TRandom3(0);
364  }
365 
367 
368  // add jets to event if not yet there
369  if (!(InputEvent()->FindListObject(fJetsName))) {
370  fJets = new TClonesArray("AliEmcalJet");
371  fJets->SetName(fJetsName);
372  ::Info("AliEmcalJetTask::ExecOnce", "Jet collection with name '%s' has been added to the event.", fJetsName.Data());
373  InputEvent()->AddObject(fJets);
374  }
375  else {
376  AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
377  return;
378  }
379 
380  // setup fj wrapper
381  fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts);
387 
388  // setting legacy mode
389  if (fLegacyMode) {
391  }
392 
393  InitUtilities();
394 
396 
397  // Setup container utils. Must be called after AliAnalysisTaskEmcal::ExecOnce() so that the
398  // containers' arrays are setup.
401 }
402 
412 void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents,
413  std::vector<fastjet::PseudoJet>& constituents_unsub, Int_t flag, TClonesArray *particles_sub)
414 {
415  Int_t nt = 0;
416  Int_t nc = 0;
417  Double_t neutralE = 0.;
418  Double_t maxCh = 0.;
419  Double_t maxNe = 0.;
420  Int_t gall = 0;
421  Int_t gemc = 0;
422  Int_t cemc = 0;
423  Int_t ncharged = 0;
424  Int_t nneutral = 0;
425  Double_t mcpt = 0.;
426  Double_t emcpt = 0.;
427 
428  Int_t uid = -1;
429 
430  jet->SetNumberOfTracks(constituents.size());
431  jet->SetNumberOfClusters(constituents.size());
432 
433  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
434 
435  if (flag == 0) {
436  uid = constituents[ic].user_index();
437  }
438  else {
439  if (constituents[ic].perp()<1.e-10) {
440  uid=-1;
441  }
442  else {
443  uid = GetIndexSub(constituents[ic].phi(), constituents_unsub);
444  }
445  if (uid==0) {
446  AliError("correspondence between un/subtracted constituent not found");
447  continue;
448  }
449  }
450 
451  AliDebug(3,Form("Processing constituent %d", uid));
452  if (uid == -1) { //ghost particle
453  ++gall;
454  if (fGeom) {
455  Double_t gphi = constituents[ic].phi();
456  if (gphi < 0) gphi += TMath::TwoPi();
457  gphi *= TMath::RadToDeg();
458  Double_t geta = constituents[ic].eta();
459  if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) &&
460  (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax()))
461  ++gemc;
462  }
463 
464  if (fFillGhost) jet->AddGhost(constituents[ic].px(),
465  constituents[ic].py(),
466  constituents[ic].pz(),
467  constituents[ic].e());
468  }
469  else if (uid >= fgkConstIndexShift) { // track constituent
470  Int_t iColl = uid / fgkConstIndexShift;
471  Int_t tid = uid - iColl * fgkConstIndexShift;
472  iColl--;
473  AliDebug(3,Form("Constituent %d is a track from collection %d and with ID %d", uid, iColl, tid));
474  AliParticleContainer* partCont = GetParticleContainer(iColl);
475  if (!partCont) {
476  AliError(Form("Could not find particle container %d",iColl));
477  continue;
478  }
479  AliVParticle *t = partCont->GetParticle(tid);
480  if (!t) {
481  AliError(Form("Could not find track %d",tid));
482  continue;
483  }
484 
485  Double_t cEta = t->Eta();
486  Double_t cPhi = t->Phi();
487  Double_t cPt = t->Pt();
488  Double_t cP = t->P();
489  if (t->Charge() == 0) {
490  neutralE += cP;
491  ++nneutral;
492  if (cPt > maxNe) maxNe = cPt;
493  } else {
494  ++ncharged;
495  if (cPt > maxCh) maxCh = cPt;
496  }
497 
498  // check if MC particle
499  if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt;
500 
501  if (fGeom) {
502  if (cPhi < 0) cPhi += TMath::TwoPi();
503  cPhi *= TMath::RadToDeg();
504  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
505  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
506  emcpt += cPt;
507  ++cemc;
508  }
509  }
510 
511  if (flag == 0 || particles_sub == 0) {
513  }
514  else {
515  Int_t part_sub_id = particles_sub->GetEntriesFast();
516  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(dynamic_cast<AliVTrack*>(t)); // SA: probably need to be fixed!!
517  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
518  jet->AddTrackAt(part_sub_id, nt);
519  }
520 
521  ++nt;
522  }
523  else if (uid <= -fgkConstIndexShift) { // cluster constituent
524  Int_t iColl = -uid / fgkConstIndexShift;
525  Int_t cid = -uid - iColl * fgkConstIndexShift;
526  iColl--;
527  AliDebug(3,Form("Constituent %d is a cluster from collection %d and with ID %d", uid, iColl, cid));
528  AliClusterContainer* clusCont = GetClusterContainer(iColl);
529  AliVCluster *c = clusCont->GetCluster(cid);
530 
531  if (!c) continue;
532 
534  clusCont->GetMomentum(nP, cid);
535 
536  Double_t cEta = nP.Eta();
537  Double_t cPhi = nP.Phi_0_2pi();
538  Double_t cPt = nP.Pt();
539  Double_t cP = nP.P();
540 
541  neutralE += cP;
542  if (cPt > maxNe) maxNe = cPt;
543 
544  // MC particle
545  if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
546 
547  if (fGeom) {
548  if (cPhi<0) cPhi += TMath::TwoPi();
549  cPhi *= TMath::RadToDeg();
550  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
551  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
552  emcpt += cPt;
553  ++cemc;
554  }
555  }
556 
557  if (flag == 0 || particles_sub == 0) {
559  }
560  else {
561  Int_t part_sub_id = particles_sub->GetEntriesFast();
562  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c);
563  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
564  jet->AddTrackAt(part_sub_id, nt);
565  }
566 
567  ++nc;
568  ++nneutral;
569  }
570  else {
571  AliError(Form("%s: No logical way to end up here (uid = %d).", GetName(), uid));
572  continue;
573  }
574  }
575 
576  jet->SetNumberOfTracks(nt);
577  jet->SetNumberOfClusters(nc);
578  jet->SetNEF(neutralE / jet->E());
579  jet->SetMaxChargedPt(maxCh);
580  jet->SetMaxNeutralPt(maxNe);
581  if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall);
582  else jet->SetAreaEmc(-1);
583  jet->SetNEmc(cemc);
584  jet->SetNumberOfCharged(ncharged);
585  jet->SetNumberOfNeutrals(nneutral);
586  jet->SetMCPt(mcpt);
587  jet->SetPtEmc(emcpt);
588  jet->SortConstituents();
589 }
590 
597 Int_t AliEmcalJetTask::GetIndexSub(Double_t phi_sub, std::vector<fastjet::PseudoJet>& constituents_unsub)
598 {
599  Double_t dphi=0;
600  Double_t phimin=100;
601  Int_t index=0;
602  for (UInt_t ii = 0; ii < constituents_unsub.size(); ii++) {
603  dphi=0;
604  Double_t phi_unsub = constituents_unsub[ii].phi();
605  dphi=phi_unsub-phi_sub;
606  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
607  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
608  if(TMath::Abs(dphi)<phimin){ phimin=TMath::Abs(dphi);
609  index=ii;} }
610  if (constituents_unsub[index].user_index()!=-1) return constituents_unsub[index].user_index();
611 
612  return 0;
613 }
614 
622 {
623  if (fLocked) {
624  AliFatal("Jet finder task is locked! Changing properties is not allowed.");
625  return kTRUE;
626  }
627  else {
628  return kFALSE;
629  }
630 }
631 
639 {
640  if(!fIsPSelSet) {
641  fIsPSelSet = kTRUE;
642  fOfflineTriggerMask = offlineTriggerMask;
643  }
644  else {
645  AliWarning("Manually setting the event selection for jet finders is not allowed! Using trigger=old_trigger | your_trigger");
646  fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask;
647  }
648 }
649 
656 {
657  if (IsLocked()) return;
658 
659  TIter nextPartColl(&fParticleCollArray);
660  AliParticleContainer* tracks = 0;
661  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
662  tracks->SetParticleEtaLimits(emi, ema);
663  }
664 }
665 
671 {
672  if (IsLocked()) return;
673 
674  TIter nextClusColl(&fClusterCollArray);
675  AliClusterContainer* clusters = 0;
676  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
677  clusters->SetClusPtCut(min);
678  }
679 }
680 
686 {
687  if (IsLocked()) return;
688 
689  TIter nextClusColl(&fClusterCollArray);
690  AliClusterContainer* clusters = 0;
691  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
692  clusters->SetClusECut(min);
693  }
694 }
695 
701 {
702  if (IsLocked()) return;
703 
704  TIter nextPartColl(&fParticleCollArray);
705  AliParticleContainer* tracks = 0;
706  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
707  tracks->SetParticlePtCut(min);
708  }
709 }
710 
717 {
718  if (IsLocked()) return;
719 
720  TIter nextPartColl(&fParticleCollArray);
721  AliParticleContainer* tracks = 0;
722  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
723  tracks->SetParticlePhiLimits(pmi, pma);
724  }
725 }
726 
733 fastjet::JetAlgorithm AliEmcalJetTask::ConvertToFJAlgo(EJetAlgo_t algo)
734 {
735  switch(algo) {
737  return fastjet::kt_algorithm;
739  return fastjet::antikt_algorithm;
741  return fastjet::cambridge_algorithm;
743  return fastjet::genkt_algorithm;
745  return fastjet::cambridge_for_passive_algorithm;
747  return fastjet::genkt_for_passive_algorithm;
749  return fastjet::plugin_algorithm;
751  return fastjet::undefined_jet_algorithm;
752 
753  default:
754  ::Error("AliEmcalJetTask::ConvertToFJAlgo", "Jet algorithm %d not recognized!!!", algo);
755  return fastjet::undefined_jet_algorithm;
756  }
757 }
758 
765 fastjet::RecombinationScheme AliEmcalJetTask::ConvertToFJRecoScheme(ERecoScheme_t reco)
766 {
767  switch(reco) {
769  return fastjet::E_scheme;
770 
772  return fastjet::pt_scheme;
773 
775  return fastjet::pt2_scheme;
776 
778  return fastjet::Et_scheme;
779 
781  return fastjet::Et2_scheme;
782 
784  return fastjet::BIpt_scheme;
785 
787  return fastjet::BIpt2_scheme;
788 
790  return fastjet::external_scheme;
791 
792  default:
793  ::Error("AliEmcalJetTask::ConvertToFJRecoScheme", "Recombination scheme %d not recognized!!!", reco);
794  return fastjet::external_scheme;
795  }
796 }
797 
803 
804  //This method has to be called after the run number is known because it needs the EMCal geometry object.
805 
806  UInt_t jetAcceptanceType = AliEmcalJet::kUser; // all jets satify the "no acceptance cut" condition
807 
808  // Check if TPC
809  if( eta < 0.9 && eta > -0.9 ) {
810  jetAcceptanceType |= AliEmcalJet::kTPC;
811  // Check if TPCfid
812  if (eta < 0.9 - r && eta > -0.9 + r)
813  jetAcceptanceType |= AliEmcalJet::kTPCfid;
814  }
815 
816  // Check if EMCAL
817  if( IsJetInEmcal(eta, phi, 0) ) {
818  jetAcceptanceType |= AliEmcalJet::kEMCAL;
819  // Check if EMCALfid
820  if( IsJetInEmcal(eta, phi, r) )
821  jetAcceptanceType |= AliEmcalJet::kEMCALfid;
822  }
823 
824  // Check if DCAL
825  if( IsJetInDcal(eta, phi, 0) ) {
826  jetAcceptanceType |= AliEmcalJet::kDCAL;
827  // Check if DCALfid
828  if( IsJetInDcal(eta, phi, r) )
829  jetAcceptanceType |= AliEmcalJet::kDCALfid;
830  }
831 
832  return jetAcceptanceType;
833 }
834 
839 {
840  if (!fGeom) return kFALSE;
841 
842  if ( eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
843  if(fRunNumber >= 177295 && fRunNumber <= 197470) {//small SM masked in 2012 and 2013
844  if ( phi < 3.135 - r && phi > 1.405 + r )
845  return kTRUE;
846  }
847  else {
848  if ( phi < fGeom->GetEMCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetArm1PhiMin() * TMath::DegToRad() + r)
849  return kTRUE;
850  }
851  }
852 
853  return kFALSE;
854 }
855 
860 {
861  if (!fGeom) return kFALSE;
862  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
863  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r)
864  return kTRUE;
865  }
866  return kFALSE;
867 }
868 
887  const TString nTracks, const TString nClusters,
888  const AliJetContainer::EJetAlgo_t jetAlgo, const Double_t radius, const AliJetContainer::EJetType_t jetType,
889  const Double_t minTrPt, const Double_t minClPt,
890  const Double_t ghostArea, const AliJetContainer::ERecoScheme_t reco,
891  const TString tag, const Double_t minJetPt,
892  const Bool_t lockTask, const Bool_t bFillGhosts
893 )
894 {
895  // Get the pointer to the existing analysis manager via the static access method
896  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
897  if (!mgr) {
898  ::Error("AddTaskEmcalJet", "No analysis manager to connect to.");
899  return 0;
900  }
901 
902  // Check the analysis type using the event handlers connected to the analysis manager
903  AliVEventHandler* handler = mgr->GetInputEventHandler();
904  if (!handler) {
905  ::Error("AddTaskEmcalJet", "This task requires an input event handler");
906  return 0;
907  }
908 
909  EDataType_t dataType = kUnknownDataType;
910 
911  if (handler->InheritsFrom("AliESDInputHandler")) {
912  dataType = kESD;
913  }
914  else if (handler->InheritsFrom("AliAODInputHandler")) {
915  dataType = kAOD;
916  }
917 
918  //-------------------------------------------------------
919  // Init the task and do settings
920  //-------------------------------------------------------
921 
922  TString trackName(nTracks);
923  TString clusName(nClusters);
924 
925  if (trackName == "usedefault") {
926  if (dataType == kESD) {
927  trackName = "Tracks";
928  }
929  else if (dataType == kAOD) {
930  trackName = "tracks";
931  }
932  else {
933  trackName = "";
934  }
935  }
936 
937  if (clusName == "usedefault") {
938  if (dataType == kESD) {
939  clusName = "CaloClusters";
940  }
941  else if (dataType == kAOD) {
942  clusName = "caloClusters";
943  }
944  else {
945  clusName = "";
946  }
947  }
948 
949  AliParticleContainer* partCont = 0;
950  if (trackName == "mcparticles") {
951  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
952  partCont = mcpartCont;
953  }
954  else if (trackName == "tracks" || trackName == "Tracks") {
955  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
956  partCont = trackCont;
957  }
958  else if (!trackName.IsNull()) {
959  partCont = new AliParticleContainer(trackName);
960  }
961  if (partCont) partCont->SetParticlePtCut(minTrPt);
962 
963  AliClusterContainer* clusCont = 0;
964  if (!clusName.IsNull()) {
965  clusCont = new AliClusterContainer(clusName);
966  clusCont->SetClusECut(0.);
967  clusCont->SetClusPtCut(0.);
968  clusCont->SetClusHadCorrEnergyCut(minClPt);
969  clusCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
970  }
971 
972  switch (jetType) {
974  if (partCont) partCont->SetCharge(AliParticleContainer::kCharged);
975  break;
977  if (partCont) partCont->SetCharge(AliParticleContainer::kNeutral);
978  break;
979  default:
980  break;
981  }
982 
983  TString name = AliJetContainer::GenerateJetName(jetType, jetAlgo, reco, radius, partCont, clusCont, tag);
984 
985  Printf("Jet task name: %s", name.Data());
986 
987  AliEmcalJetTask* mgrTask = static_cast<AliEmcalJetTask *>(mgr->GetTask(name.Data()));
988  if (mgrTask) return mgrTask;
989 
990  AliEmcalJetTask* jetTask = new AliEmcalJetTask(name);
991  jetTask->SetJetType(jetType);
992  jetTask->SetJetAlgo(jetAlgo);
993  jetTask->SetRecombScheme(reco);
994  jetTask->SetRadius(radius);
995  if (partCont) jetTask->AdoptParticleContainer(partCont);
996  if (clusCont) jetTask->AdoptClusterContainer(clusCont);
997  jetTask->SetJetsName(tag);
998  jetTask->SetMinJetPt(minJetPt);
999  jetTask->SetGhostArea(ghostArea);
1000 
1001  if (bFillGhosts) jetTask->SetFillGhost();
1002  if (lockTask) jetTask->SetLocked();
1003 
1004  // Final settings, pass to manager and set the containers
1005 
1006  mgr->AddTask(jetTask);
1007 
1008  // Create containers for input/output
1009  AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer();
1010  mgr->ConnectInput(jetTask, 0, cinput);
1011 
1012  TObjArray* cnt = mgr->GetContainers();
1013 
1014  return jetTask;
1015 }
Bool_t fTrackEfficiencyOnlyForEmbedding
void SetMaxNeutralPt(Double32_t t)
Definition: AliEmcalJet.h:184
void SetJetsName(const char *n)
void SetRecombScheme(ERecoScheme_t scheme)
fastjet::PseudoJet GetJetAreaVector(UInt_t idx) const
Double_t Area() const
Definition: AliEmcalJet.h:117
TObjArray fClusterCollArray
cluster collection array
void AddTrackAt(Int_t track, Int_t idx)
Definition: AliEmcalJet.h:197
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
DCal acceptance.
Definition: AliEmcalJet.h:62
void AdoptParticleContainer(AliParticleContainer *cont)
Double_t Eta() const
Definition: AliEmcalJet.h:108
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:104
Int_t GetNParticles() const
Container with name, TClonesArray and cuts for particles.
void FillJetConstituents(AliEmcalJet *jet, std::vector< fastjet::PseudoJet > &constituents, std::vector< fastjet::PseudoJet > &constituents_sub, Int_t flag=0, TClonesArray *particles_sub=0)
void AddClusterAt(Int_t clus, Int_t idx)
Definition: AliEmcalJet.h:196
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:178
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:64
TCanvas * c
Definition: TestFitELoss.C:172
void SetMaxRap(Double_t maxrap)
Definition: AliFJWrapper.h:100
virtual void Terminate(AliFJWrapper &fjw)=0
void SetMinJetPt(Double_t j)
void SetAreaE(Double_t a)
Definition: AliEmcalJet.h:181
AliEmcalContainerIndexMap< AliParticleContainer, AliVParticle > fParticleContainerIndexMap
! Mapping between index and particle containers
ERecoScheme_t fRecombScheme
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:96
Double_t fTrackEfficiency
void SetMCPt(Double_t p)
Definition: AliEmcalJet.h:191
Double_t E() const
Definition: AliEmcalJet.h:106
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:177
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:127
Container for particles within the EMCAL framework.
Bool_t fIsPSelSet
=true if already initialized
void SetR(Double_t r)
Definition: AliFJWrapper.h:101
void AdoptClusterContainer(AliClusterContainer *cont)
TObjArray fParticleCollArray
particle/track collection array
AliParticleContainer * GetParticleContainer(Int_t i=0) const
TPC fiducial acceptance (each eta edge narrowed by jet R)
Definition: AliEmcalJet.h:59
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
void SetJetAcceptanceType(UInt_t type)
Definition: AliEmcalJet.h:254
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:58
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
Double_t Phi_0_2pi() const
EMCal acceptance.
Definition: AliEmcalJet.h:60
void SetLegacyMode(Bool_t mode)
Definition: AliFJWrapper.h:111
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:95
static const AliEmcalContainerIndexMap< TClonesArray, AliVCluster > & GetEmcalContainerIndexMap()
Get the EMCal container utils associated with particle containers.
AliClusterContainer * GetClusterContainer(Int_t i=0) const
AliEmcalContainerIndexMap< AliClusterContainer, AliVCluster > fClusterContainerIndexMap
contituent index shift
Int_t GetNClusters() const
Double_t Phi_0_2pi() const
Definition: AliEmcalJet.h:116
void SetMaxChargedPt(Double32_t t)
Definition: AliEmcalJet.h:185
virtual void Clear(const Option_t *="")
void SetNEF(Double_t nef)
Definition: AliEmcalJet.h:186
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:192
AliVCluster * GetCluster(Int_t i) const
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:96
void SelectCollisionCandidates(UInt_t offlineTriggerMask=AliVEvent::kMB)
void SetGhostArea(Double_t gharea)
Definition: AliFJWrapper.h:99
void SetGhostArea(Double_t gharea)
EJetAlgo_t fJetAlgo
void SetNumberOfCharged(Int_t n)
Definition: AliEmcalJet.h:189
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
void SetAreaEta(Double_t a)
Definition: AliEmcalJet.h:179
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)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
void SetFillGhost(Bool_t b=kTRUE)
const AliParticleIterableMomentumContainer accepted_momentum() const
void CopyMappingFrom(const AliEmcalContainerIndexMap< U2, V > &map, U *cont)
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:188
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:63
Int_t GetIndexSub(Double_t phi_sub, std::vector< fastjet::PseudoJet > &constituents_unsub)
bool Bool_t
Definition: External.C:53
void SetPhiRange(Double_t pmi, Double_t pma)
EDataType_t
Switch for the data type.
void SetClusECut(Double_t cut)
void SetDefaultClusterEnergy(Int_t d)
void SetAreaPhi(Double_t a)
Definition: AliEmcalJet.h:180
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:97
void SetParticlePhiLimits(Double_t min, Double_t max)
void SetPtEmc(Double_t pt)
Definition: AliEmcalJet.h:193
Container structure for EMCAL clusters.
Container for MC-true particles within the EMCAL framework.
EMCal fiducial acceptance (each eta, phi edge narrowed by jet R)
Definition: AliEmcalJet.h:61
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:187
virtual void Init()=0
const Double_t phimin
void SetAxisInEmcal(Bool_t b)
Definition: AliEmcalJet.h:183
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:182
void SetClusHadCorrEnergyCut(Double_t cut)
void SetNumberOfNeutrals(Int_t n)
Definition: AliEmcalJet.h:190