AliPhysics  95775ff (95775ff)
 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 (i.e. eta-phi rectangle spanning DCal, which includes most of PHOS)
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  // Check if DCALonly (i.e. ONLY DCal, does not include any of PHOS region)
819  if( IsJetInDcalOnly(eta, phi, 0) ) {
820  jetAcceptanceType |= AliEmcalJet::kDCALonly;
821  // Check if DCALonlyfid
822  if( IsJetInDcalOnly(eta, phi, r) )
823  jetAcceptanceType |= AliEmcalJet::kDCALonlyfid;
824  }
825 
826  // Check if PHOS
827  if( IsJetInPhos(eta, phi, 0) ) {
828  jetAcceptanceType |= AliEmcalJet::kPHOS;
829  // Check if PHOSfid
830  if( IsJetInPhos(eta, phi, r) )
831  jetAcceptanceType |= AliEmcalJet::kPHOSfid;
832  }
833 
834  return jetAcceptanceType;
835 }
836 
841 {
842  if (!fGeom) return kFALSE;
843 
844  if ( eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
845  if(fRunNumber >= 177295 && fRunNumber <= 197470) {//small SM masked in 2012 and 2013
846  if ( phi < 3.135 - r && phi > 1.405 + r )
847  return kTRUE;
848  }
849  else {
850  if ( phi < fGeom->GetEMCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetArm1PhiMin() * TMath::DegToRad() + r)
851  return kTRUE;
852  }
853  }
854 
855  return kFALSE;
856 }
857 
862 {
863  if (!fGeom) return kFALSE;
864  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
865  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r)
866  return kTRUE;
867  }
868  return kFALSE;
869 }
870 
878 {
879  if (!fGeom) return kFALSE;
880 
881  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r) {
882  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r) {
883 
884  if (TMath::Abs(eta) > fGeom->GetDCALInnerExtandedEta() + r) {
885  return kTRUE;
886  }
887  if (r < 1e-6) {
888  if (phi > fGeom->GetEMCGeometry()->GetDCALStandardPhiMax() * TMath::DegToRad())
889  return kTRUE;
890  }
891 
892  }
893  }
894 
895  return kFALSE;
896 }
897 
902 {
903  Double_t etaMax = 0.130;
904  Double_t etaMin = -0.130;
905  Double_t phiMax = 320;
906  Double_t phiMin = 260; // Run 1
907  if (fRunNumber > 209121)
908  phiMin = 250; // Run 2
909 
910  if (eta < etaMax - r && eta > etaMin + r ) {
911  if (phi < phiMax * TMath::DegToRad() - r && phi > phiMin * TMath::DegToRad() + r)
912  return kTRUE;
913  }
914  return kFALSE;
915 }
916 
935  const TString nTracks, const TString nClusters,
936  const AliJetContainer::EJetAlgo_t jetAlgo, const Double_t radius, const AliJetContainer::EJetType_t jetType,
937  const Double_t minTrPt, const Double_t minClPt,
938  const Double_t ghostArea, const AliJetContainer::ERecoScheme_t reco,
939  const TString tag, const Double_t minJetPt,
940  const Bool_t lockTask, const Bool_t bFillGhosts
941 )
942 {
943  // Get the pointer to the existing analysis manager via the static access method
944  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
945  if (!mgr) {
946  ::Error("AddTaskEmcalJet", "No analysis manager to connect to.");
947  return 0;
948  }
949 
950  // Check the analysis type using the event handlers connected to the analysis manager
951  AliVEventHandler* handler = mgr->GetInputEventHandler();
952  if (!handler) {
953  ::Error("AddTaskEmcalJet", "This task requires an input event handler");
954  return 0;
955  }
956 
957  EDataType_t dataType = kUnknownDataType;
958 
959  if (handler->InheritsFrom("AliESDInputHandler")) {
960  dataType = kESD;
961  }
962  else if (handler->InheritsFrom("AliAODInputHandler")) {
963  dataType = kAOD;
964  }
965 
966  //-------------------------------------------------------
967  // Init the task and do settings
968  //-------------------------------------------------------
969 
970  TString trackName(nTracks);
971  TString clusName(nClusters);
972 
973  if (trackName == "usedefault") {
974  if (dataType == kESD) {
975  trackName = "Tracks";
976  }
977  else if (dataType == kAOD) {
978  trackName = "tracks";
979  }
980  else {
981  trackName = "";
982  }
983  }
984 
985  if (clusName == "usedefault") {
986  if (dataType == kESD) {
987  clusName = "CaloClusters";
988  }
989  else if (dataType == kAOD) {
990  clusName = "caloClusters";
991  }
992  else {
993  clusName = "";
994  }
995  }
996 
997  AliParticleContainer* partCont = 0;
998  if (trackName == "mcparticles") {
999  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
1000  partCont = mcpartCont;
1001  }
1002  else if (trackName == "tracks" || trackName == "Tracks") {
1003  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
1004  partCont = trackCont;
1005  }
1006  else if (!trackName.IsNull()) {
1007  partCont = new AliParticleContainer(trackName);
1008  }
1009  if (partCont) partCont->SetParticlePtCut(minTrPt);
1010 
1011  AliClusterContainer* clusCont = 0;
1012  if (!clusName.IsNull()) {
1013  clusCont = new AliClusterContainer(clusName);
1014  clusCont->SetClusECut(0.);
1015  clusCont->SetClusPtCut(0.);
1016  clusCont->SetClusHadCorrEnergyCut(minClPt);
1017  clusCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
1018  }
1019 
1020  switch (jetType) {
1022  if (partCont) partCont->SetCharge(AliParticleContainer::kCharged);
1023  break;
1025  if (partCont) partCont->SetCharge(AliParticleContainer::kNeutral);
1026  break;
1027  default:
1028  break;
1029  }
1030 
1031  TString name = AliJetContainer::GenerateJetName(jetType, jetAlgo, reco, radius, partCont, clusCont, tag);
1032 
1033  Printf("Jet task name: %s", name.Data());
1034 
1035  AliEmcalJetTask* mgrTask = static_cast<AliEmcalJetTask *>(mgr->GetTask(name.Data()));
1036  if (mgrTask) return mgrTask;
1037 
1038  AliEmcalJetTask* jetTask = new AliEmcalJetTask(name);
1039  jetTask->SetJetType(jetType);
1040  jetTask->SetJetAlgo(jetAlgo);
1041  jetTask->SetRecombScheme(reco);
1042  jetTask->SetRadius(radius);
1043  if (partCont) jetTask->AdoptParticleContainer(partCont);
1044  if (clusCont) jetTask->AdoptClusterContainer(clusCont);
1045  jetTask->SetJetsName(tag);
1046  jetTask->SetMinJetPt(minJetPt);
1047  jetTask->SetGhostArea(ghostArea);
1048 
1049  if (bFillGhosts) jetTask->SetFillGhost();
1050  if (lockTask) jetTask->SetLocked();
1051 
1052  // Final settings, pass to manager and set the containers
1053 
1054  mgr->AddTask(jetTask);
1055 
1056  // Create containers for input/output
1057  AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer();
1058  mgr->ConnectInput(jetTask, 0, cinput);
1059 
1060  TObjArray* cnt = mgr->GetContainers();
1061 
1062  return jetTask;
1063 }
void SetMaxNeutralPt(Double32_t t)
Definition: AliEmcalJet.h:188
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:201
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
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:200
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:182
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:100
virtual void Terminate(AliFJWrapper &fjw)=0
void SetMinJetPt(Double_t j)
void SetAreaE(Double_t a)
Definition: AliEmcalJet.h:185
ERecoScheme_t fRecombScheme
void SetRecombScheme(const fastjet::RecombinationScheme &scheme)
Definition: AliFJWrapper.h:96
Double_t fTrackEfficiency
void SetMCPt(Double_t p)
Definition: AliEmcalJet.h:195
Double_t E() const
Definition: AliEmcalJet.h:112
void SetJetAlgo(EJetAlgo_t a)
virtual ~AliEmcalJetTask()
TRandom * gRandom
static FJRecoScheme ConvertToFJRecoScheme(ERecoScheme_t reco)
void SetLabel(Int_t l)
Definition: AliEmcalJet.h:181
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: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:61
const std::vector< fastjet::PseudoJet > & GetInclusiveJets() const
Definition: AliFJWrapper.h:33
void SetJetAcceptanceType(UInt_t type)
Definition: AliEmcalJet.h:258
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: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:122
void SetMaxChargedPt(Double32_t t)
Definition: AliEmcalJet.h:189
virtual void Clear(const Option_t *="")
void SetNEF(Double_t nef)
Definition: AliEmcalJet.h:190
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:196
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:99
void SetGhostArea(Double_t gharea)
EJetAlgo_t fJetAlgo
void SetNumberOfCharged(Int_t n)
Definition: AliEmcalJet.h:193
std::vector< fastjet::PseudoJet > GetJetConstituents(UInt_t idx) const
void SetAreaEta(Double_t a)
Definition: AliEmcalJet.h:183
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:192
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
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.
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:184
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:197
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:191
virtual void Init()=0
Bool_t IsJetInPhos(Double_t eta, Double_t phi, Double_t r)
const Double_t phimin
void SetAxisInEmcal(Bool_t b)
Definition: AliEmcalJet.h:187
static const Int_t fgkConstIndexShift
fastjet wrapper
void SetAreaEmc(Double_t a)
Definition: AliEmcalJet.h:186
void SetClusHadCorrEnergyCut(Double_t cut)
void SetNumberOfNeutrals(Int_t n)
Definition: AliEmcalJet.h:194