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