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