AliPhysics  1168478 (1168478)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliEmcalCorrectionClusterTrackMatcher.cxx
Go to the documentation of this file.
1 // AliEmcalCorrectionClusterTrackMatcher
2 //
3 
5 
6 #include <TH1.h>
7 #include <TList.h>
8 
9 #include "AliClusterContainer.h"
10 #include "AliParticleContainer.h"
11 #include "AliEMCALRecoUtils.h"
12 #include "AliESDCaloCluster.h"
13 #include "AliAODCaloCluster.h"
14 #include "AliVParticle.h"
15 #include "AliEmcalParticle.h"
16 #include "AliEMCALGeometry.h"
17 #include "AliMCEvent.h"
18 
22 
23 // Actually registers the class with the base class
25 
30  AliEmcalCorrectionComponent("AliEmcalCorrectionClusterTrackMatcher"),
31  fPropDist(440),
32  fDoPropagation(kFALSE),
33  fAttemptProp(kTRUE),
34  fAttemptPropMatch(kFALSE),
35  fMaxDistance(0.1),
36  fUsePIDmass(kTRUE),
37  fUseDCA(kTRUE),
38  fUpdateTracks(kTRUE),
39  fUpdateClusters(kTRUE),
40  fClusterContainerIndexMap(),
41  fParticleContainerIndexMap(),
42  fEmcalTracks(0),
43  fEmcalClusters(0),
44  fNEmcalTracks(0),
45  fNEmcalClusters(0),
46  fHistMatchEtaAll(0),
47  fHistMatchPhiAll(0),
48  fMCGenerToAcceptForTrack(1),
49  fNMCGenerToAccept(0)
50 {
51  for(Int_t icent=0; icent<8; ++icent) {
52  for(Int_t ipt=0; ipt<9; ++ipt) {
53  for(Int_t ieta=0; ieta<2; ++ieta) {
54  fHistMatchEta[icent][ipt][ieta] = 0;
55  fHistMatchPhi[icent][ipt][ieta] = 0;
56  }
57  }
58  }
59 
60  for(Int_t j = 0; j < 5; j++) fMCGenerToAccept[j] = "";
61 }
62 
67 {
68 }
69 
74 {
75  // Initialization
77 
78  GetProperty("createHistos", fCreateHisto);
79  GetProperty("usePIDmass", fUsePIDmass);
80  GetProperty("useDCA", fUseDCA);
81  GetProperty("maxDist", fMaxDistance);
82  GetProperty("updateClusters", fUpdateClusters);
83  GetProperty("updateTracks", fUpdateTracks);
85 
86  Bool_t enableFracEMCRecalc = kFALSE;
87  GetProperty("enableFracEMCRecalc", enableFracEMCRecalc);
88  Int_t removeNMCGenerators = 0;
89  GetProperty("removeNMCGenerators", removeNMCGenerators);
90  GetProperty("enableMCGenRemovTrack", fMCGenerToAcceptForTrack);
91  std::string removeMcGen1 = "";
92  GetProperty("removeMCGen1", removeMcGen1);
93  TString removeMCGen1 = removeMcGen1.c_str();
94  std::string removeMcGen2 = "";
95  GetProperty("removeMCGen2", removeMcGen2);
96  TString removeMCGen2 = removeMcGen2.c_str();
97 
98  if(enableFracEMCRecalc){
99  if(removeNMCGenerators > 0)
100  {
101  printf("\t gen1 <%s>, gen2 <%s>, remove tracks %d\n",removeMCGen1.Data(),removeMCGen2.Data(),fMCGenerToAcceptForTrack);
102  SetNumberOfMCGeneratorsToAccept(removeNMCGenerators) ;
103  SetNameOfMCGeneratorsToAccept(0,removeMCGen1);
104  SetNameOfMCGeneratorsToAccept(1,removeMCGen2);
105  }
106  }
107 
108  return kTRUE;
109 }
110 
115 {
117 
118  // Determine all particle container array names for naming the AliEmcalParticle array
119  std::string particleContainerNames = "";
120  bool firstLoop = true;
121  AliParticleContainer * partCont = 0;
122  TIter nextPartCont(&fParticleCollArray);
123  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
124  if (firstLoop != true) {
125  particleContainerNames += "_";
126  }
127  else {
128  firstLoop = false;
129  }
130  particleContainerNames += partCont->GetArrayName();
131  }
132  // Determine all cluster container array names for naming the AliEmcalParticle array
133  std::string clusterContainerNames = "";
134  firstLoop = true;
135  AliClusterContainer * clusCont = 0;
136  TIter nextClusCont(&fClusterCollArray);
137  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
138  if (firstLoop != true) {
139  clusterContainerNames += "_";
140  }
141  else {
142  firstLoop = false;
143  }
144  clusterContainerNames += clusCont->GetArrayName();
145  }
146 
147  fEmcalTracks = new TClonesArray("AliEmcalParticle");
148  fEmcalTracks->SetName(Form("EmcalTracks_%s", particleContainerNames.c_str()));
149  fEmcalClusters = new TClonesArray("AliEmcalParticle");
150  fEmcalClusters->SetName(Form("EmcalClusters_%s", clusterContainerNames.c_str()));
151 
152  // Create my user objects.
153  if (fCreateHisto){
154  fHistMatchEtaAll = new TH1F("fHistMatchEtaAll", "fHistMatchEtaAll", 400, -0.2, 0.2);
155  fHistMatchPhiAll = new TH1F("fHistMatchPhiAll", "fHistMatchPhiAll", 400, -0.2, 0.2);
158 
159  const Int_t nCentChBins = fNcentBins * 2;
160  for(Int_t icent=0; icent<nCentChBins; ++icent) {
161  for(Int_t ipt=0; ipt<9; ++ipt) {
162  for(Int_t ieta=0; ieta<2; ++ieta) {
163  TString nameEta(Form("fHistMatchEta_%i_%i_%i",icent,ipt,ieta));
164  fHistMatchEta[icent][ipt][ieta] = new TH1F(nameEta, nameEta, 400, -0.2, 0.2);
165  fHistMatchEta[icent][ipt][ieta]->SetXTitle("#Delta#eta");
166  TString namePhi(Form("fHistMatchPhi_%i_%i_%i",icent,ipt,ieta));
167  fHistMatchPhi[icent][ipt][ieta] = new TH1F(namePhi, namePhi, 400, -0.2, 0.2);
168  fHistMatchPhi[icent][ipt][ieta]->SetXTitle("#Delta#phi");
169  fOutput->Add(fHistMatchEta[icent][ipt][ieta]);
170  fOutput->Add(fHistMatchPhi[icent][ipt][ieta]);
171  }
172  }
173  }
174  fOutput->SetOwner(kTRUE);
175  }
176 }
177 
182 {
185 }
186 
191 {
193 
195 
196  // Run the matching.
198  DoMatching();
201 
202  return kTRUE;
203 }
204 
209 {
210  Int_t pbin=-1;
211  if (p<0.5)
212  pbin=0;
213  else if (p>=0.5 && p<1.0)
214  pbin=1;
215  else if (p>=1.0 && p<1.5)
216  pbin=2;
217  else if (p>=1.5 && p<2.)
218  pbin=3;
219  else if (p>=2. && p<3.)
220  pbin=4;
221  else if (p>=3. && p<4.)
222  pbin=5;
223  else if (p>=4. && p<5.)
224  pbin=6;
225  else if (p>=5. && p<8.)
226  pbin=7;
227  else if (p>=8.)
228  pbin=8;
229 
230  return pbin;
231 }
232 
238 {
239  fEmcalTracks->Delete();
240  fEmcalClusters->Delete();
241 
242  fNEmcalTracks = 0;
243  fNEmcalClusters = 0;
244 
245  AliVCluster* cluster = 0;
246  AliVTrack* track = 0;
247 
248  AliClusterContainer * clusCont = 0;
249  TIter nextClusCont(&fClusterCollArray);
250  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
251  auto clusItCont = clusCont->accepted_momentum();
252  for (AliClusterIterableMomentumContainer::iterator clusIterator = clusItCont.begin(); clusIterator != clusItCont.end(); ++clusIterator) {
253  cluster = static_cast<AliVCluster *>(clusIterator->second);
254 
255  // Clears the matching info
256  cluster->SetEmcCpvDistance(-1);
257  cluster->SetTrackDistance(1024, 1024);
258  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
259  AliESDCaloCluster *ec = 0;
260  if (ac) {
261  const Int_t N = ac->GetNTracksMatched();
262  AliDebug(2, TString::Format("Number of matched tracks: %d", N));
263  for (Int_t i = N - 1; i >= 0; i--) {
264  TObject *ptr = ac->GetTrackMatched(i);
265  ac->RemoveTrackMatched(ptr);
266  AliDebug(2, TString::Format("N tracks matched: %i of %i", ac->GetNTracksMatched(), N));
267  }
268  }
269  else {
270  ec = dynamic_cast<AliESDCaloCluster*>(cluster);
271  TArrayI *arr = ec->GetTracksMatched();
272  if (arr) arr->Set(0);
273  }
274 
275  // Create AliEmcalParticle objects to handle the matching
276  AliEmcalParticle* emcalCluster = new ((*fEmcalClusters)[fNEmcalClusters])
277  AliEmcalParticle(cluster, fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fVertex[0], fVertex[1], fVertex[2], AliVCluster::kNonLinCorr);
278  emcalCluster->SetMatchedPtr(fEmcalTracks);
279 
280  fNEmcalClusters++;
281  }
282  }
283 
284  Double_t mass;
285  if (fUsePIDmass) {
286  // use PID-based mass, and fUseDCA to determine starting point of propagation
287  mass = -1;
288  }
289  else {
290  // use pion mass, and fUseDCA to determine starting point of propagation
291  mass = 0.1396;
292  }
293 
294  AliParticleContainer * partCont = 0;
295  TIter nextPartCont(&fParticleCollArray);
296  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
297  auto partItCont = partCont->accepted_momentum();
298  for (AliParticleIterableMomentumContainer::iterator partIterator = partItCont.begin(); partIterator != partItCont.end(); ++partIterator) {
299  track = static_cast<AliVTrack *>(partIterator->second);
300 
301  // Clears the matching info
302  track->ResetStatus(AliVTrack::kEMCALmatch);
303  track->SetEMCALcluster(-1);
304 
305  // Propagate tracks if requested
306  Bool_t propthistrack = kFALSE;
307  if (fDoPropagation) {
308  propthistrack = kTRUE;
309  }
310  else if (!track->IsExtrapolatedToEMCAL()) {
311  if (fAttemptProp) {
312  propthistrack = kTRUE;
313  }
314  else if (fAttemptPropMatch && IsTrackInEmcalAcceptance(track)) {
315  propthistrack = kTRUE;
316  }
317  }
318  if (propthistrack) {
319 
320  // Check if track comes from a particular MC generator, do not include it if it is not a selected one
321  Int_t mcLabel = TMath::Abs(track->GetLabel());
322  TString genName;
324  {
325  fMCEvent->GetCocktailGenerator(mcLabel,genName);
326 
327  Bool_t generOK = kFALSE;
328  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
329  {
330  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
331  }
332 
333  if ( !generOK ) continue;
334  }
335 
336  // Propagate the track
337  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(track, fPropDist, mass, 20, 0.35, kFALSE, fUseDCA);
338  }
339 
340  // Reset properties of the track to fix TRefArray errors which occur when AddTrackMatched(obj) is called.
341  // This particular combination is from AODTrackFilterTask
342  // Resetting just the TProcessID of the track is not sufficient!
343  // It is included here because multiple track collections can come from different files with different
344  // Process IDs (for example, when embedding). As long as the updated tracks are not stored in a file,
345  // (normally they are not) then resetting these values won't adversely change the behavior of the code
346  // (except for fixing the TRefArray errors).
347  track->SetUniqueID(0);
348  track->ResetBit(TObject::kHasUUID);
349  track->ResetBit(TObject::kIsReferenced);
350 
351  // Create AliEmcalParticle objects to handle the matching
352  AliEmcalParticle* emcalTrack = new ((*fEmcalTracks)[fNEmcalTracks]) AliEmcalParticle(track, fParticleContainerIndexMap.GlobalIndexFromLocalIndex(partCont, partIterator.current_index()));
353  emcalTrack->SetMatchedPtr(fEmcalClusters);
354 
355  AliDebug(2, Form("Now adding track %i (pT = %.3f, eta = %.3f, phi = %.3f)"
356  "Phi, Eta on EMCal = %.3f, %.3f",
357  emcalTrack->IdInCollection(), emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
358  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal()));
359 
360  fNEmcalTracks++;
361  }
362  }
363 }
364 
369 {
370  const Double_t maxd2 = fMaxDistance*fMaxDistance;
371 
372  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
373  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
374  AliVTrack* track = emcalTrack->GetTrack();
375 
376  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
377  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
378  AliVCluster* cluster = emcalCluster->GetCluster();
379 
380  Double_t deta = 999;
381  Double_t dphi = 999;
382  GetEtaPhiDiff(track, cluster, dphi, deta);
383  Double_t d2 = deta * deta + dphi * dphi;
384 
385  if (d2 > maxd2) continue;
386 
387  Double_t d = TMath::Sqrt(d2);
388  emcalCluster->AddMatchedObj(itrack, d);
389  emcalTrack->AddMatchedObj(icluster, d);
390  AliDebug(2, Form("Now matching cluster E = %.3f, pT = %.3f, eta = %.3f, phi = %.3f "
391  "with track pT = %.3f, eta = %.3f, phi = %.3f"
392  "Track eta, phi on EMCal = %.3f, %.3f, d = %.3f",
393  cluster->GetNonLinCorrEnergy(), emcalCluster->Pt(), emcalCluster->Eta(), emcalCluster->Phi(),
394  emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
395  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), d));
396 
397  if (fCreateHisto) {
398  Int_t mombin = GetMomBin(track->P());
399  Int_t centbinch = fCentBin;
400  if (track->Charge() < 0) centbinch += fNcentBins;
401  Int_t etabin = 0;
402  if(track->Eta() > 0) etabin = 1;
403 
404  fHistMatchEta[centbinch][mombin][etabin]->Fill(deta);
405  fHistMatchPhi[centbinch][mombin][etabin]->Fill(dphi);
406  fHistMatchEtaAll->Fill(deta);
407  fHistMatchPhiAll->Fill(dphi);
408  }
409  }
410  }
411 }
412 
417 {
418  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
419  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
420  const Int_t N = emcalCluster->GetNumberOfMatchedObj();
421  AliVCluster* cluster = emcalCluster->GetCluster();
422  AliDebug(3, Form("Cluster E = %.2f, eta = %.2f, phi = %.2f, Nmatch = %d", cluster->GetNonLinCorrEnergy(), emcalCluster->Eta(), emcalCluster->Phi(), N));
423 
424  if (N <= 0) continue;
425 
426  // Set the first match distance
427  const UInt_t firstMatchId = emcalCluster->GetMatchedObjId();
428  AliEmcalParticle* emcalTrackFirstMatch = static_cast<AliEmcalParticle*>(fEmcalTracks->At(firstMatchId));
429  AliVTrack* trackFirstMatch = emcalTrackFirstMatch->GetTrack();
430  Double_t deta = 999;
431  Double_t dphi = 999;
432  GetEtaPhiDiff(trackFirstMatch, cluster, dphi, deta);
433  cluster->SetTrackDistance(dphi, deta);
434 
435  // Cast into ESD/AOD objects
436  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
437  AliESDCaloCluster *ec = 0;
438  if (!ac) ec = dynamic_cast<AliESDCaloCluster*>(cluster);
439 
440  // Copy the matched tracks in the cluster. Note: different methods for ESD/AOD
441  if (ac) {
442  for (Int_t i=0; i < N; ++i) {
443  Int_t id = emcalCluster->GetMatchedObjId(i);
444  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
445 
446  AliDebug(3, Form("Pt = %.2f, eta = %.2f, phi = %.2f", emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi()));
447 
448  TObject *obj = emcalTrack->GetTrack();
449  ac->AddTrackMatched(obj);
450  }
451  }
452  else {
453  TArrayI arr(N);
454  for (Int_t i = 0; i < N; ++i) {
455  Int_t id = emcalCluster->GetMatchedObjId(i);
456  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
457  arr.AddAt(emcalTrack->IdInCollection(), i);
458  }
459  ec->AddTracksMatched(arr);
460  }
461  }
462 }
463 
468 {
469  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
470  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
471 
472  if (emcalTrack->GetNumberOfMatchedObj() <= 0) continue;
473 
474  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(emcalTrack->GetMatchedObjId()));
475 
476  AliVTrack* track = emcalTrack->GetTrack();
477  track->SetEMCALcluster(emcalCluster->IdInCollection());
478  track->SetStatus(AliVTrack::kEMCALmatch);
479  }
480 }
481 
490 {
491 
492  if (!fGeom) {
493  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
494  return kFALSE;
495  }
496 
497  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
498  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
499 
500  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
501  return kTRUE;
502  }
503  else {
504  return kFALSE;
505  }
506 }
Bool_t fUpdateClusters
update clusters with matching info
Int_t fNcentBins
How many centrality bins (this member copied from AliAnalysisTaskEmcal)
Bool_t fUseDCA
Use DCA as starting point for track propagation, rather than primary vertex.
Double_t fPropDist
distance to surface (440cm default)
AliEMCALGeometry * fGeom
! Geometry object
Bool_t fAttemptProp
if true then attempt to propagate if not done yet
double Double_t
Definition: External.C:58
Bool_t fUpdateTracks
update tracks with matching info
Int_t GetMatchedObjId(UShort_t i=0) const
Double_t Phi() const
const AliClusterIterableMomentumContainer accepted_momentum() const
UShort_t GetNumberOfMatchedObj() const
bidirectional stl iterator over the EMCAL iterable container
Double_t mass
Bool_t fAttemptPropMatch
if true then attempt to propagate if not done yet but IsEMCAL is true
void SetMatchedPtr(TObjArray *arr)
Bool_t fUsePIDmass
Use PID-based mass hypothesis for track propagation, rather than pion mass hypothesis.
AliVTrack * GetTrack() const
Int_t fCentBin
! Event centrality bin
int GlobalIndexFromLocalIndex(const U *inputObject, const int localIndex) const
Container for particles within the EMCAL framework.
TObjArray fParticleCollArray
Particle/track collection array.
TObjArray fClusterCollArray
Cluster collection array.
AliEmcalContainerIndexMap< AliClusterContainer, AliVCluster > fClusterContainerIndexMap
! Mapping between index and cluster containers
int Int_t
Definition: External.C:63
Double_t Eta() const
Double_t Pt() const
unsigned int UInt_t
Definition: External.C:33
Cluster-track matcher component in the EMCal correction framework.
Base class for correction components in the EMCal correction framework.
static const AliEmcalContainerIndexMap< TClonesArray, AliVCluster > & GetEmcalContainerIndexMap()
Get the EMCal container utils associated with particle containers.
Int_t IdInCollection() const
Bool_t fDoPropagation
if true then propagate all hybrid tracks to EMCal surface
Double_t fVertex[3]
! Event vertex
TString fMCGenerToAccept[5]
List with name of generators that should not be included.
Bool_t fMCGenerToAcceptForTrack
Activate the removal of tracks entering the track matching that come from a particular generator...
TList * fOutput
! List of output histograms
Bool_t fCreateHisto
Flag to make some basic histograms.
Int_t fNMCGenerToAccept
Number of MC generators that should not be included in analysis.
AliVCluster * GetCluster() const
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
const AliParticleIterableMomentumContainer accepted_momentum() const
void CopyMappingFrom(const AliEmcalContainerIndexMap< U2, V > &map, U *cont)
void AddMatchedObj(Int_t id, Double_t d)
Bool_t IsTrackInEmcalAcceptance(AliVParticle *part, Double_t edges=0.9) const
void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
bool Bool_t
Definition: External.C:53
Container structure for EMCAL clusters.
Double_t fMaxDistance
maximum distance to match clusters and tracks
static RegisterCorrectionComponent< AliEmcalCorrectionClusterTrackMatcher > reg
AliEmcalContainerIndexMap< AliParticleContainer, AliVParticle > fParticleContainerIndexMap
! Mapping between index and particle containers
static const AliEmcalContainerIndexMap< TClonesArray, AliVParticle > & GetEmcalContainerIndexMap()
Get the EMCal container utils associated with particle containers.
bool GetProperty(std::string propertyName, T &property, bool requiredProperty=true, std::string correctionName="")
Retrieve property.