AliPhysics  a34469b (a34469b)
 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  fUtilities(0),
71  fLocked(0),
72  fJetsName(),
73  fIsInit(0),
74  fIsPSelSet(0),
75  fIsEmcPart(0),
76  fLegacyMode(kFALSE),
77  fFillGhost(kFALSE),
78  fJets(0),
79  fFastJetWrapper("AliEmcalJetTask","AliEmcalJetTask")
80 {
81 }
82 
89  fJetsTag("Jets"),
90  fJetAlgo(AliJetContainer::antikt_algorithm),
91  fJetType(AliJetContainer::kFullJet),
92  fRecombScheme(AliJetContainer::pt_scheme),
93  fRadius(0.4),
94  fMinJetArea(0.001),
95  fMinJetPt(1.0),
96  fJetPhiMin(-10),
97  fJetPhiMax(+10),
98  fJetEtaMin(-1),
99  fJetEtaMax(+1),
100  fGhostArea(0.005),
101  fTrackEfficiency(1.),
102  fUtilities(0),
103  fLocked(0),
104  fJetsName(),
105  fIsInit(0),
106  fIsPSelSet(0),
107  fIsEmcPart(0),
108  fLegacyMode(kFALSE),
109  fFillGhost(kFALSE),
110  fJets(0),
111  fFastJetWrapper(name,name)
112 {
113 }
114 
119 {
120 }
121 
127 {
128  if (!fUtilities) fUtilities = new TObjArray();
129  if (fUtilities->FindObject(utility)) {
130  Error("AddUtility", "Jet utility %s already connected.", utility->GetName());
131  return utility;
132  }
133  fUtilities->Add(utility);
134  utility->SetJetTask(this);
135 
136  return utility;
137 }
138 
144 {
145  TIter next(fUtilities);
146  AliEmcalJetUtility *utility = 0;
147  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Init();
148 }
149 
156 {
157  TIter next(fUtilities);
158  AliEmcalJetUtility *utility = 0;
159  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Prepare(fFastJetWrapper);
160 }
161 
167 {
168  TIter next(fUtilities);
169  AliEmcalJetUtility *utility = 0;
170  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->ProcessJet(jet, ij, fFastJetWrapper);
171 }
172 
178 {
179  TIter next(fUtilities);
180  AliEmcalJetUtility *utility = 0;
181  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Terminate(fFastJetWrapper);
182 }
183 
189 {
190  // clear the jet array (normally a null operation)
191  fJets->Delete();
192 
193  Int_t n = FindJets();
194 
195  if (n == 0) return kFALSE;
196 
197  FillJetBranch();
198 
199  return kTRUE;
200 }
201 
210 {
211  if (fParticleCollArray.GetEntriesFast() == 0 && fClusterCollArray.GetEntriesFast() == 0){
212  AliError("No tracks or clusters, returning.");
213  return 0;
214  }
215 
217 
218  AliDebug(2,Form("Jet type = %d", fJetType));
219 
220  Int_t iColl = 1;
221  TIter nextPartColl(&fParticleCollArray);
222  AliParticleContainer* tracks = 0;
223  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
224  AliDebug(2,Form("Tracks from collection %d: '%s'.", iColl-1, tracks->GetName()));
226  for (AliParticleIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
227  // artificial inefficiency
228  if (fTrackEfficiency < 1.) {
229  Double_t rnd = gRandom->Rndm();
230  if (fTrackEfficiency < rnd) {
231  AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", it.current_index()));
232  continue;
233  }
234  }
235 
236  AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", it.current_index(), it->second->GetLabel(), it->first.Pt()));
237  Int_t uid = it.current_index() + fgkConstIndexShift * iColl;
238  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
239  }
240  iColl++;
241  }
242 
243  iColl = 1;
244  TIter nextClusColl(&fClusterCollArray);
245  AliClusterContainer* clusters = 0;
246  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
247  AliDebug(2,Form("Clusters from collection %d: '%s'.", iColl-1, clusters->GetName()));
249  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
250  AliDebug(2,Form("Cluster %d accepted (label = %d, energy = %.3f)", it.current_index(), it->second->GetLabel(), it->first.E()));
251  Int_t uid = -it.current_index() - fgkConstIndexShift * iColl;
252  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
253  }
254  iColl++;
255  }
256 
257  if (fFastJetWrapper.GetInputVectors().size() == 0) return 0;
258 
259  // run jet finder
261 
262  return fFastJetWrapper.GetInclusiveJets().size();
263 }
264 
271 {
273 
274  // loop over fastjet jets
275  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper.GetInclusiveJets();
276  // sort jets according to jet pt
277  static Int_t indexes[9999] = {-1};
278  GetSortedArray(indexes, jets_incl);
279 
280  AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
281  for (UInt_t ijet = 0, jetCount = 0; ijet < jets_incl.size(); ++ijet) {
282  Int_t ij = indexes[ijet];
283  AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fFastJetWrapper.GetJetArea(ij)));
284 
285  if (jets_incl[ij].perp() < fMinJetPt) continue;
286  if (fFastJetWrapper.GetJetArea(ij) < fMinJetArea) continue;
287  if ((jets_incl[ij].eta() < fJetEtaMin) || (jets_incl[ij].eta() > fJetEtaMax) ||
288  (jets_incl[ij].phi() < fJetPhiMin) || (jets_incl[ij].phi() > fJetPhiMax))
289  continue;
290 
291  AliEmcalJet *jet = new ((*fJets)[jetCount])
292  AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
293  jet->SetLabel(ij);
294 
295  fastjet::PseudoJet area(fFastJetWrapper.GetJetAreaVector(ij));
296  jet->SetArea(area.perp());
297  jet->SetAreaEta(area.eta());
298  jet->SetAreaPhi(area.phi());
299  jet->SetAreaE(area.E());
301 
302  // Fill constituent info
303  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper.GetJetConstituents(ij));
304  FillJetConstituents(jet, constituents, constituents);
305 
306  if (fGeom) {
307  if ((jet->Phi() > fGeom->GetArm1PhiMin() * TMath::DegToRad()) &&
308  (jet->Phi() < fGeom->GetArm1PhiMax() * TMath::DegToRad()) &&
309  (jet->Eta() > fGeom->GetArm1EtaMin()) &&
310  (jet->Eta() < fGeom->GetArm1EtaMax()))
311  jet->SetAxisInEmcal(kTRUE);
312  }
313 
314  ExecuteUtilities(jet, ij);
315 
316  AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), jet->GetNumberOfConstituents()));
317  jetCount++;
318  }
319 
321 }
322 
329 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
330 {
331  static Float_t pt[9999] = {0};
332 
333  const Int_t n = (Int_t)array.size();
334 
335  if (n < 1)
336  return kFALSE;
337 
338  for (Int_t i = 0; i < n; i++)
339  pt[i] = array[i].perp();
340 
341  TMath::Sort(n, pt, indexes);
342 
343  return kTRUE;
344 }
345 
352 {
353  if (fTrackEfficiency < 1.) {
354  if (gRandom) delete gRandom;
355  gRandom = new TRandom3(0);
356  }
357 
359 
360  // add jets to event if not yet there
361  if (!(InputEvent()->FindListObject(fJetsName))) {
362  fJets = new TClonesArray("AliEmcalJet");
363  fJets->SetName(fJetsName);
364  ::Info("AliEmcalJetTask::ExecOnce", "Jet collection with name '%s' has been added to the event.", fJetsName.Data());
365  InputEvent()->AddObject(fJets);
366  }
367  else {
368  AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
369  return;
370  }
371 
372  // setup fj wrapper
373  fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts);
379 
380  // setting legacy mode
381  if (fLegacyMode) {
383  }
384 
385  InitUtilities();
386 
388 }
389 
399 void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents,
400  std::vector<fastjet::PseudoJet>& constituents_unsub, Int_t flag, TClonesArray *particles_sub)
401 {
402  Int_t nt = 0;
403  Int_t nc = 0;
404  Double_t neutralE = 0.;
405  Double_t maxCh = 0.;
406  Double_t maxNe = 0.;
407  Int_t gall = 0;
408  Int_t gemc = 0;
409  Int_t cemc = 0;
410  Int_t ncharged = 0;
411  Int_t nneutral = 0;
412  Double_t mcpt = 0.;
413  Double_t emcpt = 0.;
414 
415  Int_t uid = -1;
416 
417  jet->SetNumberOfTracks(constituents.size());
418  jet->SetNumberOfClusters(constituents.size());
419 
420  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
421 
422  if (flag == 0) {
423  uid = constituents[ic].user_index();
424  }
425  else {
426  if (constituents[ic].perp()<1.e-10) {
427  uid=-1;
428  }
429  else {
430  uid = GetIndexSub(constituents[ic].phi(), constituents_unsub);
431  }
432  if (uid==0) {
433  AliError("correspondence between un/subtracted constituent not found");
434  continue;
435  }
436  }
437 
438  AliDebug(3,Form("Processing constituent %d", uid));
439  if (uid == -1) { //ghost particle
440  ++gall;
441  if (fGeom) {
442  Double_t gphi = constituents[ic].phi();
443  if (gphi < 0) gphi += TMath::TwoPi();
444  gphi *= TMath::RadToDeg();
445  Double_t geta = constituents[ic].eta();
446  if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) &&
447  (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax()))
448  ++gemc;
449  }
450 
451  if (fFillGhost) jet->AddGhost(constituents[ic].px(),
452  constituents[ic].py(),
453  constituents[ic].pz(),
454  constituents[ic].e());
455  }
456  else if (uid >= fgkConstIndexShift) { // track constituent
457  Int_t iColl = uid / fgkConstIndexShift;
458  Int_t tid = uid - iColl * fgkConstIndexShift;
459  iColl--;
460  AliDebug(3,Form("Constituent %d is a track from collection %d and with ID %d", uid, iColl, tid));
461  AliParticleContainer* partCont = GetParticleContainer(iColl);
462  if (!partCont) {
463  AliError(Form("Could not find particle container %d",iColl));
464  continue;
465  }
466  AliVParticle *t = partCont->GetParticle(tid);
467  if (!t) {
468  AliError(Form("Could not find track %d",tid));
469  continue;
470  }
471 
472  Double_t cEta = t->Eta();
473  Double_t cPhi = t->Phi();
474  Double_t cPt = t->Pt();
475  Double_t cP = t->P();
476  if (t->Charge() == 0) {
477  neutralE += cP;
478  ++nneutral;
479  if (cPt > maxNe) maxNe = cPt;
480  } else {
481  ++ncharged;
482  if (cPt > maxCh) maxCh = cPt;
483  }
484 
485  // check if MC particle
486  if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt;
487 
488  if (fGeom) {
489  if (cPhi < 0) cPhi += TMath::TwoPi();
490  cPhi *= TMath::RadToDeg();
491  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
492  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
493  emcpt += cPt;
494  ++cemc;
495  }
496  }
497 
498  if (flag == 0 || particles_sub == 0) {
499  jet->AddTrackAt(tid, nt);
500  }
501  else {
502  Int_t part_sub_id = particles_sub->GetEntriesFast();
503  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(dynamic_cast<AliVTrack*>(t)); // SA: probably need to be fixed!!
504  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
505  jet->AddTrackAt(part_sub_id, nt);
506  }
507 
508  ++nt;
509  }
510  else if (uid <= -fgkConstIndexShift) { // cluster constituent
511  Int_t iColl = -uid / fgkConstIndexShift;
512  Int_t cid = -uid - iColl * fgkConstIndexShift;
513  iColl--;
514  AliClusterContainer* clusCont = GetClusterContainer(iColl);
515  AliVCluster *c = clusCont->GetCluster(cid);
516 
517  if (!c) continue;
518 
520  clusCont->GetMomentum(nP, cid);
521 
522  Double_t cEta = nP.Eta();
523  Double_t cPhi = nP.Phi_0_2pi();
524  Double_t cPt = nP.Pt();
525  Double_t cP = nP.P();
526 
527  neutralE += cP;
528  if (cPt > maxNe) maxNe = cPt;
529 
530  // MC particle
531  if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
532 
533  if (fGeom) {
534  if (cPhi<0) cPhi += TMath::TwoPi();
535  cPhi *= TMath::RadToDeg();
536  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
537  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
538  emcpt += cPt;
539  ++cemc;
540  }
541  }
542 
543  if (flag == 0 || particles_sub == 0) {
544  jet->AddClusterAt(cid, nc);
545  }
546  else {
547  Int_t part_sub_id = particles_sub->GetEntriesFast();
548  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c);
549  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
550  jet->AddTrackAt(part_sub_id, nt);
551  }
552 
553  ++nc;
554  ++nneutral;
555  }
556  else {
557  AliError(Form("%s: No logical way to end up here (uid = %d).", GetName(), uid));
558  continue;
559  }
560  }
561 
562  jet->SetNumberOfTracks(nt);
563  jet->SetNumberOfClusters(nc);
564  jet->SetNEF(neutralE / jet->E());
565  jet->SetMaxChargedPt(maxCh);
566  jet->SetMaxNeutralPt(maxNe);
567  if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall);
568  else jet->SetAreaEmc(-1);
569  jet->SetNEmc(cemc);
570  jet->SetNumberOfCharged(ncharged);
571  jet->SetNumberOfNeutrals(nneutral);
572  jet->SetMCPt(mcpt);
573  jet->SetPtEmc(emcpt);
574  jet->SortConstituents();
575 }
576 
583 Int_t AliEmcalJetTask::GetIndexSub(Double_t phi_sub, std::vector<fastjet::PseudoJet>& constituents_unsub)
584 {
585  Double_t dphi=0;
586  Double_t phimin=100;
587  Int_t index=0;
588  for (UInt_t ii = 0; ii < constituents_unsub.size(); ii++) {
589  dphi=0;
590  Double_t phi_unsub = constituents_unsub[ii].phi();
591  dphi=phi_unsub-phi_sub;
592  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
593  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
594  if(TMath::Abs(dphi)<phimin){ phimin=TMath::Abs(dphi);
595  index=ii;} }
596  if (constituents_unsub[index].user_index()!=-1) return constituents_unsub[index].user_index();
597 
598  return 0;
599 }
600 
608 {
609  if (fLocked) {
610  AliFatal("Jet finder task is locked! Changing properties is not allowed.");
611  return kTRUE;
612  }
613  else {
614  return kFALSE;
615  }
616 }
617 
625 {
626  if(!fIsPSelSet) {
627  fIsPSelSet = kTRUE;
628  fOfflineTriggerMask = offlineTriggerMask;
629  }
630  else {
631  AliWarning("Manually setting the event selection for jet finders is not allowed! Using trigger=old_trigger | your_trigger");
632  fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask;
633  }
634 }
635 
642 {
643  if (IsLocked()) return;
644 
645  TIter nextPartColl(&fParticleCollArray);
646  AliParticleContainer* tracks = 0;
647  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
648  tracks->SetParticleEtaLimits(emi, ema);
649  }
650 }
651 
657 {
658  if (IsLocked()) return;
659 
660  TIter nextClusColl(&fClusterCollArray);
661  AliClusterContainer* clusters = 0;
662  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
663  clusters->SetClusPtCut(min);
664  }
665 }
666 
672 {
673  if (IsLocked()) return;
674 
675  TIter nextClusColl(&fClusterCollArray);
676  AliClusterContainer* clusters = 0;
677  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
678  clusters->SetClusECut(min);
679  }
680 }
681 
687 {
688  if (IsLocked()) return;
689 
690  TIter nextPartColl(&fParticleCollArray);
691  AliParticleContainer* tracks = 0;
692  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
693  tracks->SetParticlePtCut(min);
694  }
695 }
696 
703 {
704  if (IsLocked()) return;
705 
706  TIter nextPartColl(&fParticleCollArray);
707  AliParticleContainer* tracks = 0;
708  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
709  tracks->SetParticlePhiLimits(pmi, pma);
710  }
711 }
712 
719 fastjet::JetAlgorithm AliEmcalJetTask::ConvertToFJAlgo(EJetAlgo_t algo)
720 {
721  switch(algo) {
723  return fastjet::kt_algorithm;
725  return fastjet::antikt_algorithm;
727  return fastjet::cambridge_algorithm;
729  return fastjet::genkt_algorithm;
731  return fastjet::cambridge_for_passive_algorithm;
733  return fastjet::genkt_for_passive_algorithm;
735  return fastjet::plugin_algorithm;
737  return fastjet::undefined_jet_algorithm;
738 
739  default:
740  ::Error("AliEmcalJetTask::ConvertToFJAlgo", "Jet algorithm %d not recognized!!!", algo);
741  return fastjet::undefined_jet_algorithm;
742  }
743 }
744 
751 fastjet::RecombinationScheme AliEmcalJetTask::ConvertToFJRecoScheme(ERecoScheme_t reco)
752 {
753  switch(reco) {
755  return fastjet::E_scheme;
756 
758  return fastjet::pt_scheme;
759 
761  return fastjet::pt2_scheme;
762 
764  return fastjet::Et_scheme;
765 
767  return fastjet::Et2_scheme;
768 
770  return fastjet::BIpt_scheme;
771 
773  return fastjet::BIpt2_scheme;
774 
776  return fastjet::external_scheme;
777 
778  default:
779  ::Error("AliEmcalJetTask::ConvertToFJRecoScheme", "Recombination scheme %d not recognized!!!", reco);
780  return fastjet::external_scheme;
781  }
782 }
783 
789 
790  //This method has to be called after the run number is known because it needs the EMCal geometry object.
791 
792  UInt_t jetAcceptanceType = AliEmcalJet::kUser; // all jets satify the "no acceptance cut" condition
793 
794  // Check if TPC
795  if( eta < 0.9 && eta > -0.9 ) {
796  jetAcceptanceType |= AliEmcalJet::kTPC;
797  // Check if TPCfid
798  if (eta < 0.9 - r && eta > -0.9 + r)
799  jetAcceptanceType |= AliEmcalJet::kTPCfid;
800  }
801 
802  // Check if EMCAL
803  if( IsJetInEmcal(eta, phi, 0) ) {
804  jetAcceptanceType |= AliEmcalJet::kEMCAL;
805  // Check if EMCALfid
806  if( IsJetInEmcal(eta, phi, r) )
807  jetAcceptanceType |= AliEmcalJet::kEMCALfid;
808  }
809 
810  // Check if DCAL
811  if( IsJetInDcal(eta, phi, 0) ) {
812  jetAcceptanceType |= AliEmcalJet::kDCAL;
813  // Check if DCALfid
814  if( IsJetInDcal(eta, phi, r) )
815  jetAcceptanceType |= AliEmcalJet::kDCALfid;
816  }
817 
818  return jetAcceptanceType;
819 }
820 
825 {
826  if (!fGeom) return kFALSE;
827 
828  if ( eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
829  if(fRunNumber >= 177295 && fRunNumber <= 197470) {//small SM masked in 2012 and 2013
830  if ( phi < 3.135 - r && phi > 1.405 + r )
831  return kTRUE;
832  }
833  else {
834  if ( phi < fGeom->GetEMCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetArm1PhiMin() * TMath::DegToRad() + r)
835  return kTRUE;
836  }
837  }
838 
839  return kFALSE;
840 }
841 
846 {
847  if (!fGeom) return kFALSE;
848  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
849  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r)
850  return kTRUE;
851  }
852  return kFALSE;
853 }
854 
873  const TString nTracks, const TString nClusters,
874  const AliJetContainer::EJetAlgo_t jetAlgo, const Double_t radius, const AliJetContainer::EJetType_t jetType,
875  const Double_t minTrPt, const Double_t minClPt,
876  const Double_t ghostArea, const AliJetContainer::ERecoScheme_t reco,
877  const TString tag, const Double_t minJetPt,
878  const Bool_t lockTask, const Bool_t bFillGhosts
879 )
880 {
881  // Get the pointer to the existing analysis manager via the static access method
882  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
883  if (!mgr) {
884  ::Error("AddTaskEmcalJet", "No analysis manager to connect to.");
885  return 0;
886  }
887 
888  // Check the analysis type using the event handlers connected to the analysis manager
889  AliVEventHandler* handler = mgr->GetInputEventHandler();
890  if (!handler) {
891  ::Error("AddTaskEmcalJet", "This task requires an input event handler");
892  return 0;
893  }
894 
895  EDataType_t dataType = kUnknownDataType;
896 
897  if (handler->InheritsFrom("AliESDInputHandler")) {
898  dataType = kESD;
899  }
900  else if (handler->InheritsFrom("AliAODInputHandler")) {
901  dataType = kAOD;
902  }
903 
904  //-------------------------------------------------------
905  // Init the task and do settings
906  //-------------------------------------------------------
907 
908  TString trackName(nTracks);
909  TString clusName(nClusters);
910 
911  if (trackName == "usedefault") {
912  if (dataType == kESD) {
913  trackName = "Tracks";
914  }
915  else if (dataType == kAOD) {
916  trackName = "tracks";
917  }
918  else {
919  trackName = "";
920  }
921  }
922 
923  if (clusName == "usedefault") {
924  if (dataType == kESD) {
925  clusName = "CaloClusters";
926  }
927  else if (dataType == kAOD) {
928  clusName = "caloClusters";
929  }
930  else {
931  clusName = "";
932  }
933  }
934 
935  AliParticleContainer* partCont = 0;
936  if (trackName == "mcparticles") {
937  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
938  partCont = mcpartCont;
939  }
940  else if (trackName == "tracks" || trackName == "Tracks") {
941  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
942  partCont = trackCont;
943  }
944  else if (!trackName.IsNull()) {
945  partCont = new AliParticleContainer(trackName);
946  }
947  if (partCont) partCont->SetParticlePtCut(minTrPt);
948 
949  AliClusterContainer* clusCont = 0;
950  if (!clusName.IsNull()) {
951  clusCont = new AliClusterContainer(clusName);
952  clusCont->SetClusECut(0.);
953  clusCont->SetClusPtCut(0.);
954  clusCont->SetClusHadCorrEnergyCut(minClPt);
955  clusCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
956  }
957 
958  switch (jetType) {
960  if (partCont) partCont->SetCharge(AliParticleContainer::kCharged);
961  break;
963  if (partCont) partCont->SetCharge(AliParticleContainer::kNeutral);
964  break;
965  default:
966  break;
967  }
968 
969  TString name = AliJetContainer::GenerateJetName(jetType, jetAlgo, reco, radius, partCont, clusCont, tag);
970 
971  Printf("Jet task name: %s", name.Data());
972 
973  AliEmcalJetTask* mgrTask = static_cast<AliEmcalJetTask *>(mgr->GetTask(name.Data()));
974  if (mgrTask) return mgrTask;
975 
976  AliEmcalJetTask* jetTask = new AliEmcalJetTask(name);
977  jetTask->SetJetType(jetType);
978  jetTask->SetJetAlgo(jetAlgo);
979  jetTask->SetRecombScheme(reco);
980  jetTask->SetRadius(radius);
981  if (partCont) jetTask->AdoptParticleContainer(partCont);
982  if (clusCont) jetTask->AdoptClusterContainer(clusCont);
983  jetTask->SetJetsName(tag);
984  jetTask->SetMinJetPt(minJetPt);
985  jetTask->SetGhostArea(ghostArea);
986 
987  if (bFillGhosts) jetTask->SetFillGhost();
988  if (lockTask) jetTask->SetLocked();
989 
990  // Final settings, pass to manager and set the containers
991 
992  mgr->AddTask(jetTask);
993 
994  // Create containers for input/output
995  AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer();
996  mgr->ConnectInput(jetTask, 0, cinput);
997 
998  TObjArray* cnt = mgr->GetContainers();
999 
1000  return jetTask;
1001 }
void SetMaxNeutralPt(Double32_t t)
Definition: AliEmcalJet.h:182
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:195
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
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:194
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:176
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:179
ERecoScheme_t fRecombScheme
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:96
Double_t fTrackEfficiency
void SetMCPt(Double_t p)
Definition: AliEmcalJet.h:189
Double_t E() const
Definition: AliEmcalJet.h:106
void SetJetAlgo(EJetAlgo_t a)
virtual ~AliEmcalJetTask()
TRandom * gRandom
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
void SetLabel(Int_t l)
Definition: AliEmcalJet.h:175
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:252
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
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Double_t Phi_0_2pi() const
Definition: AliEmcalJet.h:116
void SetMaxChargedPt(Double32_t t)
Definition: AliEmcalJet.h:183
virtual void Clear(const Option_t *="")
void SetNEF(Double_t nef)
Definition: AliEmcalJet.h:184
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:190
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:187
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
void SetAreaEta(Double_t a)
Definition: AliEmcalJet.h:177
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
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:186
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:178
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:191
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:185
virtual void Init()=0
const Double_t phimin
void SetAxisInEmcal(Bool_t b)
Definition: AliEmcalJet.h:181
static const Int_t fgkConstIndexShift
fastjet wrapper
void SetAreaEmc(Double_t a)
Definition: AliEmcalJet.h:180
void SetClusHadCorrEnergyCut(Double_t cut)
void SetNumberOfNeutrals(Int_t n)
Definition: AliEmcalJet.h:188