]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALTracker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
avoid redefining common constants, or using magic numbers in code, by instead having...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTracker.cxx
CommitLineData
fe17d4cb 1//========================================================================
2// Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved.
3//
4// Author: The ALICE Off-line Project.
5// Contributors are mentioned in the code where appropriate.
6//
7// Permission to use, copy, modify and distribute this software and its
8// documentation strictly for non-commercial purposes is hereby granted
9// without fee, provided that the above copyright notice appears in all
10// copies and that both the copyright notice and this permission notice
11// appear in the supporting documentation. The authors make no claims
12// about the suitability of this software for any purpose. It is
13// provided "as is" without express or implied warranty.
14//========================================================================
15//
16// Class AliEMCALTracker
17// -----------------------
18// Implementation of the track matching method between barrel tracks and
19// EMCAL clusters.
20// Besides algorithm implementation, some cuts are required to be set
21// in order to define, for each track, an acceptance window where clusters
22// are searched to find best match (if any).
23// The class accepts as input an ESD container, and works directly on it,
24// simply setting, for each of its tracks, the fEMCALindex flag, for each
25// track which is matched to a cluster.
26// In order to use method, one must launch PropagateBack().
27//
28// ------------------------------------------------------------------------
29// author: A. Pulvirenti (alberto.pulvirenti@ct.infn.it)
30//=========================================================================
31
32#include <Riostream.h>
33#include <iomanip>
34
35#include <TFile.h>
36#include <TTree.h>
fe17d4cb 37#include <TList.h>
38#include <TString.h>
39#include <TVector3.h>
40#include <TClonesArray.h>
c61f0e70 41#include <TGeoMatrix.h>
fe17d4cb 42
43#include "AliLog.h"
af885e0f 44#include "AliESDEvent.h"
fe17d4cb 45#include "AliESDtrack.h"
89ffc0b0 46#include "AliESDCaloCluster.h"
fe17d4cb 47#include "AliEMCALRecPoint.h"
48#include "AliRunLoader.h"
49#include "AliEMCALTrack.h"
50#include "AliEMCALLoader.h"
c61f0e70 51#include "AliEMCALGeometry.h"
8ba062b1 52#include "AliEMCALReconstructor.h"
53#include "AliEMCALRecParam.h"
54#include "AliCDBEntry.h"
55#include "AliCDBManager.h"
c61f0e70 56
fe17d4cb 57#include "AliEMCALTracker.h"
58
59ClassImp(AliEMCALTracker)
8ba062b1 60
fe17d4cb 61//
62//------------------------------------------------------------------------------
63//
64AliEMCALTracker::AliEMCALTracker()
65 : AliTracker(),
66 fNPropSteps(0),
67 fTrackCorrMode(kTrackCorrNone),
68 fCutX(50.0),
69 fCutY(50.0),
70 fCutZ(50.0),
71 fCutAlphaMin(-200.0),
72 fCutAlphaMax(200.0),
73 fCutAngle(100.0),
5a14e691 74 fMaxDist(10.0),
dcd86c5d 75 fCutNITS(3.0),
76 fCutNTPC(20.0),
fe17d4cb 77 fRho(1.0),
78 fX0(1.0),
79 fTracks(0),
80 fClusters(0),
c61f0e70 81 fMatches(0),
82 fGeom(0)
fe17d4cb 83{
84 //
85 // Default constructor.
86 // Initializes al simple data members to default values,
87 // and all collections to NULL.
88 // Output file name is set to a default value.
89 //
0832a2bf 90 InitParameters();
fe17d4cb 91}
92//
93//------------------------------------------------------------------------------
94//
95AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
96 : AliTracker(),
97 fNPropSteps(copy.fNPropSteps),
98 fTrackCorrMode(copy.fTrackCorrMode),
99 fCutX(copy.fCutX),
100 fCutY(copy.fCutY),
101 fCutZ(copy.fCutZ),
102 fCutAlphaMin(copy.fCutAlphaMin),
103 fCutAlphaMax(copy.fCutAlphaMax),
104 fCutAngle(copy.fCutAngle),
105 fMaxDist(copy.fMaxDist),
dcd86c5d 106 fCutNITS(copy.fCutNITS),
107 fCutNTPC(copy.fCutNTPC),
fe17d4cb 108 fRho(copy.fRho),
109 fX0(copy.fX0),
110 fTracks((TObjArray*)copy.fTracks->Clone()),
111 fClusters((TObjArray*)copy.fClusters->Clone()),
c61f0e70 112 fMatches((TList*)copy.fMatches->Clone()),
113 fGeom(copy.fGeom)
fe17d4cb 114{
115 //
116 // Copy constructor
117 // Besides copying all parameters, duplicates all collections.
118 //
119}
120//
121//------------------------------------------------------------------------------
122//
123AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
124{
125 //
126 // Assignment operator.
127 // Besides copying all parameters, duplicates all collections.
128 //
129
130 fCutX = copy.fCutX;
131 fCutY = copy.fCutY;
132 fCutZ = copy.fCutZ;
133 fCutAlphaMin = copy.fCutAlphaMin;
134 fCutAlphaMax = copy.fCutAlphaMax;
135 fCutAngle = copy.fCutAngle;
136 fMaxDist = copy.fMaxDist;
dcd86c5d 137 fCutNITS = copy.fCutNITS;
138 fCutNTPC = copy.fCutNTPC;
fe17d4cb 139
140 fTracks = (TObjArray*)copy.fTracks->Clone();
141 fClusters = (TObjArray*)copy.fClusters->Clone();
142 fMatches = (TList*)copy.fMatches->Clone();
143
c61f0e70 144 fGeom = copy.fGeom;
145
fe17d4cb 146 return (*this);
147}
148//
149//------------------------------------------------------------------------------
150//
8ba062b1 151void AliEMCALTracker::InitParameters()
152{
153 //
154 // Retrieve initialization parameters
155 //
156
157 // Check if the instance of AliEMCALRecParam exists,
ba6de5ea 158 const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
159 if(!recParam){
160 AliFatal("Reconstruction parameters for EMCAL not set!");
8ba062b1 161 }
162
ba6de5ea 163 fCutX = recParam->GetTrkCutX();
164 fCutY = recParam->GetTrkCutY();
165 fCutZ = recParam->GetTrkCutZ();
166 fMaxDist = recParam->GetTrkCutR();
167 fCutAngle = recParam->GetTrkCutAngle();
168 fCutAlphaMin = recParam->GetTrkCutAlphaMin();
169 fCutAlphaMax = recParam->GetTrkCutAlphaMax();
dcd86c5d 170 fCutNITS = recParam->GetTrkCutNITS();
171 fCutNTPC = recParam->GetTrkCutNTPC();
172
8ba062b1 173}
174//
175//------------------------------------------------------------------------------
176//
c61f0e70 177TTree* AliEMCALTracker::SearchTrueMatches()
178{
b9560ae2 179 //Search through the list of
180 //track match candidates and clusters
181 //and look for true matches
182 //
183 //
c61f0e70 184 if (!fClusters) return 0;
185 if (fClusters->IsEmpty()) return 0;
186 if (!fTracks) return 0;
187 if (fTracks->IsEmpty()) return 0;
188
189 TTree *outTree = new TTree("tree", "True matches from event");
190 Int_t indexT, indexC, label;
191 outTree->Branch("indexC", &indexC, "indexC/I");
192 outTree->Branch("indexT", &indexT, "indexT/I");
193 outTree->Branch("label", &label , "label/I");
194
195 Double_t dist;
196 Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
197 Int_t it, nTracks = fTracks->GetEntries();
198
199 for (ic = 0; ic < nClusters; ic++) {
200 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
201 label = cluster->Label();
202 indexC = cluster->Index();
203 for (it = 0; it < nTracks; it++) {
204 AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
205 if (TMath::Abs(track->GetSeedLabel()) != label) continue;
206 dist = CheckPair(track, cluster);
207 if (dist <= fMaxDist) {
208 indexT = track->GetSeedIndex();
209 outTree->Fill();
210 }
211 }
212 }
213
214 return outTree;
215}
216//
217//------------------------------------------------------------------------------
218//
fe17d4cb 219void AliEMCALTracker::Clear(Option_t* option)
220{
221 //
222 // Clearing method
044225d9 223 // Deletes all objects in arrays and the arrays themselves
fe17d4cb 224 //
225
226 TString opt(option);
227 Bool_t clearTracks = opt.Contains("TRACKS");
228 Bool_t clearClusters = opt.Contains("CLUSTERS");
229 Bool_t clearMatches = opt.Contains("MATCHES");
230 if (opt.Contains("ALL")) {
231 clearTracks = kTRUE;
232 clearClusters = kTRUE;
233 clearMatches = kTRUE;
234 }
235
236 if (fTracks != 0x0 && clearTracks) {
044225d9 237 fTracks->Delete();
238 delete fTracks;
239 fTracks = 0;
fe17d4cb 240 }
241 if (fClusters != 0x0 && clearClusters) {
044225d9 242 fClusters->Delete();
243 delete fClusters;
244 fClusters = 0;
fe17d4cb 245 }
246 if (fMatches != 0x0 && clearMatches) {
044225d9 247 fMatches->Delete();
248 delete fMatches;
249 fMatches = 0;
fe17d4cb 250 }
251}
252//
253//------------------------------------------------------------------------------
254//
255Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
256{
257 //
258 // Load EMCAL clusters in the form of AliEMCALRecPoint,
259 // from simulation temporary files.
260 // (When included in reconstruction chain, this method is used automatically)
261 //
262
263 Clear("CLUSTERS");
264
265 TBranch *branch = cTree->GetBranch("EMCALECARP");
266 if (!branch) {
c61f0e70 267 AliError("Can't get the branch with the EMCAL clusters");
fe17d4cb 268 return 1;
269 }
270
9596d957 271 TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
fe17d4cb 272 branch->SetAddress(&clusters);
fe17d4cb 273
274 cTree->GetEvent(0);
9596d957 275 Int_t nClusters = (Int_t)clusters->GetEntries();
fe17d4cb 276 fClusters = new TObjArray(0);
277 for (Int_t i = 0; i < nClusters; i++) {
278 AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
279 if (!cluster) continue;
8ada0ffe 280 if (cluster->GetClusterType() != AliESDCaloCluster::kEMCALClusterv1) continue;
fe17d4cb 281 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
282 fClusters->AddLast(matchCluster);
283 }
3d9e8b15 284
2ad4424e 285 branch->SetAddress(0);
286 clusters->Delete();
287 delete clusters;
3d9e8b15 288 if (fClusters->IsEmpty())
289 AliWarning("No clusters collected");
2ad4424e 290
044225d9 291 AliInfo(Form("Collected %d clusters (RecPoints)", fClusters->GetEntries()));
fe17d4cb 292
293 return 0;
294}
295//
296//------------------------------------------------------------------------------
297//
af885e0f 298Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
fe17d4cb 299{
300 //
301 // Load EMCAL clusters in the form of AliESDCaloClusters,
302 // from an AliESD object.
303 //
304
305 // make sure that tracks/clusters collections are empty
306 Clear("CLUSTERS");
307
308 Int_t start = esd->GetFirstEMCALCluster();
309 Int_t nClustersEMC = esd->GetNumberOfEMCALClusters();
310 Int_t end = start + nClustersEMC;
311
312 fClusters = new TObjArray(0);
313
314 Int_t i;
315 for (i = start; i < end; i++) {
316 AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
317 if (!cluster) continue;
8ada0ffe 318 if (cluster->GetClusterType() != AliESDCaloCluster::kEMCALClusterv1) continue;
fe17d4cb 319 AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
320 fClusters->AddLast(matchCluster);
321 }
3d9e8b15 322 if (fClusters->IsEmpty())
323 AliWarning("No clusters collected");
fe17d4cb 324
044225d9 325 AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
fe17d4cb 326
327 return 0;
328}
329//
330//------------------------------------------------------------------------------
331//
af885e0f 332Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
fe17d4cb 333{
334 //
335 // Load ESD tracks.
336 //
337
338 Clear("TRACKS");
339
340 Int_t nTracks = esd->GetNumberOfTracks();
341 fTracks = new TObjArray(0);
342
343 Int_t i, j;
344 Bool_t isKink;
345 Double_t alpha;
346 for (i = 0; i < nTracks; i++) {
347 AliESDtrack *esdTrack = esd->GetTrack(i);
348 // set by default the value corresponding to "no match"
c61f0e70 349 esdTrack->SetEMCALcluster(kUnmatched);
fe17d4cb 350// if (esdTrack->GetLabel() < 0) continue;
351// if (!(esdTrack->GetStatus() & AliESDtrack::kTOFout)) continue;
352 isKink = kFALSE;
353 for (j = 0; j < 3; j++) {
c61f0e70 354 if (esdTrack->GetKinkIndex(j) != 0) isKink = kTRUE;
fe17d4cb 355 }
356 if (isKink) continue;
357 AliEMCALTrack *track = new AliEMCALTrack(*esdTrack);
c61f0e70 358 track->SetMass(0.13957018);
fe17d4cb 359 // check alpha and reject the tracks which fall outside EMCAL acceptance
360 alpha = track->GetAlpha() * TMath::RadToDeg();
361 if (alpha > -155.0 && alpha < 67.0) {
362 delete track;
363 continue;
364 }
c61f0e70 365// if (!PropagateToEMCAL(track)) {
366// delete track;
367// continue;
368// }
fe17d4cb 369 track->SetSeedIndex(i);
370 track->SetSeedLabel(esdTrack->GetLabel());
371 fTracks->AddLast(track);
372 }
373 if (fTracks->IsEmpty()) {
3d9e8b15 374 AliWarning("No tracks collected");
fe17d4cb 375 }
376
377 AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
378
379 return 0;
380}
381//
382//------------------------------------------------------------------------------
383//
af885e0f 384Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
fe17d4cb 385{
386 //
387 // Main operation method.
388 // Gets external AliESD containing tracks to be matched.
389 // After executing match finding, stores in the same ESD object all infos
390 // and releases the object for further reconstruction steps.
391 //
3d9e8b15 392 //
393 // Note: should always return 0=OK, because otherwise all tracking
394 // is aborted for this event
395
fe17d4cb 396 if (!esd) {
397 AliError("NULL ESD passed");
398 return 1;
399 }
400
401 // step 1:
402 // if cluster array is empty, cluster are collected
403 // from the passed ESD, and work is done with ESDCaloClusters
c61f0e70 404 Int_t okLoadClusters, nClusters;
fe17d4cb 405 if (!fClusters || (fClusters && fClusters->IsEmpty())) {
fe17d4cb 406 okLoadClusters = LoadClusters(esd);
fe17d4cb 407 }
c61f0e70 408 nClusters = fClusters->GetEntries();
fe17d4cb 409
410 // step 2:
411 // collect ESD tracks
c61f0e70 412 Int_t okLoadTracks = LoadTracks(esd), nTracks;
fe17d4cb 413 if (okLoadTracks) return 3;
c61f0e70 414 nTracks = fTracks->GetEntries();
fe17d4cb 415
416 // step 3:
417 // each track is propagated to the "R" position of each cluster.
418 // The closest cluster is assigned as match.
419 // IF no clusters lie within the maximum allowed distance, no matches are assigned.
420 Int_t nMatches = CreateMatches();
421 if (!nMatches) {
c61f0e70 422 AliInfo(Form("#clusters = %d -- #tracks = %d --> No good matches found.", nClusters, nTracks));
3d9e8b15 423 return 0;
fe17d4cb 424 }
425 else {
c61f0e70 426 AliInfo(Form("#clusters = %d -- #tracks = %d --> Found %d matches.", nClusters, nTracks, nMatches));
fe17d4cb 427 }
428
429 // step 4:
430 // when more than 1 track share the same matched cluster, only the closest one is kept.
431 Int_t nRemoved = SolveCompetitions();
432 AliInfo(Form("Removed %d duplicate matches", nRemoved));
433 if (nRemoved >= nMatches) {
434 AliError("Removed ALL matches! Check the algorithm or data. Nothing to save");
435 return 5;
436 }
437
438 // step 5:
439 // save obtained information setting the 'fEMCALindex' field of AliESDtrack object
c61f0e70 440 Int_t nSaved = 0, trackID, nGood = 0, nFake = 0;
fe17d4cb 441 TListIter iter(fMatches);
442 AliEMCALMatch *match = 0;
443 while ( (match = (AliEMCALMatch*)iter.Next()) ) {
444 if (!match->CanBeSaved()) continue;
445 AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(match->GetIndexT());
446 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(match->GetIndexC());
447 trackID = track->GetSeedIndex();
448 AliESDtrack *esdTrack = esd->GetTrack(trackID);
449 if (!esdTrack) continue;
dcd86c5d 450
451 // cut on its and tpc track hits
452 if(esdTrack->GetNcls(0)<=fCutNITS)continue;
453 if(esdTrack->GetNcls(1)<=fCutNTPC)continue;
454
c61f0e70 455 if (TMath::Abs(esdTrack->GetLabel()) == cluster->Label()) {
fe17d4cb 456 esdTrack->SetEMCALcluster(cluster->Index());
c61f0e70 457 nGood++;
fe17d4cb 458 }
459 else {
460 esdTrack->SetEMCALcluster(-cluster->Index());
c61f0e70 461 nFake++;
fe17d4cb 462 }
463 nSaved++;
464 }
465 /*
466 AliEMCALTrack *track = 0;
467 TObjArrayIter tracks(fTracks);
468 while ( (track = (AliEMCALTrack*)tracks.Next()) ) {
469 trackID = track->GetSeedIndex();
470 clusterID = track->GetMatchedClusterIndex();
471 AliESDtrack *esdTrack = esd->GetTrack(trackID);
472 if (!esdTrack) continue;
473 if (clusterID < 0) {
c61f0e70 474 esdTrack->SetEMCALcluster(kUnmatched);
fe17d4cb 475 }
476 else {
477 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(clusterID);
478 if (!cluster) continue;
479 if (esdTrack->GetLabel() == cluster->Label()) {
480 nGood++;
481 esdTrack->SetEMCALcluster(cluster->Index());
482 }
483 else {
484 esdTrack->SetEMCALcluster(-cluster->Index());
485 }
486 nSaved++;
487 }
488 }
489 */
c61f0e70 490 AliInfo(Form("Saved %d matches [%d good + %d fake]", nSaved, nGood, nFake));
fe17d4cb 491
492 return 0;
493}
494//
495//------------------------------------------------------------------------------
496//
497void AliEMCALTracker::SetTrackCorrectionMode(Option_t *option)
498{
499 //
500 // Set track correction mode
501 // gest the choice in string format and converts into
502 // internal enum
503 //
504
505 TString opt(option);
506 opt.ToUpper();
507
508 if (!opt.CompareTo("NONE")) {
509 fTrackCorrMode = kTrackCorrNone;
510 }
511 else if (!opt.CompareTo("MMB")) {
512 fTrackCorrMode = kTrackCorrMMB;
513 }
514 else if (!opt.CompareTo("FIXED")) {
515 fTrackCorrMode = kTrackCorrFixed;
516 }
517 else {
518 cerr << "E-AliEMCALTracker::SetTrackCorrectionMode '" << option << "': Unrecognized option" << endl;
519 }
520}
521
522//
523//------------------------------------------------------------------------------
524//
525Double_t AliEMCALTracker::AngleDiff(Double_t angle1, Double_t angle2)
526{
527 //
528 // [PRIVATE]
529 // Given two angles in radiants, it converts them in the range 0-2pi
530 // then computes their true difference, i.e. if the difference a1-a2
531 // results to be larger than 180 degrees, it returns 360 - diff.
532 //
533
534 if (angle1 < 0.0) angle1 += TMath::TwoPi();
535 if (angle1 > TMath::TwoPi()) angle1 -= TMath::TwoPi();
536 if (angle2 < 0.0) angle2 += TMath::TwoPi();
537 if (angle2 > TMath::TwoPi()) angle2 -= TMath::TwoPi();
538
539 Double_t diff = TMath::Abs(angle1 - angle2);
540 if (diff > TMath::Pi()) diff = TMath::TwoPi() - diff;
541
542 if (angle2 > angle1) diff = -diff;
543
544 return diff;
545}
546//
547//------------------------------------------------------------------------------
548//
549Double_t AliEMCALTracker::CheckPair
550(AliEMCALTrack *track, AliEMCALMatchCluster *cl)
551{
552 //
553 // Given a track and a cluster,
554 // propagates the first to the radius of the second.
555 // Then, checks the propagation point against all cuts.
556 // If at least a cut is not passed, a valuer equal to
557 // twice the maximum allowed distance is passed (so the value returned
558 // will not be taken into account when creating matches)
559 //
560
561 // TEMP
562 Bool_t isTrue = kFALSE;
563// if (tr->GetSeedLabel() == cl->Label()) {
564// isTrue = kTRUE;
fe17d4cb 565// }
566
567 // copy track into temporary variable
568 AliEMCALTrack *tr = new AliEMCALTrack(*track);
569
570 Double_t distance = 2.0 * fMaxDist;
571
572 // check against cut on difference 'alpha - phi'
573 Double_t phi = TMath::ATan2(cl->Y(), cl->X());
574 phi = AngleDiff(phi, tr->GetAlpha());
5a14e691 575 if (phi < fCutAlphaMin || phi > fCutAlphaMax){
576 delete tr;
577 return distance;
578 }
fe17d4cb 579
580 // try to propagate to cluster radius
581 // (return the 'distance' value if it fails)
582 Double_t pos[3], &x = pos[0], &y = pos[1], &z = pos[2];
583 Double_t x0, rho;
584 tr->GetXYZ(pos);
585 Double_t rt = TMath::Sqrt(x*x + y*y);
586 Double_t rc = TMath::Sqrt(cl->X()*cl->X() + cl->Y()*cl->Y());
587
588 if (fTrackCorrMode == kTrackCorrMMB) {
589 Double_t pos1[3], pos2[3], param[6];
590 pos1[0] = x;
591 pos1[1] = y;
592 pos1[2] = z;
593 pos2[0] = cl->X();
594 pos2[1] = cl->Y();
595 pos2[2] = cl->Z();
c3c15772 596 MeanMaterialBudget(pos1, pos2, param);
597 rho = param[0]*param[4];
fe17d4cb 598 x0 = param[1];
599 }
600 else if (fTrackCorrMode == kTrackCorrFixed) {
601 rho = fRho;
602 x0 = fX0;
603 }
604 else {
605 rho = 0.0;
606 x0 = 0.0;
607 }
608 if (fNPropSteps) {
609 Int_t i;
610 Double_t r;
611 cout.setf(ios::fixed);
612 cout.precision(5);
613 if (isTrue) cout << "Init : " << rt << ' ' << x << ' ' << y << ' ' << z << endl;
614 for (i = 0; i < fNPropSteps; i++) {
615 r = rt + (rc - rt) * ((Double_t)(i+1)/(Double_t)fNPropSteps);
5a14e691 616 if (!tr->PropagateTo(r, x0, rho)){
617 delete tr;
618 return distance;
619 }
fe17d4cb 620 tr->GetXYZ(pos);
621 if (isTrue) cout << "Step : " << r << ' ' << x << ' ' << y << ' ' << z << endl;
622 }
623 if (isTrue) cout << "Clstr: " << rc << ' ' << cl->X() << ' ' << cl->Y() << ' ' << cl->Z() << endl;
624 }
625 else {
626 // when no steps are used, no correction makes sense
c61f0e70 627 //if (!tr->PropagateTo(rc, 0.0, 0.0)) return distance;
5a14e691 628 if (!tr->PropagateToGlobal(cl->X(), cl->Y(), cl->Z(), 0.0, 0.0)){
629 delete tr;
630 return distance;
631 }
c61f0e70 632 /*
633 Bool_t propOK = kFALSE;
634 cout << "START" << endl;
635 Double_t dist, rCHK, bestDist = 10000000.0;
636 for (Double_t rTMP = rc; rTMP> rc*0.95; rTMP -= 0.1) {
637 if (!tr->PropagateTo(rTMP)) continue;
638 propOK = kTRUE;
639 tr->GetXYZ(pos);
640 rCHK = TMath::Sqrt(x*x + y*y);
641 dist = TMath::Abs(rCHK - rc);
642 cout << rCHK << " vs. " << rc << endl;
643
644 if (TMath::Abs(rCHK - rc) < 0.01) break;
645 }
646 cout << "STOP" << endl;
647 if (!propOK) return distance;
648 */
fe17d4cb 649 }
650
651 // get global propagation of track at end of propagation
652 tr->GetXYZ(pos);
653
654 // check angle cut
655 TVector3 vc(cl->X(), cl->Y(), cl->Z());
656 TVector3 vt(x, y, z);
657 Double_t angle = TMath::Abs(vc.Angle(vt)) * TMath::RadToDeg();
c61f0e70 658 // check: where is the track?
659 Double_t r, phiT, phiC;
660 r = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
661 phiT = TMath::ATan2(pos[1], pos[0]) * TMath::RadToDeg();
662 phiC = vc.Phi() * TMath::RadToDeg();
663 //cout << "Propagated R, phiT, phiC = " << r << ' ' << phiT << ' ' << phiC << endl;
664
fe17d4cb 665 if (angle > fCutAngle) {
666 //cout << "angle" << endl;
5a14e691 667 delete tr;
fe17d4cb 668 return distance;
669 }
670
671 // compute differences wr to each coordinate
672 x -= cl->X();
c61f0e70 673 if (TMath::Abs(x) > fCutX) {
fe17d4cb 674 //cout << "cut X" << endl;
5a14e691 675 delete tr;
fe17d4cb 676 return distance;
677 }
678 y -= cl->Y();
c61f0e70 679 if (TMath::Abs(y) > fCutY) {
fe17d4cb 680 //cout << "cut Y" << endl;
5a14e691 681 delete tr;
fe17d4cb 682 return distance;
683 }
684 z -= cl->Z();
c61f0e70 685 if (TMath::Abs(z) > fCutZ) {
fe17d4cb 686 //cout << "cut Z" << endl;
5a14e691 687 delete tr;
fe17d4cb 688 return distance;
689 }
690
691 // compute true distance
692 distance = TMath::Sqrt(x*x + y*y + z*z);
693 //Double_t temp = CheckPairV2(tr, cl);
694 //if (temp < distance) return temp; else
695
696 // delete temporary object
697 delete tr;
698
699 return distance;
700}
701//
702//------------------------------------------------------------------------------
703//
704Double_t AliEMCALTracker::CheckPairV2
705(AliEMCALTrack *tr, AliEMCALMatchCluster *cl)
706{
707 //
708 // Given a track and a cluster,
709 // propagates the first to the radius of the second.
710 // Then, checks the propagation point against all cuts.
711 // If at least a cut is not passed, a valuer equal to
712 // twice the maximum allowed distance is passed (so the value returned
713 // will not be taken into account when creating matches)
714 //
715
716 // TEMP
717// Bool_t isTrue = kFALSE;
718// if (tr->GetSeedLabel() == cl->Label()) {
719// isTrue = kTRUE;
720// cout << "TRUE MATCH!!!" << endl;
721// }
722
723 Double_t distance = 2.0 * fMaxDist;
724
725 Double_t x0, rho;
726 if (fTrackCorrMode == kTrackCorrMMB) {
727 Double_t pos1[3], pos2[3], param[6];
728 tr->GetXYZ(pos1);
729// pos1[0] = x;
730// pos1[1] = y;
731// pos1[2] = z;
732 pos2[0] = cl->X();
733 pos2[1] = cl->Y();
734 pos2[2] = cl->Z();
c3c15772 735 MeanMaterialBudget(pos1, pos2, param);
736 rho = param[0]*param[4];
fe17d4cb 737 x0 = param[1];
738 }
739 else if (fTrackCorrMode == kTrackCorrFixed) {
740 rho = fRho;
741 x0 = fX0;
742 }
743 else {
744 rho = 0.0;
745 x0 = 0.0;
746 }
747
748 // check against cut on difference 'alpha - phi'
749 Double_t phi = TMath::ATan2(cl->Y(), cl->X());
750 phi = AngleDiff(phi, tr->GetAlpha());
751 if (phi < fCutAlphaMin || phi > fCutAlphaMax) return distance;
752
753 // get cluster position and put them into a vector
754 TVector3 vc(cl->X(), cl->Y(), cl->Z());
755 // rotate the vector in order to put all clusters on a plane intersecting
756 // vertically the X axis; the angle depends on the sector
757 Double_t clusterRot, clusterPhi = vc.Phi() * TMath::RadToDeg();
758 if (clusterPhi < 0.0) clusterPhi += 360.0;
759 if (clusterPhi < 100.0) {
760 clusterRot = -90.0;
761 }
762 else if (clusterPhi < 120.0) {
763 clusterRot = -110.0;
764 }
765 else if (clusterPhi < 140.0) {
766 clusterRot = -130.0;
767 }
768 else if (clusterPhi < 160.0) {
769 clusterRot = -150.0;
770 }
771 else if (clusterPhi < 180.0) {
772 clusterRot = -170.0;
773 }
774 else {
775 clusterRot = -190.0;
776 }
777 vc.RotateZ(clusterRot * TMath::DegToRad());
778 // generate a track from the ESD track selected
779 AliEMCALTrack *track = new AliEMCALTrack(*tr);
780 // compute the 'phi' coordinate of the intersection point to
781 // the EMCAL surface
782 Double_t x = vc.X();
783 Double_t y;
784 track->GetYAt(vc.X(), track->GetBz(), y);
785 Double_t tmp = x*TMath::Cos(track->GetAlpha()) - y*TMath::Sin(track->GetAlpha());
786 y = x*TMath::Sin(track->GetAlpha()) + y*TMath::Cos(track->GetAlpha());
787 x = tmp;
788 Double_t trackPhi = TMath::ATan2(y, x) * TMath::RadToDeg();
789 // compute phi difference
790 Double_t dphi = trackPhi - clusterPhi;
791 if (TMath::Abs(dphi) > 180.0) {
792 dphi = 360.0 - TMath::Abs(dphi);
793 if (clusterPhi > trackPhi) dphi = -dphi;
794 }
795 // propagate track to the X position of rotated cluster
796 // and get the vector of X, Y, Z in the local ref. frame of the track
797 track->PropagateTo(vc.X(), x0, rho);
798 TVector3 vt(track->GetX(), track->GetY(), track->GetZ());
799 vt.RotateZ((clusterPhi - trackPhi) * TMath::DegToRad());
800 TVector3 vdiff = vt-vc;
801
802 // compute differences wr to each coordinate
5a14e691 803 delete track;
fe17d4cb 804 if (vdiff.X() > fCutX) return distance;
805 if (vdiff.Y() > fCutY) return distance;
806 if (vdiff.Z() > fCutZ) return distance;
807
808 // compute true distance
809 distance = vdiff.Mag();
810 return distance;
811}
812//
813//------------------------------------------------------------------------------
814//
c61f0e70 815Double_t AliEMCALTracker::CheckPairV3
816(AliEMCALTrack *track, AliEMCALMatchCluster *cl)
817{
818 //
819 // Given a track and a cluster,
820 // propagates the first to the radius of the second.
821 // Then, checks the propagation point against all cuts.
822 // If at least a cut is not passed, a valuer equal to
823 // twice the maximum allowed distance is passed (so the value returned
824 // will not be taken into account when creating matches)
825 //
826
827 AliEMCALTrack tr(*track);
828
829 Int_t sector;
830 Double_t distance = 2.0 * fMaxDist;
831 Double_t dx, dy, dz;
832 Double_t phi, alpha, slope, tgtXnum, tgtXden, sectorWidth = 20.0 * TMath::DegToRad();
e77c61fe 833 Double_t xcurr, xprop, param[6] = {0., 0., 0., 0., 0., 0.}, x0, rho, bz;
c61f0e70 834 Double_t x[3], x1[3], x2[3];
835
836 // get initial track position
837 xcurr = tr.GetX();
838
839 // evaluate the EMCAL sector number
840 phi = cl->Phi();
841 if (phi < 0.0) phi += TMath::TwoPi();
842 sector = (Int_t)(phi / sectorWidth);
843 alpha = ((Double_t)sector + 0.5) * sectorWidth;
844 // evaluate the corresponding X for track propagation
845 slope = TMath::Tan(alpha - 0.5*TMath::Pi());
846 tgtXnum = cl->Y() - slope * cl->X();
847 tgtXden = TMath::Sqrt(1.0 + slope*slope);
848 xprop = TMath::Abs(tgtXnum / tgtXden);
849
850 // propagate by small steps
851 tr.GetXYZ(x1);
852 bz = tr.GetBz();
853 if (!tr.GetXYZAt(xprop, bz, x2)) return distance;
854 //AliKalmanTrack::MeanMaterialBudget(x1, x2, param);
e77c61fe 855 rho = param[0]*param[4];
c61f0e70 856 x0 = param[1];
e77c61fe 857 if (!tr.PropagateTo(xprop, x0, rho)) return distance;
c61f0e70 858 //if (!tr.PropagateTo(xprop, 0.0, 0.0)) return distance;
859
860 // get propagated position at the end
861 tr.GetXYZ(x);
862 dx = TMath::Abs(x[0] - cl->X());
863 dy = TMath::Abs(x[1] - cl->Y());
864 dz = TMath::Abs(x[2] - cl->Z());
865 if (dx > fCutX || dy > fCutY || dz > fCutZ) return distance;
866
867 distance = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
868
869 return distance;
870}
871//
872//------------------------------------------------------------------------------
873//
874Bool_t AliEMCALTracker::PropagateToEMCAL(AliEMCALTrack *tr)
875{
876 //
877 // Propagates the track to the proximity of the EMCAL surface
878 //
879
e77c61fe 880 Double_t xcurr, xtemp, xprop = 438.0, step = 10.0, param[6], x0, rho, bz;
c61f0e70 881 Double_t x1[3], x2[3];
882
883 // get initial track position
884 xcurr = tr->GetX();
885
886 // propagate by small steps
887 for (xtemp = xcurr + step; xtemp < xprop; xtemp += step) {
888 // to compute material budget, take current position and
889 // propagated hypothesis without energy loss
890 tr->GetXYZ(x1);
891 bz = tr->GetBz();
892 if (!tr->GetXYZAt(xtemp, bz, x2)) return kFALSE;
c3c15772 893 MeanMaterialBudget(x1, x2, param);
894 rho = param[0]*param[4];
c61f0e70 895 x0 = param[1];
e77c61fe 896 if (!tr->PropagateTo(xtemp, x0, rho)) return kFALSE;
c61f0e70 897 }
898
899 return kTRUE;
900}
901//
902//------------------------------------------------------------------------------
903//
fe17d4cb 904Int_t AliEMCALTracker::CreateMatches()
905{
906 //
907 // Creation of matches between tracks and clusters.
908 // For each ESD track collected by ReadESD(), an AliEMCALTrack is made.
909 // If it finds a cluster close enough to its propagation to EMCAL,
910 // which passes all cuts, its index is stored.
911 // If many clusters are found which satisfy the criteria described above,
912 // only the closest one is stored.
913 // At this level, it is possible that two tracks share the same cluster.
914 //
915
916 // if matches collection is already present, it is deleted
917 if (fMatches) {
918 fMatches->Delete();
919 delete fMatches;
920 }
921 fMatches = new TList;
922
923 // initialize counters and indexes
924 Int_t count = 0;
925 Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
926 Int_t it, nTracks = fTracks->GetEntries();
927
928 // external loop on clusters, internal loop on tracks
929 Double_t dist;
930 for (ic = 0; ic < nClusters; ic++) {
fe17d4cb 931 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
932 for (it = 0; it < nTracks; it++) {
933 AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
934 dist = CheckPair(track, cluster);
c61f0e70 935 //cout << dist << endl;
fe17d4cb 936 if (dist <= fMaxDist) {
937 AliEMCALMatch *candidate = new AliEMCALMatch;
938 candidate->SetIndexT(it);
939 candidate->SetIndexC(ic);
940 candidate->SetDistance(dist);
6c23ffed 941 candidate->SetPt(track->GetSignedPt());
fe17d4cb 942 fMatches->Add(candidate);
943 count++;
944 }
945 }
946 }
c61f0e70 947
fe17d4cb 948 /*
949 // loop on clusters and tracks
950 Int_t icBest;
951 Double_t dist, distBest;
952 for (it = 0; it < nTracks; it++) {
953 AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
954 if (!track) continue;
955 icBest = -1;
956 distBest = fMaxDist;
957 for (ic = 0; ic < nClusters; ic++) {
958 AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
959 if (!cluster) continue;
960 dist = CheckPair(track, cluster);
961 if (dist < distBest) {
962 distBest = dist;
963 icBest = ic;
964 }
965 }
966 if (icBest >= 0) {
967 track->SetMatchedClusterIndex(icBest);
968 track->SetMatchedClusterDist(distBest);
969 count++;
970 }
971 else {
972 track->SetMatchedClusterIndex(-1);
973 }
974 }
975 */
976
977 return count;
978}
979//
980//------------------------------------------------------------------------------
981//
982Int_t AliEMCALTracker::SolveCompetitions()
983{
984 //
985 // Match selector.
986 // The match list is sorted from the best to the worst match, w.r. to the
987 // distance between track prolongation and cluster position.
988 // Based on this criterion, starting from the first (best) match, a flag
989 // is set to both the involved track and cluster, and all matches containing
990 // an already used track or cluster are removed, leaving only the best match
991 // for each cluster.
992 //
993
994 // sort matches with respect to track-cluster distance
995 fMatches->Sort(kSortAscending);
996
997 // keep track of eliminated matches
998 Int_t count = 0;
999
1000 // initialize flags to check repetitions
1001 Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
1002 Int_t it, nTracks = fTracks->GetEntries();
1003 Bool_t *usedC = new Bool_t[nClusters];
1004 Bool_t *usedT = new Bool_t[nTracks];
1005 for (ic = 0; ic < nClusters; ic++) usedC[ic] = kFALSE;
1006 for (it = 0; it < nTracks; it++) usedT[it] = kFALSE;
1007
1008 // loop on matches
1009 TListIter iter(fMatches);
1010 AliEMCALMatch *match = 0;
1011 while ( (match = (AliEMCALMatch*)iter.Next()) ) {
1012 ic = match->GetIndexC();
1013 it = match->GetIndexT();
1014 if (!usedT[it] && !usedC[ic]) {
1015 usedT[it] = kTRUE;
1016 usedC[ic] = kTRUE;
1017 match->CanBeSaved() = kTRUE;
1018 }
1019 else {
1020 count++;
1021 }
1022 }
1023
1024 /*
1025 Int_t it1, it2, nTracks = (Int_t)fTracks->GetEntries();
1026 AliEMCALTrack *track1 = 0, *track2 = 0;
1027 for (it1 = 0; it1 < nTracks; it1++) {
1028 track1 = (AliEMCALTrack*)fTracks->At(it1);
1029 if (!track1) continue;
1030 if (track1->GetMatchedClusterIndex() < 0) continue;
1031 for (it2 = it1+1; it2 < nTracks; it2++) {
1032 track2 = (AliEMCALTrack*)fTracks->At(it2);
1033 if (!track2) continue;
1034 if (track2->GetMatchedClusterIndex() < 0) continue;
1035 if (track1->GetMatchedClusterIndex() != track2->GetMatchedClusterIndex()) continue;
1036 count++;
1037 if (track1->GetMatchedClusterDist() < track2->GetMatchedClusterDist()) {
1038 track2->SetMatchedClusterIndex(-1);
1039 }
1040 else if (track2->GetMatchedClusterDist() < track1->GetMatchedClusterDist()) {
1041 track1->SetMatchedClusterIndex(-1);
1042 }
1043 }
1044 }
1045 */
1046
17773e2e 1047 delete [] usedC;
1048 delete [] usedT;
1049
fe17d4cb 1050 return count;
1051}
1052//
1053//------------------------------------------------------------------------------
1054//
1055void AliEMCALTracker::UnloadClusters()
1056{
1057 //
044225d9 1058 // Free memory from all arrays
1059 // This method is called after the local tracking step
1060 // so we can safely delete everything
fe17d4cb 1061 //
1062
044225d9 1063 Clear();
fe17d4cb 1064}
1065//
1066//------------------------------------------------------------------------------
1067//
1068AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
1069 : fIndex(index),
0f7fcbef 1070 fLabel(recPoint->GetPrimaryIndex()), //wrong! fixed below
fe17d4cb 1071 fX(0.),
1072 fY(0.),
1073 fZ(0.)
1074{
1075 //
1076 // Translates an AliEMCALRecPoint object into the internal format.
1077 // Index of passed cluster in its native array must be specified.
1078 //
1079 TVector3 clpos;
1080 recPoint->GetGlobalPosition(clpos);
1081
1082 fX = clpos.X();
1083 fY = clpos.Y();
1084 fZ = clpos.Z();
0f7fcbef 1085
1086 //AliEMCALRecPoint stores the track labels in the parents
1087 //list, sorted according to the fractional contribution to the
1088 //RecPoint. The zeroth parent gave the highest contribution
1089 Int_t multparent = 0;
1090 Int_t *parents = recPoint->GetParents(multparent);
1091 if(multparent > 0)
1092 fLabel = parents[0];
1093 else
1094 fLabel = -1;
1095
fe17d4cb 1096}
1097//
1098//------------------------------------------------------------------------------
1099//
1100AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDCaloCluster *caloCluster)
1101 : fIndex(index),
7592dfc4 1102 fLabel(caloCluster->GetLabel()),
fe17d4cb 1103 fX(0.),
1104 fY(0.),
1105 fZ(0.)
1106{
1107 //
1108 // Translates an AliESDCaloCluster object into the internal format.
1109 // Index of passed cluster in its native array must be specified.
1110 //
1111 Float_t clpos[3];
7592dfc4 1112 caloCluster->GetPosition(clpos);
fe17d4cb 1113
1114 fX = (Double_t)clpos[0];
1115 fY = (Double_t)clpos[1];
1116 fZ = (Double_t)clpos[2];
1117}
1118//
1119//------------------------------------------------------------------------------
1120//
1121Int_t AliEMCALTracker::AliEMCALMatch::Compare(const TObject *obj) const
1122{
1123 //
1124 // Tracks compared wrt their distance from matched point
1125 //
1126
1127 AliEMCALTracker::AliEMCALMatch *that = (AliEMCALTracker::AliEMCALMatch*)obj;
1128
c61f0e70 1129 Double_t thisDist = fPt;//fDistance;
1130 Double_t thatDist = that->fPt;//that->GetDistance();
fe17d4cb 1131
1132 if (thisDist > thatDist) return 1;
1133 else if (thisDist < thatDist) return -1;
1134 return 0;
1135}
1136
1137AliEMCALTracker::AliEMCALMatch::AliEMCALMatch()
1138 : TObject(),
1139 fCanBeSaved(kFALSE),
1140 fIndexC(0),
1141 fIndexT(0),
c61f0e70 1142 fDistance(0.),
1143 fPt(0.)
fe17d4cb 1144{
1145 //default constructor
1146
1147}
1148
1149AliEMCALTracker::AliEMCALMatch::AliEMCALMatch(const AliEMCALMatch& copy)
1150 : TObject(),
1151 fCanBeSaved(copy.fCanBeSaved),
1152 fIndexC(copy.fIndexC),
1153 fIndexT(copy.fIndexT),
c61f0e70 1154 fDistance(copy.fDistance),
1155 fPt(copy.fPt)
fe17d4cb 1156{
1157 //copy ctor
1158}