AliPhysics  cdeda5a (cdeda5a)
 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 
21 
22 // Actually registers the class with the base class
24 
25 //________________________________________________________________________
27  AliEmcalCorrectionComponent("AliEmcalCorrectionClusterTrackMatcher"),
28  fPropDist(440),
29  fDoPropagation(kFALSE),
30  fAttemptProp(kTRUE),
31  fAttemptPropMatch(kFALSE),
32  fMaxDistance(0.1),
33  fUpdateTracks(kTRUE),
34  fUpdateClusters(kTRUE),
35  fClusterContainerIndexMap(),
36  fParticleContainerIndexMap(),
37  fEmcalTracks(0),
38  fEmcalClusters(0),
39  fNEmcalTracks(0),
40  fNEmcalClusters(0),
41  fHistMatchEtaAll(0),
42  fHistMatchPhiAll(0)
43 {
44  // Default constructor
45  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
46 
47  for(Int_t icent=0; icent<8; ++icent) {
48  for(Int_t ipt=0; ipt<9; ++ipt) {
49  for(Int_t ieta=0; ieta<2; ++ieta) {
50  fHistMatchEta[icent][ipt][ieta] = 0;
51  fHistMatchPhi[icent][ipt][ieta] = 0;
52  }
53  }
54  }
55 
56 }
57 
58 //________________________________________________________________________
60 {
61  // Destructor
62 }
63 
64 //________________________________________________________________________
66 {
67  // Initialization
68  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
70  // Do base class initializations and if it fails -> bail out
71  //AliAnalysisTaskEmcal::ExecOnce();
72  //if (!fInitialized) return;
73 
74  GetProperty("createHistos", fCreateHisto);
75 
76  GetProperty("maxDist", fMaxDistance);
77  GetProperty("updateClusters", fUpdateClusters);
78  GetProperty("updateTracks", fUpdateTracks);
80 
81  return kTRUE;
82 }
83 
84 //________________________________________________________________________
86 {
87  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
89 
90  // Determine all particle container array names for naming the AliEmcalParticle array
91  std::string particleContainerNames = "";
92  bool firstLoop = true;
93  AliParticleContainer * partCont = 0;
94  TIter nextPartCont(&fParticleCollArray);
95  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
96  if (firstLoop != true) {
97  particleContainerNames += "_";
98  }
99  else {
100  firstLoop = false;
101  }
102  particleContainerNames += partCont->GetArrayName();
103  }
104  // Determine all cluster container array names for naming the AliEmcalParticle array
105  std::string clusterContainerNames = "";
106  firstLoop = true;
107  AliClusterContainer * clusCont = 0;
108  TIter nextClusCont(&fClusterCollArray);
109  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
110  if (firstLoop != true) {
111  clusterContainerNames += "_";
112  }
113  else {
114  firstLoop = false;
115  }
116  clusterContainerNames += clusCont->GetArrayName();
117  }
118 
119  fEmcalTracks = new TClonesArray("AliEmcalParticle");
120  fEmcalTracks->SetName(Form("EmcalTracks_%s", particleContainerNames.c_str()));
121  fEmcalClusters = new TClonesArray("AliEmcalParticle");
122  fEmcalClusters->SetName(Form("EmcalClusters_%s", clusterContainerNames.c_str()));
123 
124  // Create my user objects.
125  if (fCreateHisto){
126  fHistMatchEtaAll = new TH1F("fHistMatchEtaAll", "fHistMatchEtaAll", 400, -0.2, 0.2);
127  fHistMatchPhiAll = new TH1F("fHistMatchPhiAll", "fHistMatchPhiAll", 400, -0.2, 0.2);
130 
131  const Int_t nCentChBins = fNcentBins * 2;
132  for(Int_t icent=0; icent<nCentChBins; ++icent) {
133  for(Int_t ipt=0; ipt<9; ++ipt) {
134  for(Int_t ieta=0; ieta<2; ++ieta) {
135  TString nameEta(Form("fHistMatchEta_%i_%i_%i",icent,ipt,ieta));
136  fHistMatchEta[icent][ipt][ieta] = new TH1F(nameEta, nameEta, 400, -0.2, 0.2);
137  fHistMatchEta[icent][ipt][ieta]->SetXTitle("#Delta#eta");
138  TString namePhi(Form("fHistMatchPhi_%i_%i_%i",icent,ipt,ieta));
139  fHistMatchPhi[icent][ipt][ieta] = new TH1F(namePhi, namePhi, 400, -0.2, 0.2);
140  fHistMatchPhi[icent][ipt][ieta]->SetXTitle("#Delta#phi");
141  fOutput->Add(fHistMatchEta[icent][ipt][ieta]);
142  fOutput->Add(fHistMatchPhi[icent][ipt][ieta]);
143  }
144  }
145  }
146  fOutput->SetOwner(kTRUE);
147  }
148 }
149 
150 //________________________________________________________________________
152 {
155 }
156 
157 //________________________________________________________________________
159 {
160  // Run
161  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
163 
164  // Run the matching.
166  DoMatching();
169 
170  return kTRUE;
171 }
172 
173 //________________________________________________________________________
175 {
176  // Get momenum bin.
177 
178  Int_t pbin=-1;
179  if (p<0.5)
180  pbin=0;
181  else if (p>=0.5 && p<1.0)
182  pbin=1;
183  else if (p>=1.0 && p<1.5)
184  pbin=2;
185  else if (p>=1.5 && p<2.)
186  pbin=3;
187  else if (p>=2. && p<3.)
188  pbin=4;
189  else if (p>=3. && p<4.)
190  pbin=5;
191  else if (p>=4. && p<5.)
192  pbin=6;
193  else if (p>=5. && p<8.)
194  pbin=7;
195  else if (p>=8.)
196  pbin=8;
197 
198  return pbin;
199 }
200 
201 //________________________________________________________________________
203 {
204  // Create AliEmcalParticle collections to handle the matching efficiently.
205  // At the same time propagates tracks, if requested.
206 
207  fEmcalTracks->Delete();
208  fEmcalClusters->Delete();
209 
210  fNEmcalTracks = 0;
211  fNEmcalClusters = 0;
212 
213  AliVCluster* cluster = 0;
214  AliVTrack* track = 0;
215 
216  AliClusterContainer * clusCont = 0;
217  TIter nextClusCont(&fClusterCollArray);
218  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
219  auto clusItCont = clusCont->accepted_momentum();
220  for (AliClusterIterableMomentumContainer::iterator clusIterator = clusItCont.begin(); clusIterator != clusItCont.end(); ++clusIterator) {
221  cluster = static_cast<AliVCluster *>(clusIterator->second);
222 
223  // Clears the matching info
224  cluster->SetEmcCpvDistance(-1);
225  cluster->SetTrackDistance(1024, 1024);
226  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
227  AliESDCaloCluster *ec = 0;
228  if (ac) {
229  const Int_t N = ac->GetNTracksMatched();
230  AliDebug(2, TString::Format("Number of matched tracks: %d", N));
231  for (Int_t i = N - 1; i >= 0; i--) {
232  TObject *ptr = ac->GetTrackMatched(i);
233  ac->RemoveTrackMatched(ptr);
234  AliDebug(2, TString::Format("N tracks matched: %i of %i", ac->GetNTracksMatched(), N));
235  }
236  }
237  else {
238  ec = dynamic_cast<AliESDCaloCluster*>(cluster);
239  TArrayI *arr = ec->GetTracksMatched();
240  if (arr) arr->Set(0);
241  }
242 
243  // Create AliEmcalParticle objects to handle the matching
244  AliEmcalParticle* emcalCluster = new ((*fEmcalClusters)[fNEmcalClusters])
245  AliEmcalParticle(cluster, fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fVertex[0], fVertex[1], fVertex[2], AliVCluster::kNonLinCorr);
246  emcalCluster->SetMatchedPtr(fEmcalTracks);
247 
248  fNEmcalClusters++;
249  }
250  }
251 
252  AliParticleContainer * partCont = 0;
253  TIter nextPartCont(&fParticleCollArray);
254  while ((partCont = static_cast<AliParticleContainer*>(nextPartCont()))) {
255  auto partItCont = partCont->accepted_momentum();
256  for (AliParticleIterableMomentumContainer::iterator partIterator = partItCont.begin(); partIterator != partItCont.end(); ++partIterator) {
257  track = static_cast<AliVTrack *>(partIterator->second);
258 
259  // Clears the matching info
260  track->ResetStatus(AliVTrack::kEMCALmatch);
261  track->SetEMCALcluster(-1);
262 
263  // Propagate tracks if requested
264  Bool_t propthistrack = kFALSE;
265  if (fDoPropagation) {
266  propthistrack = kTRUE;
267  }
268  else if (!track->IsExtrapolatedToEMCAL()) {
269  if (fAttemptProp) {
270  propthistrack = kTRUE;
271  }
272  else if (fAttemptPropMatch && IsTrackInEmcalAcceptance(track)) {
273  propthistrack = kTRUE;
274  }
275  }
276  if (propthistrack) AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(track, fPropDist);
277 
278  // Create AliEmcalParticle objects to handle the matching
279  AliEmcalParticle* emcalTrack = new ((*fEmcalTracks)[fNEmcalTracks]) AliEmcalParticle(track, fParticleContainerIndexMap.GlobalIndexFromLocalIndex(partCont, partIterator.current_index()));
280  emcalTrack->SetMatchedPtr(fEmcalClusters);
281 
282  AliDebug(2, Form("Now adding track (pT = %.3f, eta = %.3f, phi = %.3f)"
283  "Phi, Eta on EMCal = %.3f, %.3f",
284  emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
285  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal()));
286 
287  fNEmcalTracks++;
288  }
289  }
290 }
291 
292 //________________________________________________________________________
294 {
295  // Set the links between tracks and clusters.
296  const Double_t maxd2 = fMaxDistance*fMaxDistance;
297 
298  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
299  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
300  AliVTrack* track = emcalTrack->GetTrack();
301 
302  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
303  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
304  AliVCluster* cluster = emcalCluster->GetCluster();
305 
306  Double_t deta = 999;
307  Double_t dphi = 999;
308  GetEtaPhiDiff(track, cluster, dphi, deta);
309  Double_t d2 = deta * deta + dphi * dphi;
310 
311  if (d2 > maxd2) continue;
312 
313  Double_t d = TMath::Sqrt(d2);
314  emcalCluster->AddMatchedObj(itrack, d);
315  emcalTrack->AddMatchedObj(icluster, d);
316  AliDebug(2, Form("Now matching cluster E = %.3f, pT = %.3f, eta = %.3f, phi = %.3f "
317  "with track pT = %.3f, eta = %.3f, phi = %.3f"
318  "Track eta, phi on EMCal = %.3f, %.3f, d = %.3f",
319  cluster->GetNonLinCorrEnergy(), emcalCluster->Pt(), emcalCluster->Eta(), emcalCluster->Phi(),
320  emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
321  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), d));
322 
323  if (fCreateHisto) {
324  Int_t mombin = GetMomBin(track->P());
325  Int_t centbinch = fCentBin;
326  if (track->Charge() < 0) centbinch += fNcentBins;
327  Int_t etabin = 0;
328  if(track->Eta() > 0) etabin = 1;
329 
330  fHistMatchEta[centbinch][mombin][etabin]->Fill(deta);
331  fHistMatchPhi[centbinch][mombin][etabin]->Fill(dphi);
332  fHistMatchEtaAll->Fill(deta);
333  fHistMatchPhiAll->Fill(dphi);
334  }
335  }
336  }
337 }
338 
339 //________________________________________________________________________
341 {
342  // Update clusters with matching info.
343 
344  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
345  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
346  const Int_t N = emcalCluster->GetNumberOfMatchedObj();
347  AliVCluster* cluster = emcalCluster->GetCluster();
348  AliDebug(3, Form("Cluster E = %.2f, eta = %.2f, phi = %.2f, Nmatch = %d", cluster->GetNonLinCorrEnergy(), emcalCluster->Eta(), emcalCluster->Phi(), N));
349 
350  if (N <= 0) continue;
351 
352  // Set the first match distance
353  const UInt_t firstMatchId = emcalCluster->GetMatchedObjId();
354  AliEmcalParticle* emcalTrackFirstMatch = static_cast<AliEmcalParticle*>(fEmcalTracks->At(firstMatchId));
355  AliVTrack* trackFirstMatch = emcalTrackFirstMatch->GetTrack();
356  Double_t deta = 999;
357  Double_t dphi = 999;
358  GetEtaPhiDiff(trackFirstMatch, cluster, dphi, deta);
359  cluster->SetTrackDistance(dphi, deta);
360 
361  // Cast into ESD/AOD objects
362  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
363  AliESDCaloCluster *ec = 0;
364  if (!ac) ec = dynamic_cast<AliESDCaloCluster*>(cluster);
365 
366  // Copy the matched tracks in the cluster. Note: different methods for ESD/AOD
367  if (ac) {
368  // TEMP
369  /*std::cout << "ProcessID for current process: " << TProcessID::GetPID()->GetName() << "/" << TProcessID::GetPID()->GetTitle() << ". Memory Address: " << TProcessID::GetPID() << std::endl;
370  std::cout << "N tracks matched: " << ac->GetNTracksMatched() << std::endl;
371  std::cout << "ProcessID for cluster: " << TProcessID::GetProcessWithUID(ac)->GetName() << "/" << TProcessID::GetProcessWithUID(ac)->GetTitle() << ". Memory Address: " << TProcessID::GetProcessWithUID(ac) << std::endl;*/
372  for (Int_t i=0; i < N; ++i) {
373  Int_t id = emcalCluster->GetMatchedObjId(i);
374  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
375 
376  AliDebug(3, Form("Pt = %.2f, eta = %.2f, phi = %.2f", emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi()));
377  // TEMP
378  //std::cout << "ProcessID for emcal track: " << TProcessID::GetProcessWithUID(emcalTrack)->GetName() << "/" << TProcessID::GetProcessWithUID(emcalTrack)->GetTitle() << ". Memory Address: " << TProcessID::GetProcessWithUID(emcalTrack) << std::endl;
379 
380  TObject *obj = emcalTrack->GetTrack();
381  // TEMP
382  // Superceded by assigning a new process ID earlier
383  //obj->SetBit(TObject::kIsReferenced);
384  //std::cout << "ProcessID for track: " << TProcessID::GetProcessWithUID(obj)->GetName() << "/" << TProcessID::GetProcessWithUID(obj)->GetTitle() << ". Memory Address: " << TProcessID::GetProcessWithUID(obj) << std::endl;
385  ac->AddTrackMatched(obj);
386  }
387  }
388  else {
389  TArrayI arr(N);
390  for (Int_t i = 0; i < N; ++i) {
391  Int_t id = emcalCluster->GetMatchedObjId(i);
392  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
393  arr.AddAt(emcalTrack->IdInCollection(), i);
394  }
395  ec->AddTracksMatched(arr);
396  }
397  }
398 }
399 
400 //________________________________________________________________________
402 {
403  // Update tracks with matching info.
404 
405  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
406  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
407 
408  if (emcalTrack->GetNumberOfMatchedObj() <= 0) continue;
409 
410  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(emcalTrack->GetMatchedObjId()));
411 
412  AliVTrack* track = emcalTrack->GetTrack();
413  track->SetEMCALcluster(emcalCluster->IdInCollection());
414  track->SetStatus(AliVTrack::kEMCALmatch);
415  }
416 }
417 
426 {
427 
428  if (!fGeom) {
429  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
430  return kFALSE;
431  }
432 
433  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
434  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
435 
436  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
437  return kTRUE;
438  }
439  else {
440  return kFALSE;
441  }
442 }
Bool_t fUpdateClusters
update clusters with matching info
Int_t fNcentBins
How many centrality bins (this member copied from AliAnalysisTaskEmcal)
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
Bool_t fAttemptPropMatch
if true then attempt to propagate if not done yet but IsEMCAL is true
void SetMatchedPtr(TObjArray *arr)
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
TList * fOutput
! List of output histograms
Bool_t fCreateHisto
Flag to make some basic histograms.
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.