AliPhysics  a0db429 (a0db429)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
AliEmcalClusTrackMatcherTask.cxx
Go to the documentation of this file.
1 //
2 // Track/cluster matcher
3 //
4 // Author: C.Loizides, S.Aiola
5 
7 
8 #include <TClonesArray.h>
9 #include <TClass.h>
10 
11 #include <AliAODCaloCluster.h>
12 #include <AliESDCaloCluster.h>
13 #include <AliLog.h>
14 #include <AliVCluster.h>
15 #include <AliVTrack.h>
16 #include <AliEMCALRecoUtils.h>
17 
18 #include "AliEmcalParticle.h"
19 #include "AliParticleContainer.h"
20 #include "AliClusterContainer.h"
21 
23 
24 //________________________________________________________________________
27  fPropDist(440),
28  fDoPropagation(kFALSE),
29  fAttemptProp(kTRUE),
30  fAttemptPropMatch(kFALSE),
31  fMaxDistance(0.1),
32  fAttachEmcalParticles(kFALSE),
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  // Constructor.
43 
44  for(Int_t icent=0; icent<8; ++icent) {
45  for(Int_t ipt=0; ipt<9; ++ipt) {
46  for(Int_t ieta=0; ieta<2; ++ieta) {
47  fHistMatchEta[icent][ipt][ieta] = 0;
48  fHistMatchPhi[icent][ipt][ieta] = 0;
49  }
50  }
51  }
52 }
53 
54 //________________________________________________________________________
56  AliAnalysisTaskEmcal(name, histo),
57  fPropDist(440),
58  fDoPropagation(kFALSE),
59  fAttemptProp(kTRUE),
60  fAttemptPropMatch(kFALSE),
61  fMaxDistance(0.1),
62  fAttachEmcalParticles(kFALSE),
63  fUpdateTracks(kTRUE),
64  fUpdateClusters(kTRUE),
65  fEmcalTracks(0),
66  fEmcalClusters(0),
67  fNEmcalTracks(0),
68  fNEmcalClusters(0),
69  fHistMatchEtaAll(0),
70  fHistMatchPhiAll(0)
71 {
72  // Standard constructor.
73 
74  for(Int_t icent=0; icent<8; ++icent) {
75  for(Int_t ipt=0; ipt<9; ++ipt) {
76  for(Int_t ieta=0; ieta<2; ++ieta) {
77  fHistMatchEta[icent][ipt][ieta] = 0;
78  fHistMatchPhi[icent][ipt][ieta] = 0;
79  }
80  }
81  }
82 }
83 
84 //________________________________________________________________________
86 {
87  // Destructor.
88 }
89 
90 //________________________________________________________________________
92 {
93  // Initialize the analysis.
94 
96  if (tracks) {
97  TClass trackClass(tracks->GetClassName());
98  if (!trackClass.InheritsFrom("AliVTrack")) {
99  tracks->SetClassName("AliVTrack"); // enforce only AliVTrack and derived classes
100  }
101  }
102 
104 
106  if (!fInitialized) return;
107 
108  TString emcalTracksName(Form("EmcalTracks_%s", tracks->GetArrayName().Data()));
109  TString emcalClustersName(Form("EmcalClusters_%s", clusters->GetArrayName().Data()));
110 
111  fEmcalTracks = new TClonesArray("AliEmcalParticle");
112  fEmcalTracks->SetName(emcalTracksName);
113  fEmcalClusters = new TClonesArray("AliEmcalParticle");
114  fEmcalClusters->SetName(emcalClustersName);
115 
116  if (fAttachEmcalParticles) {
119  }
120 }
121 
122 //________________________________________________________________________
124 {
125  // Create my user objects.
126 
127  if (!fCreateHisto) return;
128 
130 
131  const Int_t nCentChBins = fNcentBins * 2;
132 
133  fHistMatchEtaAll = new TH1F("fHistMatchEtaAll", "fHistMatchEtaAll", 400, -0.2, 0.2);
134  fHistMatchPhiAll = new TH1F("fHistMatchPhiAll", "fHistMatchPhiAll", 400, -0.2, 0.2);
137 
138  for(Int_t icent=0; icent<nCentChBins; ++icent) {
139  for(Int_t ipt=0; ipt<9; ++ipt) {
140  for(Int_t ieta=0; ieta<2; ++ieta) {
141  TString nameEta(Form("fHistMatchEta_%i_%i_%i",icent,ipt,ieta));
142  fHistMatchEta[icent][ipt][ieta] = new TH1F(nameEta, nameEta, 400, -0.2, 0.2);
143  fHistMatchEta[icent][ipt][ieta]->SetXTitle("#Delta#eta");
144  TString namePhi(Form("fHistMatchPhi_%i_%i_%i",icent,ipt,ieta));
145  fHistMatchPhi[icent][ipt][ieta] = new TH1F(namePhi, namePhi, 400, -0.2, 0.2);
146  fHistMatchPhi[icent][ipt][ieta]->SetXTitle("#Delta#phi");
147  fOutput->Add(fHistMatchEta[icent][ipt][ieta]);
148  fOutput->Add(fHistMatchPhi[icent][ipt][ieta]);
149  }
150  }
151  }
152 
153  PostData(1, fOutput);
154 }
155 
156 //________________________________________________________________________
158 {
159  // Get momenum bin.
160 
161  Int_t pbin=-1;
162  if (p<0.5)
163  pbin=0;
164  else if (p>=0.5 && p<1.0)
165  pbin=1;
166  else if (p>=1.0 && p<1.5)
167  pbin=2;
168  else if (p>=1.5 && p<2.)
169  pbin=3;
170  else if (p>=2. && p<3.)
171  pbin=4;
172  else if (p>=3. && p<4.)
173  pbin=5;
174  else if (p>=4. && p<5.)
175  pbin=6;
176  else if (p>=5. && p<8.)
177  pbin=7;
178  else if (p>=8.)
179  pbin=8;
180 
181  return pbin;
182 }
183 
184 //________________________________________________________________________
186 {
187  // Run the matching.
188 
190  DoMatching();
193 
194  return kTRUE;
195 }
196 
197 //________________________________________________________________________
199 {
200  // Create AliEmcalParticle collections to handle the matching efficiently.
201  // At the same time propagates tracks, if requested.
202 
205 
206  fEmcalTracks->Delete();
207  fEmcalClusters->Delete();
208 
209  fNEmcalTracks = 0;
210  fNEmcalClusters = 0;
211 
212  AliVCluster* cluster = 0;
213  AliVTrack* track = 0;
214 
215  clusters->ResetCurrentID();
216  while ((cluster = static_cast<AliVCluster*>(clusters->GetNextAcceptCluster()))) {
217 
218  // Clears the matching info
219  cluster->SetEmcCpvDistance(-1);
220  cluster->SetTrackDistance(1024, 1024);
221  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
222  AliESDCaloCluster *ec = 0;
223  if (ac) {
224  const Int_t N = ac->GetNTracksMatched();
225  for (Int_t i = N - 1; i >= 0; i--) {
226  TObject *ptr = ac->GetTrackMatched(i);
227  ac->RemoveTrackMatched(ptr);
228  }
229  }
230  else {
231  ec = dynamic_cast<AliESDCaloCluster*>(cluster);
232  TArrayI *arr = ec->GetTracksMatched();
233  if (arr) arr->Set(0);
234  }
235 
236  // Create AliEmcalParticle objects to handle the matching
237  AliEmcalParticle* emcalCluster = new ((*fEmcalClusters)[fNEmcalClusters])
238  AliEmcalParticle(cluster, clusters->GetCurrentID(), fVertex[0], fVertex[1], fVertex[2], AliVCluster::kNonLinCorr);
239  emcalCluster->SetMatchedPtr(fEmcalTracks);
240 
241  fNEmcalClusters++;
242  }
243 
244  tracks->ResetCurrentID();
245  while ((track = static_cast<AliVTrack*>(tracks->GetNextAcceptParticle()))) {
246 
247  // Propagate tracks if requested
248  Bool_t propthistrack = kFALSE;
249  if (fDoPropagation) {
250  propthistrack = kTRUE;
251  }
252  else if (!track->IsExtrapolatedToEMCAL()) {
253  if (fAttemptProp) {
254  propthistrack = kTRUE;
255  }
256  else if (fAttemptPropMatch && IsTrackInEmcalAcceptance(track)) {
257  propthistrack = kTRUE;
258  }
259  }
260  if (propthistrack) AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(track, fPropDist);
261 
262  // Create AliEmcalParticle objects to handle the matching
263  AliEmcalParticle* emcalTrack = new ((*fEmcalTracks)[fNEmcalTracks]) AliEmcalParticle(track, tracks->GetCurrentID());
264  emcalTrack->SetMatchedPtr(fEmcalClusters);
265 
266  AliDebug(2, Form("Now adding track (pT = %.3f, eta = %.3f, phi = %.3f)"
267  "Phi, Eta on EMCal = %.3f, %.3f",
268  emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
269  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal()));
270 
271  fNEmcalTracks++;
272  }
273 }
274 
275 //________________________________________________________________________
277 {
278  // Set the links between tracks and clusters.
279 
280  const Double_t maxd2 = fMaxDistance*fMaxDistance;
281 
282  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
283  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
284  AliVTrack* track = emcalTrack->GetTrack();
285 
286  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
287  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
288  AliVCluster* cluster = emcalCluster->GetCluster();
289 
290  Double_t deta = 999;
291  Double_t dphi = 999;
292  GetEtaPhiDiff(track, cluster, dphi, deta);
293  Double_t d2 = deta * deta + dphi * dphi;
294  if (d2 > maxd2) continue;
295 
296  Double_t d = TMath::Sqrt(d2);
297  emcalCluster->AddMatchedObj(itrack, d);
298  emcalTrack->AddMatchedObj(icluster, d);
299  AliDebug(2, Form("Now matching cluster E = %.3f, pT = %.3f, eta = %.3f, phi = %.3f "
300  "with track pT = %.3f, eta = %.3f, phi = %.3f"
301  "Track eta, phi on EMCal = %.3f, %.3f, d = %.3f",
302  cluster->GetNonLinCorrEnergy(), emcalCluster->Pt(), emcalCluster->Eta(), emcalCluster->Phi(),
303  emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi(),
304  track->GetTrackEtaOnEMCal(), track->GetTrackPhiOnEMCal(), d));
305 
306  if (fCreateHisto) {
307  Int_t mombin = GetMomBin(track->P());
308  Int_t centbinch = fCentBin;
309  if (track->Charge() < 0) centbinch += fNcentBins;
310  Int_t etabin = 0;
311  if(track->Eta() > 0) etabin = 1;
312 
313  fHistMatchEta[centbinch][mombin][etabin]->Fill(deta);
314  fHistMatchPhi[centbinch][mombin][etabin]->Fill(dphi);
315  fHistMatchEtaAll->Fill(deta);
316  fHistMatchPhiAll->Fill(dphi);
317  }
318  }
319  }
320 }
321 
322 //________________________________________________________________________
324 {
325  // Update clusters with matching info.
326 
327  for (Int_t icluster = 0; icluster < fNEmcalClusters; icluster++) {
328  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(icluster));
329  const Int_t N = emcalCluster->GetNumberOfMatchedObj();
330  AliVCluster* cluster = emcalCluster->GetCluster();
331  AliDebug(3, Form("Cluster E = %.2f, eta = %.2f, phi = %.2f, Nmatch = %d", cluster->GetNonLinCorrEnergy(), emcalCluster->Eta(), emcalCluster->Phi(), N));
332 
333  if (N <= 0) continue;
334 
335  // Set the first match distance
336  const UInt_t firstMatchId = emcalCluster->GetMatchedObjId();
337  AliEmcalParticle* emcalTrackFirstMatch = static_cast<AliEmcalParticle*>(fEmcalTracks->At(firstMatchId));
338  AliVTrack* trackFirstMatch = emcalTrackFirstMatch->GetTrack();
339  Double_t deta = 999;
340  Double_t dphi = 999;
341  GetEtaPhiDiff(trackFirstMatch, cluster, dphi, deta);
342  cluster->SetTrackDistance(deta, dphi);
343 
344  // Cast into ESD/AOD objects
345  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
346  AliESDCaloCluster *ec = 0;
347  if (!ac) ec = dynamic_cast<AliESDCaloCluster*>(cluster);
348 
349  // Copy the matched tracks in the cluster. Note: different methods for ESD/AOD
350  if (ac) {
351  for (Int_t i=0; i < N; ++i) {
352  Int_t id = emcalCluster->GetMatchedObjId(i);
353  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
354 
355  AliDebug(3, Form("Pt = %.2f, eta = %.2f, phi = %.2f", emcalTrack->Pt(), emcalTrack->Eta(), emcalTrack->Phi()));
356 
357  TObject *obj = emcalTrack->GetTrack();
358  ac->AddTrackMatched(obj);
359  }
360  }
361  else {
362  TArrayI arr(N);
363  for (Int_t i = 0; i < N; ++i) {
364  Int_t id = emcalCluster->GetMatchedObjId(i);
365  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(id));
366  arr.AddAt(emcalTrack->IdInCollection(), i);
367  }
368  ec->AddTracksMatched(arr);
369  }
370  }
371 }
372 
373 //________________________________________________________________________
375 {
376  // Update tracks with matching info.
377 
378  for (Int_t itrack = 0; itrack < fNEmcalTracks; itrack++) {
379  AliEmcalParticle* emcalTrack = static_cast<AliEmcalParticle*>(fEmcalTracks->At(itrack));
380  AliVTrack* track = emcalTrack->GetTrack();
381 
382  track->ResetStatus(AliVTrack::kEMCALmatch);
383  if (emcalTrack->GetNumberOfMatchedObj() <= 0) continue;
384 
385  AliEmcalParticle* emcalCluster = static_cast<AliEmcalParticle*>(fEmcalClusters->At(emcalTrack->GetMatchedObjId()));
386  track->SetEMCALcluster(emcalCluster->IdInCollection());
387  track->SetStatus(AliVTrack::kEMCALmatch);
388  }
389 }
Int_t GetMatchedObjId(UShort_t i=0) const
Double_t Phi() const
UShort_t GetNumberOfMatchedObj() const
void SetClassName(const char *clname)
void SetMatchedPtr(TObjArray *arr)
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
TClonesArray * fEmcalClusters
emcal tracks
TH1 * fHistMatchPhi[8][9][2]
deta distribution
Int_t fCentBin
event centrality
AliVCluster * GetNextAcceptCluster(Int_t i=-1)
AliVTrack * GetTrack() const
TList * fOutput
x-section from pythia header
const TString & GetClassName() const
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Bool_t IsTrackInEmcalAcceptance(AliVParticle *part, Double_t edges=0.9) const
Double_t Eta() const
Double_t Pt() const
Int_t GetCurrentID() const
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Int_t IdInCollection() const
AliVParticle * GetNextAcceptParticle(Int_t i=-1)
TH1 * fHistMatchEta[8][9][2]
dphi distribution
const TString & GetArrayName() const
ClassImp(AliEmcalClusTrackMatcherTask) AliEmcalClusTrackMatcherTask
Int_t fNEmcalClusters
number of emcal tracks
AliVCluster * GetCluster() const
Double_t fVertex[3]
event plane V0C
void AddMatchedObj(Int_t id, Double_t d)
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
TH1 * fHistMatchEtaAll
number of emcal clusters
void ResetCurrentID(Int_t i=0)