(Alberto) Updated track matching algorithm
authorgustavo <gustavo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 6 May 2007 18:11:34 +0000 (18:11 +0000)
committergustavo <gustavo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 6 May 2007 18:11:34 +0000 (18:11 +0000)
EMCAL/AliEMCALTrack.cxx
EMCAL/AliEMCALTrack.h
EMCAL/AliEMCALTracker.cxx
EMCAL/AliEMCALTracker.h

index 7fd8236..4d18886 100755 (executable)
@@ -27,6 +27,8 @@
 
 #include "Riostream.h"
 
+#include "TVector3.h"
+
 #include "AliLog.h"
 #include "AliTracker.h"
 #include "AliESDtrack.h" 
@@ -34,6 +36,8 @@
 #include "AliEMCALTrack.h"
 
 Bool_t AliEMCALTrack::fgUseOuterParams = kTRUE;
+Bool_t AliEMCALTrack::fgCorrectForEL   = kFALSE;
+Bool_t AliEMCALTrack::fgSortByPt       = kTRUE;
 
 ClassImp(AliEMCALTrack)
 //
@@ -42,8 +46,8 @@ ClassImp(AliEMCALTrack)
 AliEMCALTrack::AliEMCALTrack() 
   : AliExternalTrackParam(),
     fClusterIndex(-1),
-    fClusterDist(1000.0),
-    fMass(0.0),
+    fClusterDist(1000.0),          // default: extremely large distance
+    fMass(0.13957018),             // default: pion mass
     fSeedIndex(-1),
     fSeedLabel(-1)
 {
@@ -56,21 +60,6 @@ AliEMCALTrack::AliEMCALTrack()
 //
 //------------------------------------------------------------------------------
 //
-AliEMCALTrack::AliEMCALTrack(const AliEMCALTrack& t) 
-  : AliExternalTrackParam(t),
-    fClusterIndex(t.fClusterIndex),
-    fClusterDist(t.fClusterDist),
-    fMass(t.fMass),
-    fSeedIndex(t.fSeedIndex),
-    fSeedLabel(t.fSeedLabel)
-{
-       //
-       // Copy constructor.
-       //
-       }
-//
-//------------------------------------------------------------------------------
-//
 AliEMCALTrack::AliEMCALTrack(const AliESDtrack& t) 
   : AliExternalTrackParam(),
     fClusterIndex(-1),
@@ -100,7 +89,22 @@ AliEMCALTrack::AliEMCALTrack(const AliESDtrack& t)
        
        // set this accordingly
        Set(x, alpha, params, cov);
-}              
+}
+//
+//------------------------------------------------------------------------------
+//
+AliEMCALTrack::AliEMCALTrack(const AliEMCALTrack& t) 
+  : AliExternalTrackParam(t),
+    fClusterIndex(t.fClusterIndex),
+    fClusterDist(t.fClusterDist),
+    fMass(t.fMass),
+    fSeedIndex(t.fSeedIndex),
+    fSeedLabel(t.fSeedLabel)
+{
+       //
+       // Copy constructor.
+       //
+}
 //
 //------------------------------------------------------------------------------
 //
@@ -127,17 +131,31 @@ AliEMCALTrack& AliEMCALTrack::operator=(const AliEMCALTrack &t)
 Int_t AliEMCALTrack::Compare(const TObject *obj) const 
 {
        //
-       // Tracks compared wrt their distance from matched point
+       // Compare tracks.
+       // How tracks are compared depends on the static flag
+       // "fSortByPt" (boolean):
+       // true  => tracks are compared w.r. to their transverse momentum
+       // false => tracks are compared w.r. to their distance from cluster
        //
        
        AliEMCALTrack *that = (AliEMCALTrack*)obj;
        
-       Double_t thisDist = GetClusterDist();
-       Double_t thatDist = that->GetClusterDist();
+       Double_t thisP[3], thisVal, thatP[3], thatVal;
+       
+       if (fgSortByPt) {
+               this->GetPxPyPz(thisP);
+               that->GetPxPyPz(thatP); 
+               thisVal = TMath::Sqrt(thisP[0]*thisP[0] + thisP[1]*thisP[1]);
+               thatVal = TMath::Sqrt(thatP[0]*thatP[0] + thatP[1]*thatP[1]);
+       }
+       else {
+               thisVal = this->GetClusterDist();
+               thatVal = that->GetClusterDist();
+       }
        
-       if (thisDist > thatDist) return 1;
-       else if (thisDist < thatDist) return -1;
-       return 0;
+       if (thisVal > thatVal) return 1;
+       else if (thisVal < thatVal) return -1;
+       else return 0;
 }
 //
 //------------------------------------------------------------------------------
@@ -145,39 +163,136 @@ Int_t AliEMCALTrack::Compare(const TObject *obj) const
 Double_t AliEMCALTrack::GetBz() const 
 {
        //
-       // returns Bz component of the magnetic field in kG,
-       // at the point where track is actually propagated
+       // Returns Z-component of the magnetic field in kG.
+       // In case it B is not constant, its value is returned
+       // at the current position of the track (local X,Y,Z)
        //
        
-       Double_t r[3];
+       // if magnetic field is constant...
+       if (AliTracker::UniformField()) return AliTracker::GetBz();
        
-       if (AliTracker::UniformField()) {
-               return AliTracker::GetBz();
-       }
-       else {
-               GetXYZ(r);
-               return AliTracker::GetBz(r);
-       }
+       // ...else:
+       Double_t r[3];
+       GetXYZ(r);
+       return AliTracker::GetBz(r);
 }
 //
 //------------------------------------------------------------------------------
 //
-Bool_t AliEMCALTrack::PropagateTo(Double_t xk, Double_t eqWidth, Double_t x0)
+Bool_t AliEMCALTrack::PropagateTo(Double_t xk, Double_t d, Double_t x0)
 {
        //
        // Propagates a track to the plane defined by x='xk'.
-       // Second parameter represents the percent width of radiation length
-       // crossed by the track.
-       // Third parameter is the reference radiation length used (default: Si).
+       // Second parameter is the width (in units of rad. length) crossed by the track.
+       // Third parameter is the reference radiation length used.
        // Track propagation includes computing energy loss (modifies curvature)
-       // and multiple scattering perturbation (alters covariance matrix).
+       // and multiple scattering perturbation (alters covariance matrix), if requested.
        // Method returns kFALSE when something goes wrong with computations.
        //
+       // An additional operation (thanks to Yuri Belikov) is done to check
+       // when the track crosses a sector boundary. If this happens, 
+       // the local track reference frame is adjusted accordingly.
+       //
+               
+       Double_t y;
+       Double_t field = GetBz();
+       Double_t width = TMath::Pi() / 9.0; // width of TPC/TRD/EMCAL sector (= 20 deg)
+       Double_t ymax  = TMath::Abs(xk * TMath::Tan(0.5 * width)); // max allowed Y in local coords at distance xk
+       
+       // first check: try to compute the local 'y' at the distance 'xk':
+       // if this attempt fails, the propagation cannot be done
+       if (!GetYAt(xk, field, y)) return kFALSE;
+       
+       // if is -ymax < y < ymax ==> 'direct' propagation is done;
+       if (TMath::Abs(y) <= ymax) return SimplePropagation(xk, d, x0);
+       
+       // otherwise, try change a sector to find one where the propagation is ok
+       Int_t    i, incr, istart, nloops;
+       Double_t alpha = GetAlpha();
+       incr = (y > ymax) ? 1 : -1;
+       if (alpha < 0.0) alpha += TMath::TwoPi();
+       istart = (Int_t)(alpha / width);
+       for (i = istart, nloops = 0; nloops < 18; i += incr, nloops++) {
+               if (i == 18) i = 0;
+               if (i == -1) i = 17;
+               alpha = ((Double_t)i + 0.5) * width;
+               if (Rotate(alpha)) {
+                       if (GetYAt(xk, field, y)) {
+                               if (TMath::Abs(y) <= ymax) {
+                                       AliInfo(Form("Required change from sector %d to sector %d to succeed in propagation", istart, i));
+                                       return SimplePropagation(xk, d, x0);
+                               }
+                       }
+               }
+       }
+       
+       // if the routine exits from the loop and reaches this point,
+       // it means that none of the rotations succeeded
+       AliWarning("Track propagation fails in every sector. Impossible to propagate.");
+       return kFALSE;
+}
+//
+//------------------------------------------------------------------------------
+//
+Double_t AliEMCALTrack::StraightPropagateTo(Double_t xk, Double_t &x, Double_t &y, Double_t &z)
+{
+       //
+       // Does propagation with a straight line approximation.
+       // This operation does not update track parameters, but it returns a point
+       // in space, according to this propagation, which is stored in
+       // the arguments #2, #3 and #4
+       //
+       
+       Double_t oldX = GetX(), oldY = GetY(), oldZ = GetZ();
+       Double_t newPos[3];
+       
+       newPos[0] = xk;
+       newPos[1] = oldY * xk / oldX;
+       newPos[2] = oldZ * xk / oldX;
+       
+       Double_t cs = TMath::Cos(GetAlpha()), sn = TMath::Sin(GetAlpha());
+       x = newPos[0]*cs - newPos[1]*sn;
+       y = newPos[0]*sn + newPos[1]*cs;
+       z = newPos[2];
+       
+       return newPos[1];
+}
+//
+//------------------------------------------------------------------------------
+//
+Bool_t AliEMCALTrack::PropagateToGlobal(Double_t x, Double_t y, Double_t z, Double_t d, Double_t x0)
+{
+       //
+       // Propagate to a point specified with its global XYZ coordinates.
+       // Here, the correct value of the 'X' parameter to be sent to "PropagateTo" is computed.
+       //
+       
+       TVector3 vc(x, y, z);
+       Double_t width = 20.0; // width of TPC/TRD/EMCAL sector
+       Double_t phi   = vc.Phi() * TMath::RadToDeg();
+       if (phi < 0.0) phi += 360.0;
+       
+       Int_t    isector = (Int_t)(phi / width);
+       Double_t rotation = ((Double_t)isector + 0.5) * width;
+       vc.RotateZ(rotation * TMath::DegToRad());
+       
+       return PropagateTo(vc.X(), d, x0);
+}
+//
+//------------------------------------------------------------------------------
+//
+Bool_t AliEMCALTrack::SimplePropagation(Double_t xk, Double_t d, Double_t x0)
+{
+       //
+       // Recall base class method for track propagation.
+       //
        
        Double_t field = GetBz();
        
+       // propagation...
        if (!AliExternalTrackParam::PropagateTo(xk, field)) return kFALSE;
-       if (!AliExternalTrackParam::CorrectForMaterial(eqWidth, x0, GetMass())) return kFALSE;
-
-       return kTRUE;            
-}     
+       
+       // EL correction is computed only if requested...
+       if (!fgCorrectForEL) return kTRUE;
+       return AliExternalTrackParam::CorrectForMaterial(d, x0, GetMass());
+}
index beffc24..947339f 100755 (executable)
@@ -26,8 +26,8 @@ class AliEMCALTrack : public AliExternalTrackParam
 public:
 
        AliEMCALTrack();
-       AliEMCALTrack(const AliEMCALTrack& t);
        AliEMCALTrack(const AliESDtrack& t);
+       AliEMCALTrack(const AliEMCALTrack& t);
        AliEMCALTrack& operator=(const AliEMCALTrack &t);
        
        Int_t    Compare(const TObject *o) const;
@@ -38,30 +38,39 @@ public:
        Double_t GetMass() const {return fMass;}
        Int_t    GetSeedIndex() const {return fSeedIndex;}
        Int_t    GetSeedLabel() const {return fSeedLabel;}
-               
+       
        Bool_t   IsSortable() const {return kTRUE;}
-       Bool_t   PropagateTo(Double_t xr, Double_t d, Double_t x0=21.82);
+       Bool_t   PropagateTo(Double_t xr, Double_t d = 0.0, Double_t x0 = 21.82);
+       Double_t StraightPropagateTo(Double_t xk, Double_t &x, Double_t &y, Double_t &z);
+       Bool_t   PropagateToGlobal(Double_t x, Double_t y, Double_t z, Double_t d = 0.0, Double_t x0 = 21.82);
                        
-       void     SetClusterIndex(Int_t idx) {fClusterIndex=idx;}
-       void     SetClusterDist(Double_t dist) {fClusterDist=dist;}
-       void     SetMass(Double_t mass) {fMass=mass;}
-       void     SetSeedIndex(Int_t index) {fSeedIndex=index;}
-       void     SetSeedLabel(Int_t label) {fSeedLabel=label;}
+       void     SetClusterIndex(Int_t idx) {fClusterIndex = idx;}
+       void     SetClusterDist(Double_t dist) {fClusterDist = dist;}
+       void     SetMass(Double_t mass) {fMass = mass;}
+       void     SetSeedIndex(Int_t index) {fSeedIndex = index;}
+       void     SetSeedLabel(Int_t label) {fSeedLabel = label;}
                
-       static void SetUseOuterParams(Bool_t doit=kTRUE) {fgUseOuterParams=doit;}
+       static void SetUseOuterParams(Bool_t doit = kTRUE) {fgUseOuterParams = doit;}
+       static void SetCorrectForEL(Bool_t doit = kTRUE) {fgCorrectForEL = doit;}
+       static void SetSortByPt(Bool_t doit = kTRUE) {fgSortByPt = doit;}
 
 protected:
        
+       Bool_t  SimplePropagation(Double_t xk, Double_t d = 0.0, Double_t x0 = 0.0);
+       
        static  Bool_t    fgUseOuterParams;    // use outer parameters from AliESDtrack?
+       static  Bool_t    fgCorrectForEL;      // apply corrections for energy loss?
+       static  Bool_t    fgSortByPt;          // true --> sort by pt, false --> sort by track-cluster distance
                Int_t     fClusterIndex;       // index of matched cluster (if any)
                Double_t  fClusterDist;        // distance between track propagation and matched cluster
-                       Double_t  fMass;               // mass hypothesis (in GeV/c2)
+               Double_t  fMass;               // mass hypothesis (in GeV/c2)
                Int_t     fSeedIndex;          // index of imported ESD track in its owner AliESD
-               Int_t     fSeedLabel;          // label of imported ESD track
+               Int_t     fSeedLabel;          // GEANT label of imported ESD track
+
 private:
        
        ClassDef(AliEMCALTrack, 0) // track implementation for EMCAL matching
 
-};                     
+};
 
 #endif
index e12dd5d..c939483 100644 (file)
@@ -39,6 +39,7 @@
 #include <TString.h>
 #include <TVector3.h>
 #include <TClonesArray.h>
+#include <TGeoMatrix.h>
 
 #include "AliLog.h"
 #include "AliESD.h"
@@ -48,6 +49,8 @@
 #include "AliRunLoader.h"
 #include "AliEMCALTrack.h"
 #include "AliEMCALLoader.h"
+#include "AliEMCALGeometry.h"
+
 #include "AliEMCALTracker.h"
 
 ClassImp(AliEMCALTracker)
@@ -69,7 +72,8 @@ AliEMCALTracker::AliEMCALTracker()
     fX0(1.0),
     fTracks(0),
     fClusters(0),
-    fMatches(0)
+    fMatches(0),
+    fGeom(0)
 {
        //
        // Default constructor.
@@ -96,7 +100,8 @@ AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
     fX0(copy.fX0),
     fTracks((TObjArray*)copy.fTracks->Clone()),
     fClusters((TObjArray*)copy.fClusters->Clone()),
-    fMatches((TList*)copy.fMatches->Clone())
+    fMatches((TList*)copy.fMatches->Clone()),
+    fGeom(copy.fGeom)
 {
        //
        // Copy constructor
@@ -125,11 +130,50 @@ AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
        fClusters = (TObjArray*)copy.fClusters->Clone();
        fMatches = (TList*)copy.fMatches->Clone();
        
+       fGeom = copy.fGeom;
+       
        return (*this);
 }
 //
 //------------------------------------------------------------------------------
 //
+TTree* AliEMCALTracker::SearchTrueMatches()
+{
+       if (!fClusters) return 0;
+       if (fClusters->IsEmpty()) return 0;
+       if (!fTracks) return 0;
+       if (fTracks->IsEmpty()) return 0;
+       
+       TTree *outTree = new TTree("tree", "True matches from event");
+       Int_t indexT, indexC, label;
+       outTree->Branch("indexC", &indexC, "indexC/I");
+       outTree->Branch("indexT", &indexT, "indexT/I");
+       outTree->Branch("label",  &label , "label/I");
+       
+       Double_t dist;
+       Int_t ic, nClusters = (Int_t)fClusters->GetEntries();
+       Int_t it, nTracks = fTracks->GetEntries();
+       
+       for (ic = 0; ic < nClusters; ic++) {
+               AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
+               label = cluster->Label();
+               indexC = cluster->Index();
+               for (it = 0; it < nTracks; it++) {
+                       AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
+                       if (TMath::Abs(track->GetSeedLabel()) != label) continue;
+                       dist = CheckPair(track, cluster);
+                       if (dist <= fMaxDist) {
+                               indexT = track->GetSeedIndex();
+                               outTree->Fill();
+                       }
+               }
+       }
+       
+       return outTree;
+}
+//
+//------------------------------------------------------------------------------
+//
 void AliEMCALTracker::Clear(Option_t* option)
 {
        //
@@ -138,6 +182,7 @@ void AliEMCALTracker::Clear(Option_t* option)
        //
 
        TString opt(option);
+       Bool_t doDelete = opt.Contains("DELETE");
        Bool_t clearTracks = opt.Contains("TRACKS");
        Bool_t clearClusters = opt.Contains("CLUSTERS");
        Bool_t clearMatches = opt.Contains("MATCHES");
@@ -148,19 +193,34 @@ void AliEMCALTracker::Clear(Option_t* option)
        }
        
        if (fTracks != 0x0 && clearTracks) {
-               if (!fTracks->IsEmpty()) fTracks->Delete();
-               delete fTracks;
-               fTracks = 0;
+               if (!doDelete) {
+                       fTracks->Clear();
+               }
+               else {
+                       if (!fTracks->IsEmpty()) fTracks->Delete();
+                       delete fTracks;
+                       fTracks = 0;
+               }
        }
        if (fClusters != 0x0 && clearClusters) {
-               if (!fClusters->IsEmpty()) fClusters->Delete();
-               delete fClusters;
-               fClusters = 0;
+               if (!doDelete) {
+                       fClusters->Clear();
+               }
+               else {
+                       if (!fClusters->IsEmpty()) fClusters->Delete();
+                       delete fClusters;
+                       fClusters = 0;
+               }
        }
        if (fMatches != 0x0 && clearMatches) {
-               if (!fMatches->IsEmpty()) fMatches->Delete();
-               delete fMatches;
-               fMatches = 0;
+               if (!doDelete) {
+                       fMatches->Clear();
+               }
+               else {
+                       if (!fMatches->IsEmpty()) fMatches->Delete();
+                       delete fMatches;
+                       fMatches = 0;
+               }
        }
 }
 //
@@ -178,7 +238,7 @@ Int_t AliEMCALTracker::LoadClusters(TTree *cTree)
 
        TBranch *branch = cTree->GetBranch("EMCALECARP");
        if (!branch) {
-               AliError("can't get the branch with the EMCAL clusters");
+               AliError("Can't get the branch with the EMCAL clusters");
                return 1;
        }
        
@@ -261,21 +321,26 @@ Int_t AliEMCALTracker::LoadTracks(AliESD *esd)
        for (i = 0; i < nTracks; i++) {
                AliESDtrack *esdTrack = esd->GetTrack(i);
                // set by default the value corresponding to "no match"
-               esdTrack->SetEMCALcluster(-99999);
+               esdTrack->SetEMCALcluster(kUnmatched);
 //             if (esdTrack->GetLabel() < 0) continue;
 //             if (!(esdTrack->GetStatus() & AliESDtrack::kTOFout)) continue;
                isKink = kFALSE;
                for (j = 0; j < 3; j++) {
-                       if (esdTrack->GetKinkIndex(j) > 0) isKink = kTRUE;
+                       if (esdTrack->GetKinkIndex(j) != 0) isKink = kTRUE;
                }
                if (isKink) continue;
                AliEMCALTrack *track = new AliEMCALTrack(*esdTrack);
+               track->SetMass(0.13957018);
                // check alpha and reject the tracks which fall outside EMCAL acceptance
                alpha = track->GetAlpha() * TMath::RadToDeg();
                if (alpha >  -155.0 && alpha < 67.0) {
                        delete track;
                        continue;
                }
+//             if (!PropagateToEMCAL(track)) {
+//                     delete track;
+//                     continue;
+//             }
                track->SetSeedIndex(i);
                track->SetSeedLabel(esdTrack->GetLabel());
                fTracks->AddLast(track);
@@ -309,17 +374,18 @@ Int_t AliEMCALTracker::PropagateBack(AliESD* esd)
        // step 1: 
        // if cluster array is empty, cluster are collected
        // from the passed ESD, and work is done with ESDCaloClusters
-       Int_t okLoadClusters;
+       Int_t okLoadClusters, nClusters;
        if (!fClusters || (fClusters && fClusters->IsEmpty())) {
-               AliInfo("Cluster array is empty. Loading clusters...");
                okLoadClusters = LoadClusters(esd);
                if (okLoadClusters) return 2;
        }
+       nClusters = fClusters->GetEntries();
        
        // step 2:
        // collect ESD tracks
-       Int_t okLoadTracks = LoadTracks(esd);
+       Int_t okLoadTracks = LoadTracks(esd), nTracks;
        if (okLoadTracks) return 3;
+       nTracks = fTracks->GetEntries();
        
        // step 3:
        // each track is propagated to the "R" position of each cluster.
@@ -327,11 +393,11 @@ Int_t AliEMCALTracker::PropagateBack(AliESD* esd)
        // IF no clusters lie within the maximum allowed distance, no matches are assigned.
        Int_t nMatches = CreateMatches();
        if (!nMatches) {
-               AliInfo("No good matches found.");
+               AliInfo(Form("#clusters = %d -- #tracks = %d --> No good matches found.", nClusters, nTracks));
                return 4;
        }
        else {
-               AliInfo(Form("Found %d matches", nMatches));
+               AliInfo(Form("#clusters = %d -- #tracks = %d --> Found %d matches.", nClusters, nTracks, nMatches));
        }
        
        // step 4:
@@ -345,7 +411,7 @@ Int_t AliEMCALTracker::PropagateBack(AliESD* esd)
        
        // step 5:
        // save obtained information setting the 'fEMCALindex' field of AliESDtrack object
-       Int_t nSaved = 0, trackID;
+       Int_t nSaved = 0, trackID, nGood = 0, nFake = 0;
        TListIter iter(fMatches);
        AliEMCALMatch *match = 0;
        while ( (match = (AliEMCALMatch*)iter.Next()) ) {
@@ -355,11 +421,13 @@ Int_t AliEMCALTracker::PropagateBack(AliESD* esd)
                trackID = track->GetSeedIndex();
                AliESDtrack *esdTrack = esd->GetTrack(trackID);
                if (!esdTrack) continue;
-               if (esdTrack->GetLabel() == cluster->Label()) {
+               if (TMath::Abs(esdTrack->GetLabel()) == cluster->Label()) {
                        esdTrack->SetEMCALcluster(cluster->Index());
+                       nGood++;
                }
                else {
                        esdTrack->SetEMCALcluster(-cluster->Index());
+                       nFake++;
                }
                nSaved++;
        }
@@ -372,7 +440,7 @@ Int_t AliEMCALTracker::PropagateBack(AliESD* esd)
                AliESDtrack *esdTrack = esd->GetTrack(trackID);
                if (!esdTrack) continue;
                if (clusterID < 0) {
-                       esdTrack->SetEMCALcluster(-99999);
+                       esdTrack->SetEMCALcluster(kUnmatched);
                }
                else {
                        AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(clusterID);
@@ -388,7 +456,7 @@ Int_t AliEMCALTracker::PropagateBack(AliESD* esd)
                }
        }
        */
-       AliInfo(Form("Saved %d matches", nSaved));
+       AliInfo(Form("Saved %d matches [%d good + %d fake]", nSaved, nGood, nFake));
 
        return 0;
 }
@@ -463,7 +531,6 @@ Double_t AliEMCALTracker::CheckPair
        Bool_t isTrue = kFALSE;
 //     if (tr->GetSeedLabel() == cl->Label()) {
 //             isTrue = kTRUE;
-//             cout << "TRUE MATCH!!!" << endl;
 //     }
        
        // copy track into temporary variable
@@ -520,7 +587,25 @@ Double_t AliEMCALTracker::CheckPair
        }
        else {
                // when no steps are used, no correction makes sense
-               if (!tr->PropagateTo(rc, 0.0, 0.0)) return distance;
+               //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;
+               /*
+               Bool_t propOK = kFALSE;
+               cout << "START" << endl;
+               Double_t dist, rCHK, bestDist = 10000000.0;
+               for (Double_t rTMP = rc; rTMP> rc*0.95; rTMP -= 0.1) {
+                       if (!tr->PropagateTo(rTMP)) continue;
+                       propOK = kTRUE;
+                       tr->GetXYZ(pos);
+                       rCHK = TMath::Sqrt(x*x + y*y);
+                       dist = TMath::Abs(rCHK - rc);
+                       cout << rCHK << " vs. " << rc << endl;
+                       
+                       if (TMath::Abs(rCHK - rc) < 0.01) break;
+               }
+               cout << "STOP" << endl;
+               if (!propOK) return distance;
+               */
        }
        
        // get global propagation of track at end of propagation
@@ -530,6 +615,13 @@ Double_t AliEMCALTracker::CheckPair
        TVector3 vc(cl->X(), cl->Y(), cl->Z());
        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();
+       //cout << "Propagated R, phiT, phiC = " << r << ' ' << phiT << ' ' << phiC << endl;
+       
        if (angle > fCutAngle) {
                //cout << "angle" << endl;
                return distance;
@@ -537,17 +629,17 @@ Double_t AliEMCALTracker::CheckPair
                
        // compute differences wr to each coordinate
        x -= cl->X();
-       if (x > fCutX) {
+       if (TMath::Abs(x) > fCutX) {
                //cout << "cut X" << endl;
                return distance;
        }
        y -= cl->Y();
-       if (y > fCutY) {
+       if (TMath::Abs(y) > fCutY) {
                //cout << "cut Y" << endl;
                return distance;
        }
        z -= cl->Z();
-       if (z > fCutZ) {
+       if (TMath::Abs(z) > fCutZ) {
                //cout << "cut Z" << endl;
                return distance;
        }
@@ -675,6 +767,97 @@ Double_t AliEMCALTracker::CheckPairV2
 //
 //------------------------------------------------------------------------------
 //
+Double_t AliEMCALTracker::CheckPairV3
+(AliEMCALTrack *track, AliEMCALMatchCluster *cl)
+{
+       //
+       // Given a track and a cluster,
+       // propagates the first to the radius of the second.
+       // Then, checks the propagation point against all cuts.
+       // If at least a cut is not passed, a valuer equal to 
+       // twice the maximum allowed distance is passed (so the value returned
+       // will not be taken into account when creating matches)
+       //
+       
+       AliEMCALTrack tr(*track);
+       
+       Int_t    sector;
+       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.}, d, x0, rho, bz;
+       Double_t x[3], x1[3], x2[3];
+       
+       // get initial track position
+       xcurr = tr.GetX();
+       
+       // evaluate the EMCAL sector number
+       phi = cl->Phi();
+       if (phi < 0.0) phi += TMath::TwoPi();
+       sector = (Int_t)(phi / sectorWidth);
+       alpha = ((Double_t)sector + 0.5) * sectorWidth;
+       // evaluate the corresponding X for track propagation
+       slope = TMath::Tan(alpha - 0.5*TMath::Pi());
+       tgtXnum = cl->Y() - slope * cl->X();
+       tgtXden = TMath::Sqrt(1.0 + slope*slope);
+       xprop = TMath::Abs(tgtXnum / tgtXden);
+       
+       // propagate by small steps
+       tr.GetXYZ(x1);
+       bz = tr.GetBz();
+       if (!tr.GetXYZAt(xprop, bz, x2)) return distance;
+       //AliKalmanTrack::MeanMaterialBudget(x1, x2, param);
+       d = param[4];
+       rho = param[0];
+       x0 = param[1];
+       if (!tr.PropagateTo(xprop, d*rho/x0, x0)) return distance;
+       //if (!tr.PropagateTo(xprop, 0.0, 0.0)) return distance;
+       
+       // get propagated position at the end
+       tr.GetXYZ(x);
+       dx = TMath::Abs(x[0] - cl->X());
+       dy = TMath::Abs(x[1] - cl->Y());
+       dz = TMath::Abs(x[2] - cl->Z());
+       if (dx > fCutX || dy > fCutY || dz > fCutZ) return distance;
+       
+       distance = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
+       
+       return distance;
+}
+//
+//------------------------------------------------------------------------------
+//
+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], d, x0, rho, bz;
+       Double_t x1[3], x2[3];
+       
+       // get initial track position
+       xcurr = tr->GetX();
+       
+       // propagate by small steps
+       for (xtemp = xcurr + step; xtemp < xprop; xtemp += step) {
+               // to compute material budget, take current position and 
+               // propagated hypothesis without energy loss
+               tr->GetXYZ(x1);
+               bz = tr->GetBz();
+               if (!tr->GetXYZAt(xtemp, bz, x2)) return kFALSE;
+               AliKalmanTrack::MeanMaterialBudget(x1, x2, param);
+               d = param[4];
+               rho = param[0];
+               x0 = param[1];
+               if (!tr->PropagateTo(xtemp, d*rho/x0, x0)) return kFALSE;
+       }
+       
+       return kTRUE;
+}
+//
+//------------------------------------------------------------------------------
+//
 Int_t AliEMCALTracker::CreateMatches()
 {
        //
@@ -702,23 +885,23 @@ Int_t AliEMCALTracker::CreateMatches()
        // external loop on clusters, internal loop on tracks
        Double_t dist;
        for (ic = 0; ic < nClusters; ic++) {
-               cout << "\rMatching cluster " << ic+1 << " of " << nClusters << flush;
                AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
                for (it = 0; it < nTracks; it++) {
                        AliEMCALTrack *track = (AliEMCALTrack*)fTracks->At(it);
                        dist = CheckPair(track, cluster);
+                       //cout << dist << endl;
                        if (dist <= fMaxDist) {
                                AliEMCALMatch *candidate = new AliEMCALMatch;
                                candidate->SetIndexT(it);
                                candidate->SetIndexC(ic);
                                candidate->SetDistance(dist);
+                               candidate->SetPt(track->GetPt());
                                fMatches->Add(candidate);
                                count++;
                        }
                }
        }
-       cout << endl;           
-       
+               
        /*
        // loop on clusters and tracks
        Int_t icBest;
@@ -884,8 +1067,8 @@ Int_t AliEMCALTracker::AliEMCALMatch::Compare(const TObject *obj) const
        
        AliEMCALTracker::AliEMCALMatch *that = (AliEMCALTracker::AliEMCALMatch*)obj;
        
-       Double_t thisDist = fDistance;
-       Double_t thatDist = that->GetDistance();
+       Double_t thisDist = fPt;//fDistance;
+       Double_t thatDist = that->fPt;//that->GetDistance();
        
        if (thisDist > thatDist) return 1;
        else if (thisDist < thatDist) return -1;
@@ -897,7 +1080,8 @@ AliEMCALTracker::AliEMCALMatch::AliEMCALMatch()
     fCanBeSaved(kFALSE), 
     fIndexC(0), 
     fIndexT(0), 
-    fDistance(0.)
+    fDistance(0.),
+    fPt(0.)
 {
   //default constructor
 
@@ -908,7 +1092,8 @@ AliEMCALTracker::AliEMCALMatch::AliEMCALMatch(const AliEMCALMatch& copy)
     fCanBeSaved(copy.fCanBeSaved),
     fIndexC(copy.fIndexC),
     fIndexT(copy.fIndexT),
-    fDistance(copy.fDistance)
+    fDistance(copy.fDistance),
+    fPt(copy.fPt)
 {
   //copy ctor
 }
index 79b7876..ef8ea99 100644 (file)
 #include "AliTracker.h"
 
 class TList;
+class TTree;
 class TObjArray;
 class AliESD;
 class AliESDCaloCluster;
 class AliEMCALTrack;
 class AliEMCALRecPoint;
+class AliEMCALGeometry;
 
 class AliEMCALTracker : public AliTracker 
 {
@@ -50,27 +52,35 @@ public:
        virtual Int_t       RefitInward(AliESD*) {return -1;}
        virtual void        UnloadClusters();
        virtual AliCluster* GetCluster(Int_t) const {return NULL;};
+       TTree*              SearchTrueMatches();
        void                SetCorrection(Double_t rho, Double_t x0) {fRho=rho;fX0=x0;}
        void                SetCutAlpha(Double_t min, Double_t max) {fCutAlphaMin=min;fCutAlphaMax=max;}
        void                SetCutAngle(Double_t value) {fCutAngle=value;}
        void                SetCutX(Double_t value) {fCutX=value;}
        void                SetCutY(Double_t value) {fCutY=value;}
        void                SetCutZ(Double_t value) {fCutZ=value;}
+       void                SetGeometry(AliEMCALGeometry *geom) {fGeom=geom;}
        void                SetMaxDistance(Double_t value) {fMaxDist=value;}
        void                SetNumberOfSteps(Int_t n) {fNPropSteps=n;if(!n)SetTrackCorrectionMode("NONE");}
        void                SetTrackCorrectionMode(Option_t *option);
        
+       enum {
+               kUnmatched = -99999
+       };
+       
        class  AliEMCALMatchCluster : public TObject
        {
        public:
-                          AliEMCALMatchCluster(Int_t ID, AliEMCALRecPoint *recPoint);
-                                  AliEMCALMatchCluster(Int_t ID, AliESDCaloCluster *caloCluster);
-               virtual   ~AliEMCALMatchCluster() { }
+               AliEMCALMatchCluster(Int_t ID, AliEMCALRecPoint *recPoint);
+               AliEMCALMatchCluster(Int_t ID, AliESDCaloCluster *caloCluster);
+               virtual ~AliEMCALMatchCluster() { }
+               //----------------------------------------------------------------------------
                Int_t&     Index() {return fIndex;}
                Int_t&     Label() {return fLabel;}
                Double_t&  X() {return fX;}
                Double_t&  Y() {return fY;}
                Double_t&  Z() {return fZ;}
+               Double_t   Phi() {return TMath::ATan2(fY, fX);}
        private:
                Int_t     fIndex;  // index of cluster in its native container (ESD or TClonesArray)
                Int_t     fLabel;  // track label of assigned cluster
@@ -82,23 +92,27 @@ public:
        class  AliEMCALMatch : public TObject
        {
        public:
-                         AliEMCALMatch();
-                         AliEMCALMatch(const AliEMCALMatch& copy);
-                        ~AliEMCALMatch() { }
+               AliEMCALMatch();
+               AliEMCALMatch(const AliEMCALMatch& copy);
+               virtual ~AliEMCALMatch() { }
+               //----------------------------------------------------------------------------
                Bool_t&   CanBeSaved() {return fCanBeSaved;}
                Int_t     Compare(const TObject *obj) const;
                Double_t  GetDistance() {return fDistance;}
+               Double_t  GetPt() {return fPt;}
                Int_t     GetIndexC() {return fIndexC;}
                Int_t     GetIndexT() {return fIndexT;}
                Bool_t    IsSortable() const {return kTRUE;}
                void      SetIndexC(Int_t icl) {fIndexC=icl;}
                void      SetIndexT(Int_t itr) {fIndexT=itr;}
                void      SetDistance(Double_t dist) {fDistance=dist;}
+               void      SetPt(Double_t pt) {fPt=pt;}
        private:
                Bool_t     fCanBeSaved;  // when true, this match can be saved, otherwise it will not be
                Int_t      fIndexC;      // cluster index in 'fClusters' array
                Int_t      fIndexT;      // track index in 'fTracks' array
                Double_t   fDistance;    // track - cluster distance
+               Double_t   fPt;          // track pt
        };
 
 private:
@@ -106,7 +120,9 @@ private:
        Double_t   AngleDiff(Double_t angle1, Double_t angle2);
        Double_t   CheckPair(AliEMCALTrack *tr, AliEMCALMatchCluster *cluster);
        Double_t   CheckPairV2(AliEMCALTrack *tr, AliEMCALMatchCluster *cluster);
+       Double_t   CheckPairV3(AliEMCALTrack *tr, AliEMCALMatchCluster *cluster);
        Int_t      CreateMatches();
+       Bool_t     PropagateToEMCAL(AliEMCALTrack *tr);
        Int_t      SolveCompetitions();
        
        enum ETrackCorr { 
@@ -129,10 +145,12 @@ private:
        Double_t    fRho;             // energy correction: density
        Double_t    fX0;              // energy correction: radiation length
 
-       TObjArray  *fTracks;          //! collection of ESD tracks
+       TObjArray  *fTracks;          //! collection of tracks
        TObjArray  *fClusters;        //! collection of EMCAL clusters (ESDCaloCluster or EMCALRecPoint)
        TList      *fMatches;         //! collection of matches between tracks and clusters
        
+       AliEMCALGeometry *fGeom;      //! EMCAL geometry
+       
        ClassDef(AliEMCALTracker, 1)  // EMCAL "tracker"
 };