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