AliPhysics  63e47e1 (63e47e1)
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  fNMCGenerToAccept(0),
49  fMCGenerToAcceptForTrack(1)
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("usePIDmass", fUsePIDmass);
79  GetProperty("useDCA", fUseDCA);
80  GetProperty("maxDist", fMaxDistance);
81  GetProperty("updateClusters", fUpdateClusters);
82  GetProperty("updateTracks", fUpdateTracks);
84 
85  Bool_t enableFracEMCRecalc = kFALSE;
86  GetProperty("enableFracEMCRecalc", enableFracEMCRecalc);
87  Int_t removeNMCGenerators = 0;
88  GetProperty("removeNMCGenerators", removeNMCGenerators);
89  GetProperty("enableMCGenRemovTrack", fMCGenerToAcceptForTrack);
90  std::string removeMcGen1 = "";
91  GetProperty("removeMCGen1", removeMcGen1);
92  TString removeMCGen1 = removeMcGen1.c_str();
93  std::string removeMcGen2 = "";
94  GetProperty("removeMCGen2", removeMcGen2);
95  TString removeMCGen2 = removeMcGen2.c_str();
96 
97  if(enableFracEMCRecalc){
98  if(removeNMCGenerators > 0)
99  {
100  printf("\t gen1 <%s>, gen2 <%s>, remove tracks %d\n",removeMCGen1.Data(),removeMCGen2.Data(),fMCGenerToAcceptForTrack);
101  SetNumberOfMCGeneratorsToAccept(removeNMCGenerators) ;
102  SetNameOfMCGeneratorsToAccept(0,removeMCGen1);
103  SetNameOfMCGeneratorsToAccept(1,removeMCGen2);
104  }
105  }
106 
107  return kTRUE;
108 }
109 
114 {
116 
117  // Determine all particle container array names for naming the AliEmcalParticle array
118  std::string particleContainerNames = "";
119  bool firstLoop = true;
120  AliParticleContainer * partCont = 0;
121  TIter nextPartCont(&fParticleCollArray);
122  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
123  if (firstLoop != true) {
124  particleContainerNames += "_";
125  }
126  else {
127  firstLoop = false;
128  }
129  particleContainerNames += partCont->GetArrayName();
130  }
131  // Determine all cluster container array names for naming the AliEmcalParticle array
132  std::string clusterContainerNames = "";
133  firstLoop = true;
134  AliClusterContainer * clusCont = 0;
135  TIter nextClusCont(&fClusterCollArray);
136  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
137  if (firstLoop != true) {
138  clusterContainerNames += "_";
139  }
140  else {
141  firstLoop = false;
142  }
143  clusterContainerNames += clusCont->GetArrayName();
144  }
145 
146  fEmcalTracks = new TClonesArray("AliEmcalParticle");
147  fEmcalTracks->SetName(Form("EmcalTracks_%s", particleContainerNames.c_str()));
148  fEmcalClusters = new TClonesArray("AliEmcalParticle");
149  fEmcalClusters->SetName(Form("EmcalClusters_%s", clusterContainerNames.c_str()));
150 
151  // Create my user objects.
152  if (fCreateHisto){
153  fHistMatchEtaAll = new TH1F("fHistMatchEtaAll", "fHistMatchEtaAll", 400, -0.2, 0.2);
154  fHistMatchPhiAll = new TH1F("fHistMatchPhiAll", "fHistMatchPhiAll", 400, -0.2, 0.2);
157 
158  const Int_t nCentChBins = fNcentBins * 2;
159  for(Int_t icent=0; icent<nCentChBins; ++icent) {
160  for(Int_t ipt=0; ipt<9; ++ipt) {
161  for(Int_t ieta=0; ieta<2; ++ieta) {
162  TString nameEta(Form("fHistMatchEta_%i_%i_%i",icent,ipt,ieta));
163  fHistMatchEta[icent][ipt][ieta] = new TH1F(nameEta, nameEta, 400, -0.2, 0.2);
164  fHistMatchEta[icent][ipt][ieta]->SetXTitle("#Delta#eta");
165  TString namePhi(Form("fHistMatchPhi_%i_%i_%i",icent,ipt,ieta));
166  fHistMatchPhi[icent][ipt][ieta] = new TH1F(namePhi, namePhi, 400, -0.2, 0.2);
167  fHistMatchPhi[icent][ipt][ieta]->SetXTitle("#Delta#phi");
168  fOutput->Add(fHistMatchEta[icent][ipt][ieta]);
169  fOutput->Add(fHistMatchPhi[icent][ipt][ieta]);
170  }
171  }
172  }
173  fOutput->SetOwner(kTRUE);
174  }
175 }
176 
181 {
184 }
185 
190 {
192 
194 
195  // Run the matching.
197  DoMatching();
200 
201  return kTRUE;
202 }
203 
208 {
209  Int_t pbin=-1;
210  if (p<0.5)
211  pbin=0;
212  else if (p>=0.5 && p<1.0)
213  pbin=1;
214  else if (p>=1.0 && p<1.5)
215  pbin=2;
216  else if (p>=1.5 && p<2.)
217  pbin=3;
218  else if (p>=2. && p<3.)
219  pbin=4;
220  else if (p>=3. && p<4.)
221  pbin=5;
222  else if (p>=4. && p<5.)
223  pbin=6;
224  else if (p>=5. && p<8.)
225  pbin=7;
226  else if (p>=8.)
227  pbin=8;
228 
229  return pbin;
230 }
231 
237 {
238  fEmcalTracks->Delete();
239  fEmcalClusters->Delete();
240 
241  fNEmcalTracks = 0;
242  fNEmcalClusters = 0;
243 
244  AliVCluster* cluster = 0;
245  AliVTrack* track = 0;
246 
247  AliClusterContainer * clusCont = 0;
248  TIter nextClusCont(&fClusterCollArray);
249  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
250  auto clusItCont = clusCont->accepted_momentum();
251  for (AliClusterIterableMomentumContainer::iterator clusIterator = clusItCont.begin(); clusIterator != clusItCont.end(); ++clusIterator) {
252  cluster = static_cast<AliVCluster *>(clusIterator->second);
253 
254  // Clears the matching info
255  cluster->SetEmcCpvDistance(-1);
256  cluster->SetTrackDistance(1024, 1024);
257  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
258  AliESDCaloCluster *ec = 0;
259  if (ac) {
260  const Int_t N = ac->GetNTracksMatched();
261  AliDebug(2, TString::Format("Number of matched tracks: %d", N));
262  for (Int_t i = N - 1; i >= 0; i--) {
263  TObject *ptr = ac->GetTrackMatched(i);
264  ac->RemoveTrackMatched(ptr);
265  AliDebug(2, TString::Format("N tracks matched: %i of %i", ac->GetNTracksMatched(), N));
266  }
267  }
268  else {
269  ec = dynamic_cast<AliESDCaloCluster*>(cluster);
270  TArrayI *arr = ec->GetTracksMatched();
271  if (arr) arr->Set(0);
272  }
273 
274  // Create AliEmcalParticle objects to handle the matching
275  AliEmcalParticle* emcalCluster = new ((*fEmcalClusters)[fNEmcalClusters])
276  AliEmcalParticle(cluster, fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fVertex[0], fVertex[1], fVertex[2], AliVCluster::kNonLinCorr);
277  emcalCluster->SetMatchedPtr(fEmcalTracks);
278 
279  fNEmcalClusters++;
280  }
281  }
282 
283  Double_t mass;
284  if (fUsePIDmass) {
285  // use PID-based mass, and fUseDCA to determine starting point of propagation
286  mass = -1;
287  }
288  else {
289  // use pion mass, and fUseDCA to determine starting point of propagation
290  mass = 0.1396;
291  }
292 
293  AliParticleContainer * partCont = 0;
294  TIter nextPartCont(&fParticleCollArray);
295  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
296  auto partItCont = partCont->accepted_momentum();
297  for (AliParticleIterableMomentumContainer::iterator partIterator = partItCont.begin(); partIterator != partItCont.end(); ++partIterator) {
298  track = static_cast<AliVTrack *>(partIterator->second);
299 
300  // Clears the matching info
301  track->ResetStatus(AliVTrack::kEMCALmatch);
302  track->SetEMCALcluster(-1);
303 
304  // Propagate tracks if requested
305  Bool_t propthistrack = kFALSE;
306  if (fDoPropagation) {
307  propthistrack = kTRUE;
308  }
309  else if (!track->IsExtrapolatedToEMCAL()) {
310  if (fAttemptProp) {
311  propthistrack = kTRUE;
312  }
313  else if (fAttemptPropMatch && IsTrackInEmcalAcceptance(track)) {
314  propthistrack = kTRUE;
315  }
316  }
317  if (propthistrack) {
318 
319  // Check if track comes from a particular MC generator, do not include it if it is not a selected one
320  Int_t mcLabel = TMath::Abs(track->GetLabel());
321  TString genName;
323  {
324  fMCEvent->GetCocktailGenerator(mcLabel,genName);
325 
326  Bool_t generOK = kFALSE;
327  for(Int_t ig = 0; ig < fNMCGenerToAccept; ig++)
328  {
329  if ( genName.Contains(fMCGenerToAccept[ig]) ) generOK = kTRUE;
330  }
331 
332  if ( !generOK ) continue;
333  }
334 
335  // Propagate the track
336  AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(track, fPropDist, mass, 20, 0.35, kFALSE, fUseDCA);
337  }
338 
339  // Reset properties of the track to fix TRefArray errors which occur when AddTrackMatched(obj) is called.
340  // This particular combination is from AODTrackFilterTask
341  // Resetting just the TProcessID of the track is not sufficient!
342  // It is included here because multiple track collections can come from different files with different
343  // Process IDs (for example, when embedding). As long as the updated tracks are not stored in a file,
344  // (normally they are not) then resetting these values won't adversely change the behavior of the code
345  // (except for fixing the TRefArray errors).
346  track->SetUniqueID(0);
347  track->ResetBit(TObject::kHasUUID);
348  track->ResetBit(TObject::kIsReferenced);
349 
350  // Create AliEmcalParticle objects to handle the matching
351  AliEmcalParticle* emcalTrack = new ((*fEmcalTracks)[fNEmcalTracks]) AliEmcalParticle(track, fParticleContainerIndexMap.GlobalIndexFromLocalIndex(partCont, partIterator.current_index()));
352  emcalTrack->SetMatchedPtr(fEmcalClusters);
353 
354  AliDebug(2, Form("Now adding track %i (pT = %.3f, eta = %.3f, phi = %.3f)"
355  "Phi, Eta on EMCal = %.3f, %.3f",
356  emcalTrack->IdInCollection(), emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
357  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal()));
358 
359  fNEmcalTracks++;
360  }
361  }
362 }
363 
368 {
369  const Double_t maxd2 = fMaxDistance*fMaxDistance;
370 
371  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
372  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
373  AliVTrack* track = emcalTrack->GetTrack();
374 
375  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
376  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
377  AliVCluster* cluster = emcalCluster->GetCluster();
378 
379  Double_t deta = 999;
380  Double_t dphi = 999;
381  GetEtaPhiDiff(track, cluster, dphi, deta);
382  Double_t d2 = deta * deta + dphi * dphi;
383 
384  if (d2 > maxd2) continue;
385 
386  Double_t d = TMath::Sqrt(d2);
387  emcalCluster->AddMatchedObj(itrack, d);
388  emcalTrack->AddMatchedObj(icluster, d);
389  AliDebug(2, Form("Now matching cluster E = %.3f, pT = %.3f, eta = %.3f, phi = %.3f "
390  "with track pT = %.3f, eta = %.3f, phi = %.3f"
391  "Track eta, phi on EMCal = %.3f, %.3f, d = %.3f",
392  cluster->GetNonLinCorrEnergy(), emcalCluster->Pt(), emcalCluster->Eta(), emcalCluster->Phi(),
393  emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
394  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), d));
395 
396  if (fCreateHisto) {
397  Int_t mombin = GetMomBin(track->P());
398  Int_t centbinch = fCentBin;
399  if (track->Charge() < 0) centbinch += fNcentBins;
400  Int_t etabin = 0;
401  if(track->Eta() > 0) etabin = 1;
402 
403  fHistMatchEta[centbinch][mombin][etabin]->Fill(deta);
404  fHistMatchPhi[centbinch][mombin][etabin]->Fill(dphi);
405  fHistMatchEtaAll->Fill(deta);
406  fHistMatchPhiAll->Fill(dphi);
407  }
408  }
409  }
410 }
411 
416 {
417  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
418  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
419  const Int_t N = emcalCluster->GetNumberOfMatchedObj();
420  AliVCluster* cluster = emcalCluster->GetCluster();
421  AliDebug(3, Form("Cluster E = %.2f, eta = %.2f, phi = %.2f, Nmatch = %d", cluster->GetNonLinCorrEnergy(), emcalCluster->Eta(), emcalCluster->Phi(), N));
422 
423  if (N <= 0) continue;
424 
425  // Set the first match distance
426  const UInt_t firstMatchId = emcalCluster->GetMatchedObjId();
427  AliEmcalParticle* emcalTrackFirstMatch = static_cast<AliEmcalParticle*>(fEmcalTracks->At(firstMatchId));
428  AliVTrack* trackFirstMatch = emcalTrackFirstMatch->GetTrack();
429  Double_t deta = 999;
430  Double_t dphi = 999;
431  GetEtaPhiDiff(trackFirstMatch, cluster, dphi, deta);
432  cluster->SetTrackDistance(dphi, deta);
433 
434  // Cast into ESD/AOD objects
435  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
436  AliESDCaloCluster *ec = 0;
437  if (!ac) ec = dynamic_cast<AliESDCaloCluster*>(cluster);
438 
439  // Copy the matched tracks in the cluster. Note: different methods for ESD/AOD
440  if (ac) {
441  for (Int_t i=0; i < N; ++i) {
442  Int_t id = emcalCluster->GetMatchedObjId(i);
443  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
444 
445  AliDebug(3, Form("Pt = %.2f, eta = %.2f, phi = %.2f", emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi()));
446 
447  TObject *obj = emcalTrack->GetTrack();
448  ac->AddTrackMatched(obj);
449  }
450  }
451  else {
452  TArrayI arr(N);
453  for (Int_t i = 0; i < N; ++i) {
454  Int_t id = emcalCluster->GetMatchedObjId(i);
455  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
456  arr.AddAt(emcalTrack->IdInCollection(), i);
457  }
458  ec->AddTracksMatched(arr);
459  }
460  }
461 }
462 
467 {
468  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
469  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
470 
471  if (emcalTrack->GetNumberOfMatchedObj() <= 0) continue;
472 
473  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(emcalTrack->GetMatchedObjId()));
474 
475  AliVTrack* track = emcalTrack->GetTrack();
476  track->SetEMCALcluster(emcalCluster->IdInCollection());
477  track->SetStatus(AliVTrack::kEMCALmatch);
478  }
479 }
480 
489 {
490 
491  if (!fGeom) {
492  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
493  return kFALSE;
494  }
495 
496  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
497  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
498 
499  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
500  return kTRUE;
501  }
502  else {
503  return kFALSE;
504  }
505 }
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
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.