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