AliPhysics  95775ff (95775ff)
 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  fEmcalTracks(0),
36  fEmcalClusters(0),
37  fNEmcalTracks(0),
38  fNEmcalClusters(0),
39  fHistMatchEtaAll(0),
40  fHistMatchPhiAll(0)
41 {
42  // Default constructor
43  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
44 
45  for(Int_t icent=0; icent<8; ++icent) {
46  for(Int_t ipt=0; ipt<9; ++ipt) {
47  for(Int_t ieta=0; ieta<2; ++ieta) {
48  fHistMatchEta[icent][ipt][ieta] = 0;
49  fHistMatchPhi[icent][ipt][ieta] = 0;
50  }
51  }
52  }
53 
54 }
55 
56 //________________________________________________________________________
58 {
59  // Destructor
60 }
61 
62 //________________________________________________________________________
64 {
65  // Initialization
66  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
68  // Do base class initializations and if it fails -> bail out
69  //AliAnalysisTaskEmcal::ExecOnce();
70  //if (!fInitialized) return;
71 
72  GetProperty("createHistos", fCreateHisto);
73 
74  GetProperty("maxDist", fMaxDistance);
75  GetProperty("updateClusters", fUpdateClusters);
76  GetProperty("updateTracks", fUpdateTracks);
78 
79  return kTRUE;
80 }
81 
82 //________________________________________________________________________
84 {
85  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
87 
88  fEmcalTracks = new TClonesArray("AliEmcalParticle");
89  fEmcalTracks->SetName(Form("EmcalTracks_%s", fPartCont->GetArrayName().Data()));
90  fEmcalClusters = new TClonesArray("AliEmcalParticle");
91  fEmcalClusters->SetName(Form("EmcalClusters_%s", fClusCont->GetArrayName().Data()));
92 
93  // Create my user objects.
94  if (fCreateHisto){
95  fHistMatchEtaAll = new TH1F("fHistMatchEtaAll", "fHistMatchEtaAll", 400, -0.2, 0.2);
96  fHistMatchPhiAll = new TH1F("fHistMatchPhiAll", "fHistMatchPhiAll", 400, -0.2, 0.2);
99 
100  const Int_t nCentChBins = fNcentBins * 2;
101  for(Int_t icent=0; icent<nCentChBins; ++icent) {
102  for(Int_t ipt=0; ipt<9; ++ipt) {
103  for(Int_t ieta=0; ieta<2; ++ieta) {
104  TString nameEta(Form("fHistMatchEta_%i_%i_%i",icent,ipt,ieta));
105  fHistMatchEta[icent][ipt][ieta] = new TH1F(nameEta, nameEta, 400, -0.2, 0.2);
106  fHistMatchEta[icent][ipt][ieta]->SetXTitle("#Delta#eta");
107  TString namePhi(Form("fHistMatchPhi_%i_%i_%i",icent,ipt,ieta));
108  fHistMatchPhi[icent][ipt][ieta] = new TH1F(namePhi, namePhi, 400, -0.2, 0.2);
109  fHistMatchPhi[icent][ipt][ieta]->SetXTitle("#Delta#phi");
110  fOutput->Add(fHistMatchEta[icent][ipt][ieta]);
111  fOutput->Add(fHistMatchPhi[icent][ipt][ieta]);
112  }
113  }
114  }
115  fOutput->SetOwner(kTRUE);
116  }
117 }
118 
119 //________________________________________________________________________
121 {
122  // Run
123  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
125 
126  // Run the matching.
128  DoMatching();
131 
132  return kTRUE;
133 }
134 
135 //________________________________________________________________________
137 {
138  // Get momenum bin.
139 
140  Int_t pbin=-1;
141  if (p<0.5)
142  pbin=0;
143  else if (p>=0.5 && p<1.0)
144  pbin=1;
145  else if (p>=1.0 && p<1.5)
146  pbin=2;
147  else if (p>=1.5 && p<2.)
148  pbin=3;
149  else if (p>=2. && p<3.)
150  pbin=4;
151  else if (p>=3. && p<4.)
152  pbin=5;
153  else if (p>=4. && p<5.)
154  pbin=6;
155  else if (p>=5. && p<8.)
156  pbin=7;
157  else if (p>=8.)
158  pbin=8;
159 
160  return pbin;
161 }
162 
163 //________________________________________________________________________
165 {
166  // Create AliEmcalParticle collections to handle the matching efficiently.
167  // At the same time propagates tracks, if requested.
168 
169  fEmcalTracks->Delete();
170  fEmcalClusters->Delete();
171 
172  fNEmcalTracks = 0;
173  fNEmcalClusters = 0;
174 
175  AliVCluster* cluster = 0;
176  AliVTrack* track = 0;
177 
178  fClusCont->ResetCurrentID();
179  while ((cluster = static_cast<AliVCluster*>(fClusCont->GetNextAcceptCluster()))) {
180 
181  // Clears the matching info
182  cluster->SetEmcCpvDistance(-1);
183  cluster->SetTrackDistance(1024, 1024);
184  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
185  AliESDCaloCluster *ec = 0;
186  if (ac) {
187  const Int_t N = ac->GetNTracksMatched();
188  AliDebug(2, TString::Format("Number of matched tracks: %d", N));
189  for (Int_t i = N - 1; i >= 0; i--) {
190  TObject *ptr = ac->GetTrackMatched(i);
191  ac->RemoveTrackMatched(ptr);
192  AliDebug(2, TString::Format("N tracks matched: %i of %i", ac->GetNTracksMatched(), N));
193  }
194  }
195  else {
196  ec = dynamic_cast<AliESDCaloCluster*>(cluster);
197  TArrayI *arr = ec->GetTracksMatched();
198  if (arr) arr->Set(0);
199  }
200 
201  // Create AliEmcalParticle objects to handle the matching
202  AliEmcalParticle* emcalCluster = new ((*fEmcalClusters)[fNEmcalClusters])
203  AliEmcalParticle(cluster, fClusCont->GetCurrentID(), fVertex[0], fVertex[1], fVertex[2], AliVCluster::kNonLinCorr);
204  emcalCluster->SetMatchedPtr(fEmcalTracks);
205 
206  fNEmcalClusters++;
207  }
208 
209  fPartCont->ResetCurrentID();
210  while ((track = static_cast<AliVTrack*>(fPartCont->GetNextAcceptParticle()))) {
211 
212  // Clears the matching info
213  track->ResetStatus(AliVTrack::kEMCALmatch);
214  track->SetEMCALcluster(-1);
215 
216  // Propagate tracks if requested
217  Bool_t propthistrack = kFALSE;
218  if (fDoPropagation) {
219  propthistrack = kTRUE;
220  }
221  else if (!track->IsExtrapolatedToEMCAL()) {
222  if (fAttemptProp) {
223  propthistrack = kTRUE;
224  }
225  else if (fAttemptPropMatch && IsTrackInEmcalAcceptance(track)) {
226  propthistrack = kTRUE;
227  }
228  }
229  if (propthistrack) AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(track, fPropDist);
230 
231  // Create AliEmcalParticle objects to handle the matching
232  AliEmcalParticle* emcalTrack = new ((*fEmcalTracks)[fNEmcalTracks]) AliEmcalParticle(track, fPartCont->GetCurrentID());
233  emcalTrack->SetMatchedPtr(fEmcalClusters);
234 
235  AliDebug(2, Form("Now adding track (pT = %.3f, eta = %.3f, phi = %.3f)"
236  "Phi, Eta on EMCal = %.3f, %.3f",
237  emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
238  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal()));
239 
240  fNEmcalTracks++;
241  }
242 }
243 
244 //________________________________________________________________________
246 {
247  // Set the links between tracks and clusters.
248  const Double_t maxd2 = fMaxDistance*fMaxDistance;
249 
250  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
251  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
252  AliVTrack* track = emcalTrack->GetTrack();
253 
254  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
255  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
256  AliVCluster* cluster = emcalCluster->GetCluster();
257 
258  Double_t deta = 999;
259  Double_t dphi = 999;
260  GetEtaPhiDiff(track, cluster, dphi, deta);
261  Double_t d2 = deta * deta + dphi * dphi;
262 
263  if (d2 > maxd2) continue;
264 
265  Double_t d = TMath::Sqrt(d2);
266  emcalCluster->AddMatchedObj(itrack, d);
267  emcalTrack->AddMatchedObj(icluster, d);
268  AliDebug(2, Form("Now matching cluster E = %.3f, pT = %.3f, eta = %.3f, phi = %.3f "
269  "with track pT = %.3f, eta = %.3f, phi = %.3f"
270  "Track eta, phi on EMCal = %.3f, %.3f, d = %.3f",
271  cluster->GetNonLinCorrEnergy(), emcalCluster->Pt(), emcalCluster->Eta(), emcalCluster->Phi(),
272  emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
273  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), d));
274 
275  if (fCreateHisto) {
276  Int_t mombin = GetMomBin(track->P());
277  Int_t centbinch = fCentBin;
278  if (track->Charge() < 0) centbinch += fNcentBins;
279  Int_t etabin = 0;
280  if(track->Eta() > 0) etabin = 1;
281 
282  fHistMatchEta[centbinch][mombin][etabin]->Fill(deta);
283  fHistMatchPhi[centbinch][mombin][etabin]->Fill(dphi);
284  fHistMatchEtaAll->Fill(deta);
285  fHistMatchPhiAll->Fill(dphi);
286  }
287  }
288  }
289 }
290 
291 //________________________________________________________________________
293 {
294  // Update clusters with matching info.
295 
296  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
297  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
298  const Int_t N = emcalCluster->GetNumberOfMatchedObj();
299  AliVCluster* cluster = emcalCluster->GetCluster();
300  AliDebug(3, Form("Cluster E = %.2f, eta = %.2f, phi = %.2f, Nmatch = %d", cluster->GetNonLinCorrEnergy(), emcalCluster->Eta(), emcalCluster->Phi(), N));
301 
302  if (N <= 0) continue;
303 
304  // Set the first match distance
305  const UInt_t firstMatchId = emcalCluster->GetMatchedObjId();
306  AliEmcalParticle* emcalTrackFirstMatch = static_cast<AliEmcalParticle*>(fEmcalTracks->At(firstMatchId));
307  AliVTrack* trackFirstMatch = emcalTrackFirstMatch->GetTrack();
308  Double_t deta = 999;
309  Double_t dphi = 999;
310  GetEtaPhiDiff(trackFirstMatch, cluster, dphi, deta);
311  cluster->SetTrackDistance(dphi, deta);
312 
313  // Cast into ESD/AOD objects
314  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
315  AliESDCaloCluster *ec = 0;
316  if (!ac) ec = dynamic_cast<AliESDCaloCluster*>(cluster);
317 
318  // Copy the matched tracks in the cluster. Note: different methods for ESD/AOD
319  if (ac) {
320  // TEMP
321  /*std::cout << "ProcessID for current process: " << TProcessID::GetPID()->GetName() << "/" << TProcessID::GetPID()->GetTitle() << ". Memory Address: " << TProcessID::GetPID() << std::endl;
322  std::cout << "N tracks matched: " << ac->GetNTracksMatched() << std::endl;
323  std::cout << "ProcessID for cluster: " << TProcessID::GetProcessWithUID(ac)->GetName() << "/" << TProcessID::GetProcessWithUID(ac)->GetTitle() << ". Memory Address: " << TProcessID::GetProcessWithUID(ac) << std::endl;*/
324  for (Int_t i=0; i < N; ++i) {
325  Int_t id = emcalCluster->GetMatchedObjId(i);
326  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
327 
328  AliDebug(3, Form("Pt = %.2f, eta = %.2f, phi = %.2f", emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi()));
329  // TEMP
330  //std::cout << "ProcessID for emcal track: " << TProcessID::GetProcessWithUID(emcalTrack)->GetName() << "/" << TProcessID::GetProcessWithUID(emcalTrack)->GetTitle() << ". Memory Address: " << TProcessID::GetProcessWithUID(emcalTrack) << std::endl;
331 
332  TObject *obj = emcalTrack->GetTrack();
333  // TEMP
334  // Superceded by assigning a new process ID earlier
335  //obj->SetBit(TObject::kIsReferenced);
336  //std::cout << "ProcessID for track: " << TProcessID::GetProcessWithUID(obj)->GetName() << "/" << TProcessID::GetProcessWithUID(obj)->GetTitle() << ". Memory Address: " << TProcessID::GetProcessWithUID(obj) << std::endl;
337  ac->AddTrackMatched(obj);
338  }
339  }
340  else {
341  TArrayI arr(N);
342  for (Int_t i = 0; i < N; ++i) {
343  Int_t id = emcalCluster->GetMatchedObjId(i);
344  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
345  arr.AddAt(emcalTrack->IdInCollection(), i);
346  }
347  ec->AddTracksMatched(arr);
348  }
349  }
350 }
351 
352 //________________________________________________________________________
354 {
355  // Update tracks with matching info.
356 
357  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
358  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
359 
360  if (emcalTrack->GetNumberOfMatchedObj() <= 0) continue;
361 
362  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(emcalTrack->GetMatchedObjId()));
363 
364  AliVTrack* track = emcalTrack->GetTrack();
365  track->SetEMCALcluster(emcalCluster->IdInCollection());
366  track->SetStatus(AliVTrack::kEMCALmatch);
367  }
368 }
369 
378 {
379 
380  if (!fGeom) {
381  AliWarning(Form("%s - AliAnalysisTaskEmcal::IsTrackInEmcalAcceptance - Geometry is not available!", GetName()));
382  return kFALSE;
383  }
384 
385  Double_t minPhi = fGeom->GetArm1PhiMin() - edges;
386  Double_t maxPhi = fGeom->GetArm1PhiMax() + edges;
387 
388  if (part->Phi() > minPhi && part->Phi() < maxPhi) {
389  return kTRUE;
390  }
391  else {
392  return kFALSE;
393  }
394 }
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
virtual AliVParticle * GetNextAcceptParticle()
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
UShort_t GetNumberOfMatchedObj() const
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
AliClusterContainer * fClusCont
Pointer to the cluster container.
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.
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
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
AliParticleContainer * fPartCont
Pointer to the track/particle container.
Double_t fMaxDistance
maximum distance to match clusters and tracks
AliVCluster * GetNextAcceptCluster()
static RegisterCorrectionComponent< AliEmcalCorrectionClusterTrackMatcher > reg
bool GetProperty(std::string propertyName, T &property, bool requiredProperty=true, std::string correctionName="")
Retrieve property.