AliPhysics  7f2a7c4 (7f2a7c4)
 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 }
154 {
155  TIter next(fUtilities);
156  AliEmcalJetUtility *utility = 0;
157  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->InitEvent(fFastJetWrapper);
158 }
159 
166 {
167  TIter next(fUtilities);
168  AliEmcalJetUtility *utility = 0;
169  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Prepare(fFastJetWrapper);
170 }
171 
177 {
178  TIter next(fUtilities);
179  AliEmcalJetUtility *utility = 0;
180  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->ProcessJet(jet, ij, fFastJetWrapper);
181 }
182 
188 {
189  TIter next(fUtilities);
190  AliEmcalJetUtility *utility = 0;
191  while ((utility=static_cast<AliEmcalJetUtility*>(next()))) utility->Terminate(fFastJetWrapper);
192 }
193 
199 {
200  InitEvent();
201  // clear the jet array (normally a null operation)
202  fJets->Delete();
203  Int_t n = FindJets();
204 
205  if (n == 0) return kFALSE;
206 
207  FillJetBranch();
208 
209  return kTRUE;
210 }
211 
220 {
221  if (fParticleCollArray.GetEntriesFast() == 0 && fClusterCollArray.GetEntriesFast() == 0){
222  AliError("No tracks or clusters, returning.");
223  return 0;
224  }
225 
227 
228  AliDebug(2,Form("Jet type = %d", fJetType));
229 
230  Int_t iColl = 1;
231  TIter nextPartColl(&fParticleCollArray);
232  AliParticleContainer* tracks = 0;
233  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
234  AliDebug(2,Form("Tracks from collection %d: '%s'.", iColl-1, tracks->GetName()));
236  for (AliParticleIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
237  // artificial inefficiency
238  if (fTrackEfficiency < 1.) {
239  Double_t rnd = gRandom->Rndm();
240  if (fTrackEfficiency < rnd) {
241  AliDebug(2,Form("Track %d rejected due to artificial tracking inefficiency", it.current_index()));
242  continue;
243  }
244  }
245 
246  AliDebug(2,Form("Track %d accepted (label = %d, pt = %f)", it.current_index(), it->second->GetLabel(), it->first.Pt()));
247  Int_t uid = it.current_index() + fgkConstIndexShift * iColl;
248  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
249  }
250  iColl++;
251  }
252 
253  iColl = 1;
254  TIter nextClusColl(&fClusterCollArray);
255  AliClusterContainer* clusters = 0;
256  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
257  AliDebug(2,Form("Clusters from collection %d: '%s'.", iColl-1, clusters->GetName()));
259  for (AliClusterIterableMomentumContainer::iterator it = itcont.begin(); it != itcont.end(); it++) {
260  AliDebug(2,Form("Cluster %d accepted (label = %d, energy = %.3f)", it.current_index(), it->second->GetLabel(), it->first.E()));
261  Int_t uid = -it.current_index() - fgkConstIndexShift * iColl;
262  fFastJetWrapper.AddInputVector(it->first.Px(), it->first.Py(), it->first.Pz(), it->first.E(), uid);
263  }
264  iColl++;
265  }
266 
267  if (fFastJetWrapper.GetInputVectors().size() == 0) return 0;
268 
269  // run jet finder
271 
272  return fFastJetWrapper.GetInclusiveJets().size();
273 }
274 
281 {
283 
284  // loop over fastjet jets
285  std::vector<fastjet::PseudoJet> jets_incl = fFastJetWrapper.GetInclusiveJets();
286  // sort jets according to jet pt
287  static Int_t indexes[9999] = {-1};
288  GetSortedArray(indexes, jets_incl);
289 
290  AliDebug(1,Form("%d jets found", (Int_t)jets_incl.size()));
291  for (UInt_t ijet = 0, jetCount = 0; ijet < jets_incl.size(); ++ijet) {
292  Int_t ij = indexes[ijet];
293  AliDebug(3,Form("Jet pt = %f, area = %f", jets_incl[ij].perp(), fFastJetWrapper.GetJetArea(ij)));
294 
295  if (jets_incl[ij].perp() < fMinJetPt) continue;
296  if (fFastJetWrapper.GetJetArea(ij) < fMinJetArea) continue;
297  if ((jets_incl[ij].eta() < fJetEtaMin) || (jets_incl[ij].eta() > fJetEtaMax) ||
298  (jets_incl[ij].phi() < fJetPhiMin) || (jets_incl[ij].phi() > fJetPhiMax))
299  continue;
300 
301  AliEmcalJet *jet = new ((*fJets)[jetCount])
302  AliEmcalJet(jets_incl[ij].perp(), jets_incl[ij].eta(), jets_incl[ij].phi(), jets_incl[ij].m());
303  jet->SetLabel(ij);
304 
305  fastjet::PseudoJet area(fFastJetWrapper.GetJetAreaVector(ij));
306  jet->SetArea(area.perp());
307  jet->SetAreaEta(area.eta());
308  jet->SetAreaPhi(area.phi());
309  jet->SetAreaE(area.E());
311 
312  // Fill constituent info
313  std::vector<fastjet::PseudoJet> constituents(fFastJetWrapper.GetJetConstituents(ij));
314  FillJetConstituents(jet, constituents, constituents);
315 
316  if (fGeom) {
317  if ((jet->Phi() > fGeom->GetArm1PhiMin() * TMath::DegToRad()) &&
318  (jet->Phi() < fGeom->GetArm1PhiMax() * TMath::DegToRad()) &&
319  (jet->Eta() > fGeom->GetArm1EtaMin()) &&
320  (jet->Eta() < fGeom->GetArm1EtaMax()))
321  jet->SetAxisInEmcal(kTRUE);
322  }
323 
324  ExecuteUtilities(jet, ij);
325 
326  AliDebug(2,Form("Added jet n. %d, pt = %f, area = %f, constituents = %d", jetCount, jet->Pt(), jet->Area(), jet->GetNumberOfConstituents()));
327  jetCount++;
328  }
329 
331 }
332 
339 Bool_t AliEmcalJetTask::GetSortedArray(Int_t indexes[], std::vector<fastjet::PseudoJet> array) const
340 {
341  static Float_t pt[9999] = {0};
342 
343  const Int_t n = (Int_t)array.size();
344 
345  if (n < 1)
346  return kFALSE;
347 
348  for (Int_t i = 0; i < n; i++)
349  pt[i] = array[i].perp();
350 
351  TMath::Sort(n, pt, indexes);
352 
353  return kTRUE;
354 }
355 
362 {
363  if (fTrackEfficiency < 1.) {
364  if (gRandom) delete gRandom;
365  gRandom = new TRandom3(0);
366  }
367 
369 
370  // add jets to event if not yet there
371  if (!(InputEvent()->FindListObject(fJetsName))) {
372  fJets = new TClonesArray("AliEmcalJet");
373  fJets->SetName(fJetsName);
374  ::Info("AliEmcalJetTask::ExecOnce", "Jet collection with name '%s' has been added to the event.", fJetsName.Data());
375  InputEvent()->AddObject(fJets);
376  }
377  else {
378  AliError(Form("%s: Object with name %s already in event! Returning", GetName(), fJetsName.Data()));
379  return;
380  }
381 
382  // setup fj wrapper
383  fFastJetWrapper.SetAreaType(fastjet::active_area_explicit_ghosts);
389 
390 
391  // setting legacy mode
392  if (fLegacyMode) {
394  }
395 
396  InitUtilities();
397 
398 
400 }
401 
411 void AliEmcalJetTask::FillJetConstituents(AliEmcalJet *jet, std::vector<fastjet::PseudoJet>& constituents,
412  std::vector<fastjet::PseudoJet>& constituents_unsub, Int_t flag, TClonesArray *particles_sub)
413 {
414  Int_t nt = 0;
415  Int_t nc = 0;
416  Double_t neutralE = 0.;
417  Double_t maxCh = 0.;
418  Double_t maxNe = 0.;
419  Int_t gall = 0;
420  Int_t gemc = 0;
421  Int_t cemc = 0;
422  Int_t ncharged = 0;
423  Int_t nneutral = 0;
424  Double_t mcpt = 0.;
425  Double_t emcpt = 0.;
426 
427  Int_t uid = -1;
428 
429  jet->SetNumberOfTracks(constituents.size());
430  jet->SetNumberOfClusters(constituents.size());
431 
432  for (UInt_t ic = 0; ic < constituents.size(); ++ic) {
433 
434  if (flag == 0) {
435  uid = constituents[ic].user_index();
436  }
437  else {
438  if (constituents[ic].perp()<1.e-10) {
439  uid=-1;
440  }
441  else {
442  uid = constituents[ic].user_index();
443  }
444  if (uid==0) {
445  AliError("correspondence between un/subtracted constituent not found");
446  continue;
447  }
448  }
449 
450  AliDebug(3,Form("Processing constituent %d", uid));
451  if (uid == -1) { //ghost particle
452  ++gall;
453  if (fGeom) {
454  Double_t gphi = constituents[ic].phi();
455  if (gphi < 0) gphi += TMath::TwoPi();
456  gphi *= TMath::RadToDeg();
457  Double_t geta = constituents[ic].eta();
458  if ((gphi > fGeom->GetArm1PhiMin()) && (gphi < fGeom->GetArm1PhiMax()) &&
459  (geta > fGeom->GetArm1EtaMin()) && (geta < fGeom->GetArm1EtaMax()))
460  ++gemc;
461  }
462 
463  if (fFillGhost) jet->AddGhost(constituents[ic].px(),
464  constituents[ic].py(),
465  constituents[ic].pz(),
466  constituents[ic].e());
467  }
468  else if (uid >= fgkConstIndexShift) { // track constituent
469  Int_t iColl = uid / fgkConstIndexShift;
470  Int_t tid = uid - iColl * fgkConstIndexShift;
471  iColl--;
472  AliDebug(3,Form("Constituent %d is a track from collection %d and with ID %d", uid, iColl, tid));
473  AliParticleContainer* partCont = GetParticleContainer(iColl);
474  if (!partCont) {
475  AliError(Form("Could not find particle container %d",iColl));
476  continue;
477  }
478  AliVParticle *t = partCont->GetParticle(tid);
479  if (!t) {
480  AliError(Form("Could not find track %d",tid));
481  continue;
482  }
483 
484  Double_t cEta = t->Eta();
485  Double_t cPhi = t->Phi();
486  Double_t cPt = t->Pt();
487  Double_t cP = t->P();
488  if (t->Charge() == 0) {
489  neutralE += cP;
490  ++nneutral;
491  if (cPt > maxNe) maxNe = cPt;
492  } else {
493  ++ncharged;
494  if (cPt > maxCh) maxCh = cPt;
495  }
496 
497  // check if MC particle
498  if (TMath::Abs(t->GetLabel()) > fMinMCLabel) mcpt += cPt;
499 
500  if (fGeom) {
501  if (cPhi < 0) cPhi += TMath::TwoPi();
502  cPhi *= TMath::RadToDeg();
503  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
504  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
505  emcpt += cPt;
506  ++cemc;
507  }
508  }
509 
510  if (flag == 0 || particles_sub == 0) {
511  jet->AddTrackAt(tid, nt);
512  }
513  else {
514  Int_t part_sub_id = particles_sub->GetEntriesFast();
515  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(dynamic_cast<AliVTrack*>(t)); // SA: probably need to be fixed!!
516  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
517  jet->AddTrackAt(part_sub_id, nt);
518  }
519 
520  ++nt;
521  }
522  else if (uid <= -fgkConstIndexShift) { // cluster constituent
523  Int_t iColl = -uid / fgkConstIndexShift;
524  Int_t cid = -uid - iColl * fgkConstIndexShift;
525  iColl--;
526  AliClusterContainer* clusCont = GetClusterContainer(iColl);
527  AliVCluster *c = clusCont->GetCluster(cid);
528 
529  if (!c) continue;
530 
532  clusCont->GetMomentum(nP, cid);
533 
534  Double_t cEta = nP.Eta();
535  Double_t cPhi = nP.Phi_0_2pi();
536  Double_t cPt = nP.Pt();
537  Double_t cP = nP.P();
538 
539  neutralE += cP;
540  if (cPt > maxNe) maxNe = cPt;
541 
542  // MC particle
543  if (TMath::Abs(c->GetLabel()) > fMinMCLabel) mcpt += c->GetMCEnergyFraction() > 1e-6 ? cPt * c->GetMCEnergyFraction() : cPt;
544 
545  if (fGeom) {
546  if (cPhi<0) cPhi += TMath::TwoPi();
547  cPhi *= TMath::RadToDeg();
548  if ((cPhi > fGeom->GetArm1PhiMin()) && (cPhi < fGeom->GetArm1PhiMax()) &&
549  (cEta > fGeom->GetArm1EtaMin()) && (cEta < fGeom->GetArm1EtaMax())) {
550  emcpt += cPt;
551  ++cemc;
552  }
553  }
554 
555  if (flag == 0 || particles_sub == 0) {
556  jet->AddClusterAt(cid, nc);
557  }
558  else {
559  Int_t part_sub_id = particles_sub->GetEntriesFast();
560  AliEmcalParticle* part_sub = new ((*particles_sub)[part_sub_id]) AliEmcalParticle(c);
561  part_sub->SetPtEtaPhiM(constituents[ic].perp(),constituents[ic].eta(),constituents[ic].phi(),constituents[ic].m());
562  jet->AddTrackAt(part_sub_id, nt);
563  }
564 
565  ++nc;
566  ++nneutral;
567  }
568  else {
569  AliError(Form("%s: No logical way to end up here (uid = %d).", GetName(), uid));
570  continue;
571  }
572  }
573 
574  jet->SetNumberOfTracks(nt);
575  jet->SetNumberOfClusters(nc);
576  jet->SetNEF(neutralE / jet->E());
577  jet->SetMaxChargedPt(maxCh);
578  jet->SetMaxNeutralPt(maxNe);
579  if (gall > 0) jet->SetAreaEmc(jet->Area() * gemc / gall);
580  else jet->SetAreaEmc(-1);
581  jet->SetNEmc(cemc);
582  jet->SetNumberOfCharged(ncharged);
583  jet->SetNumberOfNeutrals(nneutral);
584  jet->SetMCPt(mcpt);
585  jet->SetPtEmc(emcpt);
586  jet->SortConstituents();
587 }
588 
596 {
597  if (fLocked) {
598  AliFatal("Jet finder task is locked! Changing properties is not allowed.");
599  return kTRUE;
600  }
601  else {
602  return kFALSE;
603  }
604 }
605 
613 {
614  if(!fIsPSelSet) {
615  fIsPSelSet = kTRUE;
616  fOfflineTriggerMask = offlineTriggerMask;
617  }
618  else {
619  AliWarning("Manually setting the event selection for jet finders is not allowed! Using trigger=old_trigger | your_trigger");
620  fOfflineTriggerMask = fOfflineTriggerMask | offlineTriggerMask;
621  }
622 }
623 
630 {
631  if (IsLocked()) return;
632 
633  TIter nextPartColl(&fParticleCollArray);
634  AliParticleContainer* tracks = 0;
635  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
636  tracks->SetParticleEtaLimits(emi, ema);
637  }
638 }
639 
645 {
646  if (IsLocked()) return;
647 
648  TIter nextClusColl(&fClusterCollArray);
649  AliClusterContainer* clusters = 0;
650  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
651  clusters->SetClusPtCut(min);
652  }
653 }
654 
660 {
661  if (IsLocked()) return;
662 
663  TIter nextClusColl(&fClusterCollArray);
664  AliClusterContainer* clusters = 0;
665  while ((clusters = static_cast<AliClusterContainer*>(nextClusColl()))) {
666  clusters->SetClusECut(min);
667  }
668 }
669 
675 {
676  if (IsLocked()) return;
677 
678  TIter nextPartColl(&fParticleCollArray);
679  AliParticleContainer* tracks = 0;
680  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
681  tracks->SetParticlePtCut(min);
682  }
683 }
684 
691 {
692  if (IsLocked()) return;
693 
694  TIter nextPartColl(&fParticleCollArray);
695  AliParticleContainer* tracks = 0;
696  while ((tracks = static_cast<AliParticleContainer*>(nextPartColl()))) {
697  tracks->SetParticlePhiLimits(pmi, pma);
698  }
699 }
700 
707 fastjet::JetAlgorithm AliEmcalJetTask::ConvertToFJAlgo(EJetAlgo_t algo)
708 {
709  switch(algo) {
711  return fastjet::kt_algorithm;
713  return fastjet::antikt_algorithm;
715  return fastjet::cambridge_algorithm;
717  return fastjet::genkt_algorithm;
719  return fastjet::cambridge_for_passive_algorithm;
721  return fastjet::genkt_for_passive_algorithm;
723  return fastjet::plugin_algorithm;
725  return fastjet::undefined_jet_algorithm;
726 
727  default:
728  ::Error("AliEmcalJetTask::ConvertToFJAlgo", "Jet algorithm %d not recognized!!!", algo);
729  return fastjet::undefined_jet_algorithm;
730  }
731 }
732 
739 fastjet::RecombinationScheme AliEmcalJetTask::ConvertToFJRecoScheme(ERecoScheme_t reco)
740 {
741  switch(reco) {
743  return fastjet::E_scheme;
744 
746  return fastjet::pt_scheme;
747 
749  return fastjet::pt2_scheme;
750 
752  return fastjet::Et_scheme;
753 
755  return fastjet::Et2_scheme;
756 
758  return fastjet::BIpt_scheme;
759 
761  return fastjet::BIpt2_scheme;
762 
764  return fastjet::external_scheme;
765 
766  default:
767  ::Error("AliEmcalJetTask::ConvertToFJRecoScheme", "Recombination scheme %d not recognized!!!", reco);
768  return fastjet::external_scheme;
769  }
770 }
771 
777 
778  //This method has to be called after the run number is known because it needs the EMCal geometry object.
779 
780  UInt_t jetAcceptanceType = AliEmcalJet::kUser; // all jets satify the "no acceptance cut" condition
781 
782  // Check if TPC
783  if( eta < 0.9 && eta > -0.9 ) {
784  jetAcceptanceType |= AliEmcalJet::kTPC;
785  // Check if TPCfid
786  if (eta < 0.9 - r && eta > -0.9 + r)
787  jetAcceptanceType |= AliEmcalJet::kTPCfid;
788  }
789 
790  // Check if EMCAL
791  if( IsJetInEmcal(eta, phi, 0) ) {
792  jetAcceptanceType |= AliEmcalJet::kEMCAL;
793  // Check if EMCALfid
794  if( IsJetInEmcal(eta, phi, r) )
795  jetAcceptanceType |= AliEmcalJet::kEMCALfid;
796  }
797 
798  // Check if DCAL (i.e. eta-phi rectangle spanning DCal, which includes most of PHOS)
799  if( IsJetInDcal(eta, phi, 0) ) {
800  jetAcceptanceType |= AliEmcalJet::kDCAL;
801  // Check if DCALfid
802  if( IsJetInDcal(eta, phi, r) )
803  jetAcceptanceType |= AliEmcalJet::kDCALfid;
804  }
805 
806  // Check if DCALonly (i.e. ONLY DCal, does not include any of PHOS region)
807  if( IsJetInDcalOnly(eta, phi, 0) ) {
808  jetAcceptanceType |= AliEmcalJet::kDCALonly;
809  // Check if DCALonlyfid
810  if( IsJetInDcalOnly(eta, phi, r) )
811  jetAcceptanceType |= AliEmcalJet::kDCALonlyfid;
812  }
813 
814  // Check if PHOS
815  if( IsJetInPhos(eta, phi, 0) ) {
816  jetAcceptanceType |= AliEmcalJet::kPHOS;
817  // Check if PHOSfid
818  if( IsJetInPhos(eta, phi, r) )
819  jetAcceptanceType |= AliEmcalJet::kPHOSfid;
820  }
821 
822  return jetAcceptanceType;
823 }
824 
829 {
830  if (!fGeom) return kFALSE;
831 
832  if ( eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
833  if(fRunNumber >= 177295 && fRunNumber <= 197470) {//small SM masked in 2012 and 2013
834  if ( phi < 3.135 - r && phi > 1.405 + r )
835  return kTRUE;
836  }
837  else {
838  if ( phi < fGeom->GetEMCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetArm1PhiMin() * TMath::DegToRad() + r)
839  return kTRUE;
840  }
841  }
842 
843  return kFALSE;
844 }
845 
850 {
851  if (!fGeom) return kFALSE;
852  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r ) {
853  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r)
854  return kTRUE;
855  }
856  return kFALSE;
857 }
858 
866 {
867  if (!fGeom) return kFALSE;
868 
869  if (eta < fGeom->GetArm1EtaMax() - r && eta > fGeom->GetArm1EtaMin() + r) {
870  if ( phi < fGeom->GetDCALPhiMax() * TMath::DegToRad() - r && phi > fGeom->GetDCALPhiMin() * TMath::DegToRad() + r) {
871 
872  if (TMath::Abs(eta) > fGeom->GetDCALInnerExtandedEta() + r) {
873  return kTRUE;
874  }
875  if (r < 1e-6) {
876  if (phi > fGeom->GetEMCGeometry()->GetDCALStandardPhiMax() * TMath::DegToRad())
877  return kTRUE;
878  }
879 
880  }
881  }
882 
883  return kFALSE;
884 }
885 
890 {
891  Double_t etaMax = 0.130;
892  Double_t etaMin = -0.130;
893  Double_t phiMax = 320;
894  Double_t phiMin = 260; // Run 1
895  if (fRunNumber > 209121)
896  phiMin = 250; // Run 2
897 
898  if (eta < etaMax - r && eta > etaMin + r ) {
899  if (phi < phiMax * TMath::DegToRad() - r && phi > phiMin * TMath::DegToRad() + r)
900  return kTRUE;
901  }
902  return kFALSE;
903 }
904 
923  const TString nTracks, const TString nClusters,
924  const AliJetContainer::EJetAlgo_t jetAlgo, const Double_t radius, const AliJetContainer::EJetType_t jetType,
925  const Double_t minTrPt, const Double_t minClPt,
926  const Double_t ghostArea, const AliJetContainer::ERecoScheme_t reco,
927  const TString tag, const Double_t minJetPt,
928  const Bool_t lockTask, const Bool_t bFillGhosts
929 )
930 {
931  // Get the pointer to the existing analysis manager via the static access method
932  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
933  if (!mgr) {
934  ::Error("AddTaskEmcalJet", "No analysis manager to connect to.");
935  return 0;
936  }
937 
938  // Check the analysis type using the event handlers connected to the analysis manager
939  AliVEventHandler* handler = mgr->GetInputEventHandler();
940  if (!handler) {
941  ::Error("AddTaskEmcalJet", "This task requires an input event handler");
942  return 0;
943  }
944 
945  EDataType_t dataType = kUnknownDataType;
946 
947  if (handler->InheritsFrom("AliESDInputHandler")) {
948  dataType = kESD;
949  }
950  else if (handler->InheritsFrom("AliAODInputHandler")) {
951  dataType = kAOD;
952  }
953 
954  //-------------------------------------------------------
955  // Init the task and do settings
956  //-------------------------------------------------------
957 
958  TString trackName(nTracks);
959  TString clusName(nClusters);
960 
961  if (trackName == "usedefault") {
962  if (dataType == kESD) {
963  trackName = "Tracks";
964  }
965  else if (dataType == kAOD) {
966  trackName = "tracks";
967  }
968  else {
969  trackName = "";
970  }
971  }
972 
973  if (clusName == "usedefault") {
974  if (dataType == kESD) {
975  clusName = "CaloClusters";
976  }
977  else if (dataType == kAOD) {
978  clusName = "caloClusters";
979  }
980  else {
981  clusName = "";
982  }
983  }
984 
985  AliParticleContainer* partCont = 0;
986  if (trackName == "mcparticles") {
987  AliMCParticleContainer* mcpartCont = new AliMCParticleContainer(trackName);
988  partCont = mcpartCont;
989  }
990  else if (trackName == "tracks" || trackName == "Tracks") {
991  AliTrackContainer* trackCont = new AliTrackContainer(trackName);
992  partCont = trackCont;
993  }
994  else if (!trackName.IsNull()) {
995  partCont = new AliParticleContainer(trackName);
996  }
997  if (partCont) partCont->SetParticlePtCut(minTrPt);
998 
999  AliClusterContainer* clusCont = 0;
1000  if (!clusName.IsNull()) {
1001  clusCont = new AliClusterContainer(clusName);
1002  clusCont->SetClusECut(0.);
1003  clusCont->SetClusPtCut(0.);
1004  clusCont->SetClusHadCorrEnergyCut(minClPt);
1005  clusCont->SetDefaultClusterEnergy(AliVCluster::kHadCorr);
1006  }
1007 
1008  switch (jetType) {
1010  if (partCont) partCont->SetCharge(AliParticleContainer::kCharged);
1011  break;
1013  if (partCont) partCont->SetCharge(AliParticleContainer::kNeutral);
1014  break;
1015  default:
1016  break;
1017  }
1018 
1019  TString name = AliJetContainer::GenerateJetName(jetType, jetAlgo, reco, radius, partCont, clusCont, tag);
1020 
1021  Printf("Jet task name: %s", name.Data());
1022 
1023  AliEmcalJetTask* mgrTask = static_cast<AliEmcalJetTask *>(mgr->GetTask(name.Data()));
1024  if (mgrTask) return mgrTask;
1025 
1026  AliEmcalJetTask* jetTask = new AliEmcalJetTask(name);
1027  jetTask->SetJetType(jetType);
1028  jetTask->SetJetAlgo(jetAlgo);
1029  jetTask->SetRecombScheme(reco);
1030  jetTask->SetRadius(radius);
1031  if (partCont) jetTask->AdoptParticleContainer(partCont);
1032  if (clusCont) jetTask->AdoptClusterContainer(clusCont);
1033  jetTask->SetJetsName(tag);
1034  jetTask->SetMinJetPt(minJetPt);
1035  jetTask->SetGhostArea(ghostArea);
1036 
1037  if (bFillGhosts) jetTask->SetFillGhost();
1038  if (lockTask) jetTask->SetLocked();
1039 
1040  // Final settings, pass to manager and set the containers
1041 
1042  mgr->AddTask(jetTask);
1043 
1044  // Create containers for input/output
1045  AliAnalysisDataContainer* cinput = mgr->GetCommonInputContainer();
1046  mgr->ConnectInput(jetTask, 0, cinput);
1047 
1048  TObjArray* cnt = mgr->GetContainers();
1049 
1050  return jetTask;
1051 }
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:107
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:103
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:108
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:34
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:118
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:102
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:106
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
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
virtual void InitEvent(AliFJWrapper &fjw)=0
void SetAreaType(const fastjet::AreaType &atype)
Definition: AliFJWrapper.h:104
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)
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