]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - EMCAL/AliEMCALTracker.cxx
Next round of fixes (Jochen)
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTracker.cxx
index f3ef163b898107a2ba719c80f6cc3e9aa74c2f1f..c381ff377d464098ad9e12756e0f583a1f331210 100644 (file)
 #include "AliEMCALTrack.h"
 #include "AliEMCALLoader.h"
 #include "AliEMCALGeometry.h"
+#include "AliEMCALReconstructor.h"
+#include "AliEMCALRecParam.h"
+#include "AliCDBEntry.h"
+#include "AliCDBManager.h"
+#include "AliEMCALReconstructor.h"
 
 #include "AliEMCALTracker.h"
 
 ClassImp(AliEMCALTracker)
+
 //
 //------------------------------------------------------------------------------
 //
@@ -66,7 +72,9 @@ AliEMCALTracker::AliEMCALTracker()
     fCutAlphaMin(-200.0),
     fCutAlphaMax(200.0),
     fCutAngle(100.0),
-    fMaxDist(100.0),
+    fMaxDist(10.0),
+    fCutNITS(3.0),
+    fCutNTPC(20.0),
     fRho(1.0),
     fX0(1.0),
     fTracks(0),
@@ -80,6 +88,7 @@ AliEMCALTracker::AliEMCALTracker()
        // and all collections to NULL.
        // Output file name is set to a default value.
        //
+  InitParameters();
 }
 //
 //------------------------------------------------------------------------------
@@ -95,6 +104,8 @@ AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
     fCutAlphaMax(copy.fCutAlphaMax),
     fCutAngle(copy.fCutAngle),
     fMaxDist(copy.fMaxDist),
+    fCutNITS(copy.fCutNITS),
+    fCutNTPC(copy.fCutNTPC),
     fRho(copy.fRho),
     fX0(copy.fX0),
     fTracks((TObjArray*)copy.fTracks->Clone()),
@@ -124,6 +135,8 @@ AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
        fCutAlphaMax = copy.fCutAlphaMax;
        fCutAngle = copy.fCutAngle;
        fMaxDist = copy.fMaxDist;
+       fCutNITS = copy.fCutNITS;
+       fCutNTPC = copy.fCutNTPC;
        
        fTracks = (TObjArray*)copy.fTracks->Clone();
        fClusters = (TObjArray*)copy.fClusters->Clone();
@@ -136,8 +149,41 @@ AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
 //
 //------------------------------------------------------------------------------
 //
+void AliEMCALTracker::InitParameters()
+{
+       //
+       // Retrieve initialization parameters
+       //
+       
+  // Check if the instance of AliEMCALRecParam exists, 
+  const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
+
+  if(!recParam){
+    AliFatal("Reconstruction parameters for EMCAL not set!");
+  }
+  else{
+    fCutX =  recParam->GetTrkCutX();
+    fCutY =  recParam->GetTrkCutY();
+    fCutZ =  recParam->GetTrkCutZ();
+    fMaxDist =  recParam->GetTrkCutR();
+    fCutAngle =  recParam->GetTrkCutAngle();
+    fCutAlphaMin =  recParam->GetTrkCutAlphaMin();
+    fCutAlphaMax =  recParam->GetTrkCutAlphaMax();
+    fCutNITS = recParam->GetTrkCutNITS();
+    fCutNTPC = recParam->GetTrkCutNTPC();
+  }
+       
+}
+//
+//------------------------------------------------------------------------------
+//
 TTree* AliEMCALTracker::SearchTrueMatches()
 {
+  //Search through the list of
+  //track match candidates and clusters
+  //and look for true matches
+  //
+  //
        if (!fClusters) return 0;
        if (fClusters->IsEmpty()) return 0;
        if (!fTracks) return 0;
@@ -149,7 +195,7 @@ TTree* AliEMCALTracker::SearchTrueMatches()
        outTree->Branch("indexT", &indexT, "indexT/I");
        outTree->Branch("label",  &label , "label/I");
        
-       Double_t dist;
+       Double_t dist=0.;
        Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
        Int_t it, nTracks = fTracks->GetEntries();
        
@@ -219,6 +265,9 @@ Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
        
        Clear("CLUSTERS");
 
+       cTree->SetBranchStatus("*",0); //disable all branches
+       cTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
+
        TBranch *branch = cTree->GetBranch("EMCALECARP");
        if (!branch) {
                AliError("Can't get the branch with the EMCAL clusters");
@@ -228,13 +277,15 @@ Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
        TClonesArray *clusters = new TClonesArray("AliEMCALRecPoint", 1000);
        branch->SetAddress(&clusters);
        
-       cTree->GetEvent(0);
+       //cTree->GetEvent(0);
+       branch->GetEntry(0);
        Int_t nClusters = (Int_t)clusters->GetEntries();
-       fClusters = new TObjArray(0);
+       if(fClusters) fClusters->Delete();
+       else fClusters = new TObjArray(0);
        for (Int_t i = 0; i < nClusters; i++) {
                AliEMCALRecPoint *cluster = (AliEMCALRecPoint*)clusters->At(i);
                if (!cluster) continue;
-               if (cluster->GetClusterType() != AliESDCaloCluster::kClusterv1) continue;
+               if (cluster->GetClusterType() != AliVCluster::kEMCALClusterv1) continue;
                AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
                fClusters->AddLast(matchCluster);
        }
@@ -243,9 +294,9 @@ Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
         clusters->Delete();
         delete clusters;
         if (fClusters->IsEmpty())
-           AliWarning("No clusters collected");
+           AliDebug(1,"No clusters collected");
 
-       AliInfo(Form("Collected %d clusters (RecPoints)", fClusters->GetEntries()));
+       AliDebug(1,Form("Collected %d clusters (RecPoints)", fClusters->GetEntries()));
 
        return 0;
 }
@@ -262,8 +313,8 @@ Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
        // make sure that tracks/clusters collections are empty
        Clear("CLUSTERS");
        
-       Int_t start = esd->GetFirstEMCALCluster();
-       Int_t nClustersEMC = esd->GetNumberOfEMCALClusters();
+       Int_t start = 0;
+       Int_t nClustersEMC = esd->GetNumberOfCaloClusters();
        Int_t end = start + nClustersEMC;
        
        fClusters = new TObjArray(0);
@@ -272,14 +323,14 @@ Int_t AliEMCALTracker::LoadClusters(AliESDEvent *esd)
        for (i = start; i < end; i++) {
                AliESDCaloCluster *cluster = esd->GetCaloCluster(i);
                if (!cluster) continue;
-               if (cluster->GetClusterType() != AliESDCaloCluster::kClusterv1) continue;
+        if (!cluster->IsEMCAL()) continue ; 
                AliEMCALMatchCluster *matchCluster = new AliEMCALMatchCluster(i, cluster);
                fClusters->AddLast(matchCluster);
        }
         if (fClusters->IsEmpty())
-           AliWarning("No clusters collected");
+           AliDebug(1,"No clusters collected");
        
-       AliInfo(Form("Collected %d clusters from ESD", fClusters->GetEntries()));
+       AliDebug(1,Form("Collected %d clusters from ESD", fClusters->GetEntries()));
 
        return 0;
 }
@@ -297,9 +348,9 @@ Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
        Int_t nTracks = esd->GetNumberOfTracks();
        fTracks = new TObjArray(0);
        
-       Int_t i, j;
-       Bool_t isKink;
-       Double_t alpha; 
+       Int_t i=0, j=0;
+       Bool_t isKink=kFALSE;
+       Double_t alpha=0.
        for (i = 0; i < nTracks; i++) {
                AliESDtrack *esdTrack = esd->GetTrack(i);
                // set by default the value corresponding to "no match"
@@ -328,10 +379,10 @@ Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
                fTracks->AddLast(track);
        }
        if (fTracks->IsEmpty()) {
-               AliWarning("No tracks collected");
+               AliDebug(1,"No tracks collected");
        }
        
-       AliInfo(Form("Collected %d tracks", fTracks->GetEntries()));
+       AliDebug(1,Form("Collected %d tracks", fTracks->GetEntries()));
 
        return 0;
 }
@@ -376,17 +427,17 @@ Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
        // IF no clusters lie within the maximum allowed distance, no matches are assigned.
        Int_t nMatches = CreateMatches();
        if (!nMatches) {
-               AliInfo(Form("#clusters = %d -- #tracks = %d --> No good matches found.", nClusters, nTracks));
+               AliDebug(1,Form("#clusters = %d -- #tracks = %d --> No good matches found.", nClusters, nTracks));
                return 0;
        }
        else {
-               AliInfo(Form("#clusters = %d -- #tracks = %d --> Found %d matches.", nClusters, nTracks, nMatches));
+               AliDebug(1,Form("#clusters = %d -- #tracks = %d --> Found %d matches.", nClusters, nTracks, nMatches));
        }
        
        // step 4:
        // when more than 1 track share the same matched cluster, only the closest one is kept.
        Int_t nRemoved = SolveCompetitions();
-       AliInfo(Form("Removed %d duplicate matches", nRemoved));
+       AliDebug(1,Form("Removed %d duplicate matches", nRemoved));
        if (nRemoved >= nMatches) {
                AliError("Removed ALL matches! Check the algorithm or data. Nothing to save");
                return 5;
@@ -394,7 +445,7 @@ Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
        
        // step 5:
        // save obtained information setting the 'fEMCALindex' field of AliESDtrack object
-       Int_t nSaved = 0, trackID, nGood = 0, nFake = 0;
+       Int_t nSaved = 0, trackID;
        TListIter iter(fMatches);
        AliEMCALMatch *match = 0;
        while ( (match = (AliEMCALMatch*)iter.Next()) ) {
@@ -404,14 +455,12 @@ Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
                trackID = track->GetSeedIndex();
                AliESDtrack *esdTrack = esd->GetTrack(trackID);
                if (!esdTrack) continue;
-               if (TMath::Abs(esdTrack->GetLabel()) == cluster->Label()) {
-                       esdTrack->SetEMCALcluster(cluster->Index());
-                       nGood++;
-               }
-               else {
-                       esdTrack->SetEMCALcluster(-cluster->Index());
-                       nFake++;
-               }
+
+               // cut on its and tpc track hits
+               if(esdTrack->GetNcls(0)<=fCutNITS)continue;
+               if(esdTrack->GetNcls(1)<=fCutNTPC)continue;
+             
+               esdTrack->SetEMCALcluster(cluster->Index());
                nSaved++;
        }
        /*
@@ -428,18 +477,13 @@ Int_t AliEMCALTracker::PropagateBack(AliESDEvent* esd)
                else {
                        AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(clusterID);
                        if (!cluster) continue;
-                       if (esdTrack->GetLabel() == cluster->Label()) {
-                               nGood++;
-                               esdTrack->SetEMCALcluster(cluster->Index());
-                       }
-                       else {
-                               esdTrack->SetEMCALcluster(-cluster->Index());
-                       }
+               
+                       esdTrack->SetEMCALcluster(cluster->Index());
                        nSaved++;
                }
        }
        */
-       AliInfo(Form("Saved %d matches [%d good + %d fake]", nSaved, nGood, nFake));
+       AliDebug(1,Form("Saved %d matches", nSaved));
 
        return 0;
 }
@@ -511,7 +555,7 @@ Double_t AliEMCALTracker::CheckPair
        //
        
        // TEMP
-       Bool_t isTrue = kFALSE;
+       //Bool_t isTrue = kFALSE;
 //     if (tr->GetSeedLabel() == cl->Label()) {
 //             isTrue = kTRUE;
 //     }
@@ -524,7 +568,10 @@ Double_t AliEMCALTracker::CheckPair
        // check against cut on difference 'alpha - phi'
        Double_t phi = TMath::ATan2(cl->Y(), cl->X());
        phi = AngleDiff(phi, tr->GetAlpha());
-       if (phi < fCutAlphaMin || phi > fCutAlphaMax) return distance;
+       if (phi < fCutAlphaMin || phi > fCutAlphaMax){
+         delete tr;
+         return distance;
+       }
        
        // try to propagate to cluster radius
        // (return the 'distance' value if it fails)
@@ -555,23 +602,29 @@ Double_t AliEMCALTracker::CheckPair
                x0 = 0.0;
        }
        if (fNPropSteps) {
-               Int_t i;
-               Double_t r;
+               Int_t i=0;
+               Double_t r=0.;
                cout.setf(ios::fixed);
                cout.precision(5);
-               if (isTrue) cout << "Init : " << rt << ' ' << x << ' ' << y << ' ' << z << endl;
+               //if (isTrue) cout << "Init : " << rt << ' ' << x << ' ' << y << ' ' << z << endl;
                for (i = 0; i < fNPropSteps; i++) {
                        r = rt + (rc - rt) * ((Double_t)(i+1)/(Double_t)fNPropSteps);
-                       if (!tr->PropagateTo(r, x0, rho)) return distance;
+                       if (!tr->PropagateTo(r, x0, rho)){
+                         delete tr;
+                         return distance;
+                       }
                        tr->GetXYZ(pos);
-                       if (isTrue) cout << "Step : " << r << ' ' << x << ' ' << y << ' ' << z << endl;
+               //      if (isTrue) cout << "Step : " << r << ' ' << x << ' ' << y << ' ' << z << endl;
                }
-               if (isTrue) cout << "Clstr: " << rc << ' ' << cl->X() << ' ' << cl->Y() << ' ' << cl->Z() << endl;
+               //if (isTrue) cout << "Clstr: " << rc << ' ' << cl->X() << ' ' << cl->Y() << ' ' << cl->Z() << endl;
        }
        else {
                // when no steps are used, no correction makes sense
                //if (!tr->PropagateTo(rc, 0.0, 0.0)) return distance;
-               if (!tr->PropagateToGlobal(cl->X(), cl->Y(), cl->Z(), 0.0, 0.0)) return distance;
+               if (!tr->PropagateToGlobal(cl->X(), cl->Y(), cl->Z(), 0.0, 0.0)){
+                 delete tr;
+                 return distance;
+               }
                /*
                Bool_t propOK = kFALSE;
                cout << "START" << endl;
@@ -599,14 +652,14 @@ Double_t AliEMCALTracker::CheckPair
        TVector3 vt(x, y, z);
        Double_t angle = TMath::Abs(vc.Angle(vt)) * TMath::RadToDeg();
        // check: where is the track?
-       Double_t r, phiT, phiC;
-       r = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
-       phiT = TMath::ATan2(pos[1], pos[0]) * TMath::RadToDeg();
-       phiC = vc.Phi() * TMath::RadToDeg();
+//     Double_t r    = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
+//     Double_t phiT = TMath::ATan2(pos[1], pos[0]) * TMath::RadToDeg();
+//     Double_t phiC = vc.Phi() * TMath::RadToDeg();
        //cout << "Propagated R, phiT, phiC = " << r << ' ' << phiT << ' ' << phiC << endl;
        
        if (angle > fCutAngle) {
                //cout << "angle" << endl;
+         delete tr;
                return distance;
        }
                
@@ -614,16 +667,19 @@ Double_t AliEMCALTracker::CheckPair
        x -= cl->X();
        if (TMath::Abs(x) > fCutX) {
                //cout << "cut X" << endl;
+         delete tr;
                return distance;
        }
        y -= cl->Y();
        if (TMath::Abs(y) > fCutY) {
                //cout << "cut Y" << endl;
+         delete tr;
                return distance;
        }
        z -= cl->Z();
        if (TMath::Abs(z) > fCutZ) {
                //cout << "cut Z" << endl;
+         delete tr;
                return distance;
        }
        
@@ -661,7 +717,7 @@ Double_t AliEMCALTracker::CheckPairV2
        
        Double_t distance = 2.0 * fMaxDist;
        
-       Double_t x0, rho;
+       Double_t x0=0., rho=0.;
        if (fTrackCorrMode == kTrackCorrMMB) {
                Double_t pos1[3], pos2[3], param[6];
                tr->GetXYZ(pos1);
@@ -693,7 +749,7 @@ Double_t AliEMCALTracker::CheckPairV2
        TVector3 vc(cl->X(), cl->Y(), cl->Z());
        // rotate the vector in order to put all clusters on a plane intersecting 
        // vertically the X axis; the angle depends on the sector
-       Double_t clusterRot, clusterPhi = vc.Phi() * TMath::RadToDeg();
+       Double_t clusterRot=0., clusterPhi = vc.Phi() * TMath::RadToDeg();
        if (clusterPhi < 0.0) clusterPhi += 360.0;
        if (clusterPhi < 100.0) {
                clusterRot = -90.0;
@@ -719,7 +775,7 @@ Double_t AliEMCALTracker::CheckPairV2
        // compute the 'phi' coordinate of the intersection point to 
        // the EMCAL surface
        Double_t x = vc.X();
-       Double_t y;
+       Double_t y = 0.;
        track->GetYAt(vc.X(), track->GetBz(), y);
        Double_t tmp = x*TMath::Cos(track->GetAlpha()) - y*TMath::Sin(track->GetAlpha());
        y = x*TMath::Sin(track->GetAlpha()) + y*TMath::Cos(track->GetAlpha());
@@ -739,6 +795,7 @@ Double_t AliEMCALTracker::CheckPairV2
        TVector3 vdiff = vt-vc;
                
        // compute differences wr to each coordinate
+       delete track;
        if (vdiff.X() > fCutX) return distance;
        if (vdiff.Y() > fCutY) return distance;
        if (vdiff.Z() > fCutZ) return distance;
@@ -764,12 +821,12 @@ Double_t AliEMCALTracker::CheckPairV3
        
        AliEMCALTrack tr(*track);
        
-       Int_t    sector;
+       Int_t    sector=-1;
        Double_t distance = 2.0 * fMaxDist;
-       Double_t dx, dy, dz;
-       Double_t phi, alpha, slope, tgtXnum, tgtXden, sectorWidth = 20.0 * TMath::DegToRad();
-       Double_t xcurr, xprop, param[6] = {0., 0., 0., 0., 0., 0.}, x0, rho, bz;
-       Double_t x[3], x1[3], x2[3];
+       Double_t dx=0., dy=0., dz=0.;
+       Double_t phi=0., alpha=0., slope=0., tgtXnum=0., tgtXden=0., sectorWidth = 20.0 * TMath::DegToRad();
+       Double_t xcurr=0., xprop=0., param[6] = {0., 0., 0., 0., 0., 0.}, x0=0., rho=0., bz=0.;
+       Double_t x[3]= {0., 0., 0.}, x1[3]= {0., 0., 0.}, x2[3]= {0., 0., 0.};
        
        // get initial track position
        xcurr = tr.GetX();
@@ -791,7 +848,7 @@ Double_t AliEMCALTracker::CheckPairV3
        if (!tr.GetXYZAt(xprop, bz, x2)) return distance;
        //AliKalmanTrack::MeanMaterialBudget(x1, x2, param);
        rho = param[0]*param[4];
-       x0 = param[1];
+       x0  = param[1];
        if (!tr.PropagateTo(xprop, x0, rho)) return distance;
        //if (!tr.PropagateTo(xprop, 0.0, 0.0)) return distance;
        
@@ -815,8 +872,8 @@ Bool_t AliEMCALTracker::PropagateToEMCAL(AliEMCALTrack *tr)
        // Propagates the track to the proximity of the EMCAL surface
        //
        
-       Double_t xcurr, xtemp, xprop = 438.0, step = 10.0, param[6], x0, rho, bz;
-       Double_t x1[3], x2[3];
+       Double_t xcurr=0., xtemp=0., xprop = 438.0, step = 10.0, param[6]= {0., 0., 0., 0., 0., 0.}, x0=0., rho=0., bz=0.;
+       Double_t x1[3]= {0., 0., 0.}, x2[3]= {0., 0., 0.};
        
        // get initial track position
        xcurr = tr->GetX();
@@ -860,11 +917,11 @@ Int_t AliEMCALTracker::CreateMatches()
        
        // initialize counters and indexes
        Int_t count = 0;
-       Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
-       Int_t it, nTracks = fTracks->GetEntries();
+       Int_t ic=0, nClusters = (Int_t)fClusters->GetEntries();
+       Int_t it=0, nTracks = fTracks->GetEntries();
        
        // external loop on clusters, internal loop on tracks
-       Double_t dist;
+       Double_t dist=0.;
        for (ic = 0; ic < nClusters; ic++) {
                AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
                for (it = 0; it < nTracks; it++) {
@@ -882,35 +939,6 @@ Int_t AliEMCALTracker::CreateMatches()
                        }
                }
        }
-               
-       /*
-       // loop on clusters and tracks
-       Int_t icBest;
-       Double_t dist, distBest;
-       for (it = 0; it < nTracks; it++) {
-               AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
-               if (!track) continue;
-               icBest = -1;
-               distBest = fMaxDist;
-               for (ic = 0; ic < nClusters; ic++) {
-                       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
-                       if (!cluster) continue;
-                       dist = CheckPair(track, cluster);
-                       if (dist < distBest) {
-                               distBest = dist;
-                               icBest = ic;
-                       }
-               }
-               if (icBest >= 0) {
-                       track->SetMatchedClusterIndex(icBest);
-                       track->SetMatchedClusterDist(distBest);
-                       count++;
-               }
-               else {
-                       track->SetMatchedClusterIndex(-1);
-               }
-       }
-       */
        
        return count;
 }
@@ -936,8 +964,8 @@ Int_t AliEMCALTracker::SolveCompetitions()
        Int_t count = 0;
        
        // initialize flags to check repetitions
-       Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
-       Int_t it, nTracks = fTracks->GetEntries();
+       Int_t ic=0, nClusters = (Int_t)fClusters->GetEntries();
+       Int_t it=0, nTracks = fTracks->GetEntries();
        Bool_t *usedC = new Bool_t[nClusters];
        Bool_t *usedT = new Bool_t[nTracks];
        for (ic = 0; ic < nClusters; ic++) usedC[ic] = kFALSE;
@@ -959,29 +987,6 @@ Int_t AliEMCALTracker::SolveCompetitions()
                }
        }
        
-       /*
-       Int_t it1, it2, nTracks = (Int_t)fTracks->GetEntries();
-       AliEMCALTrack *track1 = 0, *track2 = 0;
-       for (it1 = 0; it1 < nTracks; it1++) {
-               track1 = (AliEMCALTrack*)fTracks->At(it1);
-               if (!track1) continue;
-               if (track1->GetMatchedClusterIndex() < 0) continue;
-               for (it2 = it1+1; it2 < nTracks; it2++) {
-                       track2 = (AliEMCALTrack*)fTracks->At(it2);
-                       if (!track2) continue;
-                       if (track2->GetMatchedClusterIndex() < 0) continue;
-                       if (track1->GetMatchedClusterIndex() != track2->GetMatchedClusterIndex()) continue;
-                       count++;
-                       if (track1->GetMatchedClusterDist() < track2->GetMatchedClusterDist()) {
-                               track2->SetMatchedClusterIndex(-1);
-                       }
-                       else if (track2->GetMatchedClusterDist() < track1->GetMatchedClusterDist()) {
-                               track1->SetMatchedClusterIndex(-1);
-                       }
-               }
-       }
-       */
-       
        delete [] usedC;
        delete [] usedT;
 
@@ -1000,12 +1005,29 @@ void AliEMCALTracker::UnloadClusters()
        
        Clear();
 }
+//
+//------------------------------------------------------------------------------
+//
+TVector3 AliEMCALTracker::FindExtrapolationPoint(Double_t x,Double_t y,Double_t z, AliESDtrack *track)
+{
+  //Method to determine extrapolation point of track at location x,y,z
+  AliEMCALTrack *tr = new AliEMCALTrack(*track);
+  TVector3 error(-100.,-100.,-100.);
+  if (!tr->PropagateToGlobal(x,y,z, 0.0, 0.0)) {
+    return error;
+  }
+  Double_t pos[3]= {0., 0., 0.};
+  tr->GetXYZ(pos);
+  TVector3 ExTrPos(pos[0],pos[1],pos[2]);
+  return ExTrPos;
+}
+
 //
 //------------------------------------------------------------------------------
 //
 AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCALRecPoint *recPoint)
   : fIndex(index),
-    fLabel(recPoint->GetPrimaryIndex()),
+    fLabel(recPoint->GetPrimaryIndex()), //wrong!  fixed below
     fX(0.),
     fY(0.),
     fZ(0.)
@@ -1020,6 +1042,17 @@ AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliEMCA
        fX = clpos.X();
        fY = clpos.Y();
        fZ = clpos.Z();
+
+       //AliEMCALRecPoint stores the track labels in the parents
+       //list, sorted according to the fractional contribution to the
+       //RecPoint.  The zeroth parent gave the highest contribution
+       Int_t multparent = 0;
+        Int_t *parents = recPoint->GetParents(multparent);
+        if(multparent > 0)
+          fLabel = parents[0];
+        else
+          fLabel = -1;
+
 }
 //
 //------------------------------------------------------------------------------
@@ -1035,7 +1068,7 @@ AliEMCALTracker::AliEMCALMatchCluster::AliEMCALMatchCluster(Int_t index, AliESDC
        // Translates an AliESDCaloCluster object into the internal format.
        // Index of passed cluster in its native array must be specified.
        //
-       Float_t clpos[3];
+       Float_t clpos[3]= {0., 0., 0.};
        caloCluster->GetPosition(clpos);
        
        fX = (Double_t)clpos[0];