Rongrong: 1. All the extrapolation methods are implemented in AliEMCALTracker; 2...
authormploskon <mploskon@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 19 Oct 2011 15:39:53 +0000 (15:39 +0000)
committermploskon <mploskon@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 19 Oct 2011 15:39:53 +0000 (15:39 +0000)
EMCAL/AliEMCALRecParam.cxx
EMCAL/AliEMCALRecoUtils.cxx
EMCAL/AliEMCALRecoUtils.h
EMCAL/AliEMCALReconstructor.cxx
EMCAL/AliEMCALTracker.cxx
EMCAL/AliEMCALTracker.h

index 5c165c0..9787724 100644 (file)
@@ -411,6 +411,7 @@ AliEMCALRecParam* AliEMCALRecParam::GetLowFluxParam()
   params->SetName("Low Flux - p+p");
   params->SetTitle("Low Flux - p+p");
   params->SetEventSpecie(AliRecoParam::kLowMult);
+  params->SetExtrapolateStep(1);
   
   
   //PID parameters for pp  implemented 
index 91bdf69..a2504a0 100644 (file)
@@ -59,6 +59,7 @@
 #include "AliEMCALTrack.h"
 #include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
 #include "AliEMCALPIDUtils.h"
+#include "AliEMCALTracker.h"
 
 ClassImp(AliEMCALRecoUtils)
   
@@ -75,8 +76,9 @@ AliEMCALRecoUtils::AliEMCALRecoUtils():
   fRejectExoticCluster(kFALSE),           fPIDUtils(),                            fAODFilterMask(32),
   fMatchedTrackIndex(0x0),                fMatchedClusterIndex(0x0), 
   fResidualEta(0x0), fResidualPhi(0x0),   fCutEtaPhiSum(kTRUE),                   fCutEtaPhiSeparate(kFALSE), 
-  fCutR(0.1),                             fCutEta(0.025),                         fCutPhi(0.05), 
-  fMass(0.139),                           fStep(10),
+  fCutR(0.1),                             fCutEta(0.025),                         fCutPhi(0.05),
+  fClusterWindow(100),
+  fMass(0.139),                           fStepSurface(10.),                      fStepCluster(5.),
   fTrackCutsType(kLooseCut),              fCutMinTrackPt(0),                      fCutMinNClusterTPC(-1), 
   fCutMinNClusterITS(-1),                 fCutMaxChi2PerClusterTPC(1e10),         fCutMaxChi2PerClusterITS(1e10),
   fCutRequireTPCRefit(kFALSE),            fCutRequireITSRefit(kFALSE),            fCutAcceptKinkDaughters(kFALSE),
@@ -148,7 +150,8 @@ AliEMCALRecoUtils::AliEMCALRecoUtils(const AliEMCALRecoUtils & reco)
   fResidualPhi(        reco.fResidualPhi?        new TArrayF(*reco.fResidualPhi):0x0),
   fCutEtaPhiSum(reco.fCutEtaPhiSum),                         fCutEtaPhiSeparate(reco.fCutEtaPhiSeparate), 
   fCutR(reco.fCutR),        fCutEta(reco.fCutEta),           fCutPhi(reco.fCutPhi),
-  fMass(reco.fMass),        fStep(reco.fStep),
+  fClusterWindow(reco.fClusterWindow),
+  fMass(reco.fMass),        fStepSurface(reco.fStepSurface), fStepCluster(reco.fStepCluster),
   fTrackCutsType(reco.fTrackCutsType),                       fCutMinTrackPt(reco.fCutMinTrackPt), 
   fCutMinNClusterTPC(reco.fCutMinNClusterTPC),               fCutMinNClusterITS(reco.fCutMinNClusterITS), 
   fCutMaxChi2PerClusterTPC(reco.fCutMaxChi2PerClusterTPC),   fCutMaxChi2PerClusterITS(reco.fCutMaxChi2PerClusterITS),
@@ -214,8 +217,10 @@ AliEMCALRecoUtils & AliEMCALRecoUtils::operator = (const AliEMCALRecoUtils & rec
   fCutR                      = reco.fCutR;
   fCutEta                    = reco.fCutEta;
   fCutPhi                    = reco.fCutPhi;
+  fClusterWindow             = reco.fClusterWindow;
   fMass                      = reco.fMass;
-  fStep                      = reco.fStep;
+  fStepSurface               = reco.fStepSurface;
+  fStepCluster               = reco.fStepCluster;
   fRejectExoticCluster       = reco.fRejectExoticCluster;
 
   fTrackCutsType             = reco.fTrackCutsType;
@@ -940,7 +945,7 @@ void AliEMCALRecoUtils::RecalibrateCellTime(const Int_t absId, const Int_t bc, D
   
 }
   
-//________________________________________________________________________________________________________________
+//__________________________________________________
 void AliEMCALRecoUtils::RecalculateClusterPosition(AliEMCALGeometry *geom, AliVCaloCells* cells, AliVCluster* clu)
 {
   //For a given CaloCluster recalculates the position for a given set of misalignment shifts and puts it again in the CaloCluster.
@@ -1371,6 +1376,18 @@ void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  Al
   
   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
+
+  TObjArray *clusterArray = 0x0;
+  if(!clusterArr)
+    {
+      clusterArray = new TObjArray(event->GetNumberOfCaloClusters());
+      for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
+       {
+         AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
+         if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
+         clusterArray->AddAt(cluster,icl);
+       }
+    }
   
   Int_t    matched=0;
   Double_t cv[21];
@@ -1380,22 +1397,29 @@ void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  Al
     AliExternalTrackParam *trackParam = 0;
     
     //If the input event is ESD, the starting point for extrapolation is TPCOut, if available, or TPCInner 
+    AliESDtrack *esdTrack = 0;
+    AliAODTrack *aodTrack = 0;
     if(esdevent)
     {
-      AliESDtrack *esdTrack = esdevent->GetTrack(itr);
-      if(!esdTrack || !IsAccepted(esdTrack)) continue;
+      esdTrack = esdevent->GetTrack(itr);
+      if(!esdTrack) continue;
+      if(!IsAccepted(esdTrack)) continue;
       if(esdTrack->Pt()<fCutMinTrackPt) continue;
+      Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
+      if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
       trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
     }
     
     //If the input event is AOD, the starting point for extrapolation is at vertex
-    //AOD tracks are selected according to its bit.
+    //AOD tracks are selected according to its filterbit.
     else if(aodevent)
     {
-      AliAODTrack *aodTrack = aodevent->GetTrack(itr);
+      aodTrack = aodevent->GetTrack(itr);
       if(!aodTrack) continue;
       if(!aodTrack->TestFilterMask(fAODFilterMask)) continue; //Select AOD tracks that fulfill GetStandardITSTPCTrackCuts2010()
       if(aodTrack->Pt()<fCutMinTrackPt) continue;
+      Double_t phi = aodTrack->Phi()*TMath::RadToDeg();
+      if(TMath::Abs(aodTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
       Double_t pos[3],mom[3];
       aodTrack->GetXYZ(pos);
       aodTrack->GetPxPyPz(mom);
@@ -1407,162 +1431,160 @@ void AliEMCALRecoUtils::FindMatches(AliVEvent *event,TObjArray * clusterArr,  Al
     else
     {
       printf("Wrong input data type! Should be \"AOD\" or \"ESD\"\n");
+      if(clusterArray)
+       {
+         clusterArray->Clear();
+         delete clusterArray;
+       }
       return;
     }
     
     if(!trackParam) continue;
-    
-    Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
+
+    //Extrapolate the track to EMCal surface
+    AliExternalTrackParam emcalParam(*trackParam);
+    Double_t eta, phi;
+    if(!AliEMCALTracker::ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) 
+      {
+       if(aodevent && trackParam) delete trackParam;
+       continue;
+      }
+
+    if(esdevent)
+      {
+       esdTrack->SetOuterParam(&emcalParam,AliExternalTrackParam::kMultSec);
+      }
+
+    if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad())
+      {
+       if(aodevent && trackParam) delete trackParam;
+       continue;
+      }
+
+
+    //Find matched clusters
     Int_t index = -1;
-    if(!clusterArr){// get clusters from event
-      for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
+    Float_t dEta = -999, dPhi = -999;
+    if(!clusterArr)
       {
-        AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
-        if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
-        AliExternalTrackParam trkPamTmp(*trackParam);//Retrieve the starting point every time before the extrapolation 
-        Float_t tmpEta=-999, tmpPhi=-999;
-        if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
-        if(fCutEtaPhiSum)
-        {
-          Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
-          if(tmpR<dRMax)
-          {
-            dRMax=tmpR;
-            dEtaMax=tmpEta;
-            dPhiMax=tmpPhi;
-            index=icl;
-          }
-        }
-        else if(fCutEtaPhiSeparate)
-        {
-          if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
-          {
-            dEtaMax = tmpEta;
-            dPhiMax = tmpPhi;
-            index=icl;
-          }
-        }
-        else
-        {
-          printf("Error: please specify your cut criteria\n");
-          printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
-          printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
-          if(aodevent && trackParam) delete trackParam;
-          return;
-        }
-      }//cluster loop
-    } 
-    else { // external cluster array, not from ESD event
-      for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
+       index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArray, dEta, dPhi);  
+      }
+    else
       {
-        AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
-        if(!cluster){ 
-          AliInfo("Cluster not found!!!");
-          continue;
-        }
-        if(!cluster->IsEMCAL()) continue;
-        AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
-        Float_t tmpEta=-999, tmpPhi=-999;
-        if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
-        if(fCutEtaPhiSum)
-        {
-          Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
-          if(tmpR<dRMax)
-          {
-            dRMax=tmpR;
-            dEtaMax=tmpEta;
-            dPhiMax=tmpPhi;
-            index=icl;
-          }
-        }
-        else if(fCutEtaPhiSeparate)
-        {
-          if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
-          {
-            dEtaMax = tmpEta;
-            dPhiMax = tmpPhi;
-            index=icl;
-          }
-        }
-        else
-        {
-          printf("Error: please specify your cut criteria\n");
-          printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
-          printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
-          if(aodevent && trackParam) delete trackParam;
-          return;
-        }
-      }//cluster loop
-    }// external list of clusters
+       index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);  
+      }  
     
     if(index>-1)
     {
-      fMatchedTrackIndex  ->AddAt(itr,matched);
-      fMatchedClusterIndex->AddAt(index,matched);
-      fResidualEta          ->AddAt(dEtaMax,matched);
-      fResidualPhi          ->AddAt(dPhiMax,matched);
+      fMatchedTrackIndex   ->AddAt(itr,matched);
+      fMatchedClusterIndex ->AddAt(index,matched);
+      fResidualEta         ->AddAt(dEta,matched);
+      fResidualPhi         ->AddAt(dPhi,matched);
       matched++;
     }
     if(aodevent && trackParam) delete trackParam;
   }//track loop
+
+  if(clusterArray)
+    {
+      clusterArray->Clear();
+      delete clusterArray;
+    }
   
   AliDebug(2,Form("Number of matched pairs = %d !\n",matched));
   
-  fMatchedTrackIndex  ->Set(matched);
-  fMatchedClusterIndex->Set(matched);
-  fResidualPhi          ->Set(matched);
-  fResidualEta          ->Set(matched);
+  fMatchedTrackIndex   ->Set(matched);
+  fMatchedClusterIndex ->Set(matched);
+  fResidualPhi         ->Set(matched);
+  fResidualEta         ->Set(matched);
 }
 
 //________________________________________________________________________________
-Int_t AliEMCALRecoUtils::FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom)
+Int_t AliEMCALRecoUtils::FindMatchedClusterInEvent(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom, Float_t &dEta, Float_t &dPhi)
 {
   //
   // This function returns the index of matched cluster to input track
   // Returns -1 if no match is found
-  
-  
-  Float_t dRMax = fCutR, dEtaMax = fCutEta, dPhiMax = fCutPhi;
   Int_t index = -1;
-  
+  Double_t phiV = track->Phi()*TMath::RadToDeg();
+  if(TMath::Abs(track->Eta())>0.8 || phiV <= 20 || phiV >= 240 ) return index;
   AliExternalTrackParam *trackParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
-  
-  if(!trackParam) return index;          
+  if(!trackParam) return index;
+  AliExternalTrackParam emcalParam(*trackParam);
+  Double_t eta, phi;
+  if(!AliEMCALTracker::ExtrapolateTrackToEMCalSurface(&emcalParam, 430., fMass, fStepSurface, eta, phi)) return index;
+  if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
+
+  TObjArray *clusterArr = new TObjArray(event->GetNumberOfCaloClusters());
+
   for(Int_t icl=0; icl<event->GetNumberOfCaloClusters(); icl++)
   {
     AliVCluster *cluster = (AliVCluster*) event->GetCaloCluster(icl);
-    if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;  
-    AliExternalTrackParam trkPamTmp (*trackParam);//Retrieve the starting point every time before the extrapolation
-    Float_t tmpEta=-999, tmpPhi=-999;
-    if(!ExtrapolateTrackToCluster(&trkPamTmp, cluster, tmpEta, tmpPhi)) continue;
-    if(fCutEtaPhiSum)
+    if(geom && !IsGoodCluster(cluster,geom,(AliVCaloCells*)event->GetEMCALCells())) continue;
+    clusterArr->AddAt(cluster,icl);
+  }
+
+  index = FindMatchedClusterInClusterArr(&emcalParam, &emcalParam, clusterArr, dEta, dPhi);  
+  clusterArr->Clear();
+  delete clusterArr;
+  
+  return index;
+}
+
+//________________________________________________________________________________
+Int_t  AliEMCALRecoUtils::FindMatchedClusterInClusterArr(AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, TObjArray * clusterArr, Float_t &dEta, Float_t &dPhi)
+{
+  dEta=-999, dPhi=-999;
+  Float_t dRMax = fCutR, dEtaMax=fCutEta, dPhiMax=fCutPhi;
+  Int_t index = -1;
+  Double_t tmpEta=-999, tmpPhi=-999;
+
+  Double_t exPos[3] = {0.,0.,0.};
+  if(!emcalParam->GetXYZ(exPos)) return index;
+
+  Float_t clsPos[3] = {0.,0.,0.};
+  for(Int_t icl=0; icl<clusterArr->GetEntriesFast(); icl++)
     {
-      Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
-      if(tmpR<dRMax)
+      AliVCluster *cluster = dynamic_cast<AliVCluster*> (clusterArr->At(icl)) ;
+      if(!cluster || !cluster->IsEMCAL()) continue;
+      cluster->GetPosition(clsPos);
+      Double_t dR = TMath::Sqrt(TMath::Power(exPos[0]-clsPos[0],2)+TMath::Power(exPos[1]-clsPos[1],2)+TMath::Power(exPos[2]-clsPos[2],2));
+      if(dR > fClusterWindow) continue;
+
+      AliExternalTrackParam trkPamTmp (*trkParam);//Retrieve the starting point every time before the extrapolation
+      if(!AliEMCALTracker::ExtrapolateTrackToCluster(&trkPamTmp, cluster, fMass, fStepCluster, tmpEta, tmpPhi)) continue;
+      if(fCutEtaPhiSum)
+        {
+          Float_t tmpR=TMath::Sqrt(tmpEta*tmpEta + tmpPhi*tmpPhi);
+          if(tmpR<dRMax)
            {
              dRMax=tmpR;
              dEtaMax=tmpEta;
              dPhiMax=tmpPhi;
              index=icl;
            }
-    }
-    else if(fCutEtaPhiSeparate)
-    {
-      if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
+        }
+      else if(fCutEtaPhiSeparate)
+        {
+          if(TMath::Abs(tmpEta)<TMath::Abs(dEtaMax) && TMath::Abs(tmpPhi)<TMath::Abs(dPhiMax))
            {
              dEtaMax = tmpEta;
              dPhiMax = tmpPhi;
              index=icl;
            }
+        }
+      else
+        {
+          printf("Error: please specify your cut criteria\n");
+          printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
+          printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
+          return index;
+        }
     }
-    else
-    {
-      printf("Error: please specify your cut criteria\n");
-      printf("To cut on sqrt(dEta^2+dPhi^2), use: SwitchOnCutEtaPhiSum()\n");
-      printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
-      return -1;
-    }
-  }//cluster loop
+
+  dEta=dEtaMax;
+  dPhi=dPhiMax;
+
   return index;
 }
 
@@ -1572,24 +1594,11 @@ Bool_t  AliEMCALRecoUtils::ExtrapolateTrackToCluster(AliExternalTrackParam *trkP
   //
   //Return the residual by extrapolating a track to a cluster
   //
-  if(!trkParam || !cluster) return kFALSE;
-  Float_t clsPos[3];
-  Double_t trkPos[3];
-  cluster->GetPosition(clsPos); //Has been recalculated
-  TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
-  Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
-  vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
-  if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), fMass, fStep,kTRUE, 0.8, -1)) return kFALSE;
-  trkParam->GetXYZ(trkPos); //Get the extrapolated global position
-
-  TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
-  TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
-
-  // track cluster matching
-  tmpPhi = clsPosVec.DeltaPhi(trkPosVec); // tmpPhi is between -pi and pi
-  tmpEta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
-
-  return kTRUE;
+
+  Double_t dEta = -999, dPhi = -999;
+  Bool_t result = AliEMCALTracker::ExtrapolateTrackToCluster(trkParam, cluster, fMass, fStepCluster, dEta, dPhi);
+  tmpEta=dEta, tmpPhi=dPhi;
+  return result;
 }
 
 //________________________________________________________________________________
@@ -1906,11 +1915,11 @@ void AliEMCALRecoUtils::Print(const Option_t *) const
   printf("Matching criteria: ");
   if(fCutEtaPhiSum)
     {
-      printf("sqrt(dEta^2+dPhi^2)<%2.2f\n",fCutR);
+      printf("sqrt(dEta^2+dPhi^2)<%4.3f\n",fCutR);
     }
   else if(fCutEtaPhiSeparate)
     {
-      printf("dEta<%2.2f, dPhi<%2.2f\n",fCutEta,fCutPhi);
+      printf("dEta<%4.3f, dPhi<%4.3f\n",fCutEta,fCutPhi);
     }
   else
     {
@@ -1920,7 +1929,8 @@ void AliEMCALRecoUtils::Print(const Option_t *) const
       printf("To cut on dEta and dPhi separately, use: SwitchOnCutEtaPhiSeparate()\n");
     }
 
-  printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step = %2.2f[cm]\n",fMass,fStep);
+  printf("Mass hypothesis = %2.3f [GeV/c^2], extrapolation step to surface = %2.2f[cm], step to cluster = %2.2f[cm]\n",fMass,fStepSurface, fStepCluster);
+  printf("Cluster selection window: dR < %2.0f\n",fClusterWindow);
 
   printf("Track cuts: \n");
   printf("Minimum track pT: %1.2f\n",fCutMinTrackPt);
index 21b2cb2..9b83e7c 100644 (file)
@@ -245,7 +245,9 @@ public:
   Bool_t   ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Float_t &tmpEta, Float_t &tmpPhi);
 
   void     FindMatches(AliVEvent *event, TObjArray * clusterArr=0x0, AliEMCALGeometry *geom=0x0);
-  Int_t    FindMatchedCluster(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom);
+  Int_t    FindMatchedClusterInEvent(AliESDtrack *track, AliVEvent *event, AliEMCALGeometry *geom, Float_t &dEta, Float_t &dPhi);
+  Int_t    FindMatchedClusterInClusterArr(AliExternalTrackParam *emcalParam, AliExternalTrackParam *trkParam, TObjArray * clusterArr, Float_t &dEta, Float_t &dPhi);
+
   UInt_t   FindMatchedPosForCluster(Int_t clsIndex) const;
   UInt_t   FindMatchedPosForTrack(Int_t trkIndex)   const;
   
@@ -266,15 +268,19 @@ public:
   Float_t  GetCutR()                            const { return fCutR                  ; }
   Float_t  GetCutEta()                          const { return fCutEta                ; }
   Float_t  GetCutPhi()                          const { return fCutPhi                ; }
+  Double_t GetClusterWindow()                   const { return fClusterWindow         ; }
   void     SetCutR(Float_t cutR)                      { fCutR   = cutR                ; }
   void     SetCutEta(Float_t cutEta)                  { fCutEta = cutEta              ; }
   void     SetCutPhi(Float_t cutPhi)                  { fCutPhi = cutPhi              ; }
+  void     SetClusterWindow(Double_t window)          { fClusterWindow = window       ; }
   void     SetCutZ(Float_t cutZ)                      { printf("Obsolete fucntion of cutZ=%1.1f\n",cutZ) ; } //Obsolete
 
   Double_t GetMass()                            const { return fMass                  ; }
-  Double_t GetStep()                            const { return fStep                  ; }
+  Double_t GetStep()                            const { return fStepCluster           ; }
+  Double_t GetStepSurface()                     const { return fStepSurface           ; }
   void     SetMass(Double_t mass)                     { fMass = mass                  ; }
-  void     SetStep(Double_t step)                     { fStep = step                  ; }
+  void     SetStep(Double_t step)                     { fStepCluster = step           ; }
+  void     SetStepSurface(Double_t step)              { fStepSurface = step           ; }
  
   //Cluster cut
   Bool_t   IsGoodCluster(AliVCluster *cluster, AliEMCALGeometry *geom, AliVCaloCells* cells);
@@ -319,8 +325,7 @@ public:
   Bool_t   GetDCAToVertex2D()                  const { return fCutDCAToVertex2D        ; }
 
 
-private:
-  
+private:  
   //Position recalculation
   Float_t    fMisalTransShift[15];       // Shift parameters
   Float_t    fMisalRotShift[15];         // Shift parameters
@@ -377,8 +382,10 @@ private:
   Float_t    fCutR;                      // sqrt(dEta^2+dPhi^2) cut on matching
   Float_t    fCutEta;                    // dEta cut on matching
   Float_t    fCutPhi;                    // dPhi cut on matching
+  Double_t   fClusterWindow;             // Select clusters in the window to be matched
   Double_t   fMass;                      // Mass hypothesis of the track
-  Double_t   fStep;                      // Length of each step used in extrapolation in the unit of cm.
+  Double_t   fStepSurface;               // Length of step to extrapolate tracks to EMCal surface
+  Double_t   fStepCluster;               // Length of step to extrapolate tracks to clusters
 
   // Track cuts  
   Int_t      fTrackCutsType;             // Esd track cuts type for matching
@@ -392,9 +399,9 @@ private:
   Bool_t     fCutAcceptKinkDaughters;    // Accepting kink daughters?
   Float_t    fCutMaxDCAToVertexXY;       // Track-to-vertex cut in max absolute distance in xy-plane
   Float_t    fCutMaxDCAToVertexZ;        // Track-to-vertex cut in max absolute distance in z-plane
-  Bool_t     fCutDCAToVertex2D;          // If true a 2D DCA cut is made. Tracks are accepted if sqrt((DCAXY / fCutMaxDCAToVertexXY)^2 + (DCAZ / fCutMaxDCAToVertexZ)^2) < 1 AND sqrt((DCAXY / fCutMinDCAToVertexXY)^2 + (DCAZ / fCutMinDCAToVertexZ)^2) > 1
+  Bool_t     fCutDCAToVertex2D;          // If true a 2D DCA cut is made.
   
-  ClassDef(AliEMCALRecoUtils, 14)
+  ClassDef(AliEMCALRecoUtils, 15)
   
 };
 
index feaaecd..15d62ba 100644 (file)
@@ -49,6 +49,7 @@
 #include "AliEMCALClusterizerNxN.h"
 #include "AliEMCALRecPoint.h"
 #include "AliEMCALPID.h"
+#include "AliEMCALTracker.h"
 #include "AliRawReader.h"
 #include "AliCDBEntry.h"
 #include "AliCDBManager.h"
@@ -699,6 +700,9 @@ Bool_t AliEMCALReconstructor::CalculateResidual(AliESDtrack *track, AliESDCaloCl
 
   // If the esdFriend is available, use the TPCOuter point as the starting point of extrapolation
   // Otherwise use the TPCInner point
+
+  dEta = -999, dPhi = -999;
+
   AliExternalTrackParam *trkParam = 0;
   const AliESDfriendTrack*  friendTrack = track->GetFriendTrack();
   if(friendTrack && friendTrack->GetTPCOut())
@@ -707,28 +711,8 @@ Bool_t AliEMCALReconstructor::CalculateResidual(AliESDtrack *track, AliESDCaloCl
     trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
   if(!trkParam) return kFALSE;
 
-  //Perform extrapolation
-  Double_t trkPos[3];
-  Float_t  clsPos[3];
-
   AliExternalTrackParam trkParamTmp (*trkParam);
-  cluster->GetPosition(clsPos);
-  TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
-  Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
-  //Rotate to proper local coordinate
-  vec.RotateZ(-alpha); 
-  //extrapolation is done here
-  if(!AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, vec.X(), track->GetMass(), GetRecParam()->GetExtrapolateStep(), kTRUE, 0.8, -1)) 
-    return kFALSE; 
-
-  //Calculate the residuals
-  trkParamTmp.GetXYZ(trkPos); 
-   
-  TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
-  TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
-
-  dPhi = clsPosVec.DeltaPhi(trkPosVec);
-  dEta = clsPosVec.Eta()-trkPosVec.Eta();
+  if(!AliEMCALTracker::ExtrapolateTrackToCluster(&trkParamTmp, cluster, track->GetMass(), GetRecParam()->GetExtrapolateStep(), dEta, dPhi)) return kFALSE;
 
   return kTRUE;
 }
index a405668..81daa8e 100644 (file)
@@ -64,17 +64,18 @@ ClassImp(AliEMCALTracker)
 //------------------------------------------------------------------------------
 //
 AliEMCALTracker::AliEMCALTracker() 
-  : AliTracker(),
-    fCutPt(0),
-    fCutNITS(0),
-    fCutNTPC(50),
-    fStep(50),
-    fTrackCorrMode(kTrackCorrMMB),     
-    fCutEta(0.025),
-    fCutPhi(0.05),
-    fTracks(0),
-    fClusters(0),
-    fGeom(0)
+: AliTracker(),
+  fCutPt(0),
+  fCutNITS(0),
+  fCutNTPC(50),
+  fStep(10),
+  fTrackCorrMode(kTrackCorrMMB),
+  fClusterWindow(50),
+  fCutEta(0.025),
+  fCutPhi(0.05),
+  fTracks(0),
+  fClusters(0),
+  fGeom(0)
 {
   //
   // Default constructor.
@@ -94,6 +95,7 @@ AliEMCALTracker::AliEMCALTracker(const AliEMCALTracker& copy)
     fCutNTPC(copy.fCutNTPC),
     fStep(copy.fStep),
     fTrackCorrMode(copy.fTrackCorrMode),
+    fClusterWindow(copy.fClusterWindow),
     fCutEta(copy.fCutEta),
     fCutPhi(copy.fCutPhi),
     fTracks((TObjArray*)copy.fTracks->Clone()),
@@ -116,6 +118,7 @@ AliEMCALTracker& AliEMCALTracker::operator=(const AliEMCALTracker& copy)
   //
 
   fCutPt  = copy.fCutPt;
+  fClusterWindow = copy.fClusterWindow;
   fCutEta = copy.fCutEta;
   fCutPhi = copy.fCutPhi;      
   fStep = copy.fStep;
@@ -149,10 +152,10 @@ void AliEMCALTracker::InitParameters()
  
     fCutEta  =  recParam->GetMthCutEta();
     fCutPhi  =  recParam->GetMthCutPhi();
-    fStep    =   recParam->GetExtrapolateStep();
+    fStep    =  recParam->GetExtrapolateStep();
     fCutPt   =  recParam->GetTrkCutPt();
-    fCutNITS = recParam->GetTrkCutNITS();
-    fCutNTPC = recParam->GetTrkCutNTPC();
+    fCutNITS =  recParam->GetTrkCutNITS();
+    fCutNTPC =  recParam->GetTrkCutNTPC();
   }
        
 }
@@ -272,25 +275,18 @@ Int_t AliEMCALTracker::LoadTracks(AliESDEvent *esd)
   fTracks = new TObjArray(0);
        
   Int_t nTracks = esd->GetNumberOfTracks();
-  Bool_t isKink=kFALSE;
+  //Bool_t isKink=kFALSE;
   for (Int_t i = 0; i < nTracks; i++) 
     {
       AliESDtrack *esdTrack = esd->GetTrack(i);
       // set by default the value corresponding to "no match"
       esdTrack->SetEMCALcluster(kUnmatched);
+      esdTrack->ResetStatus(AliESDtrack::kEMCALmatch);
 
       //Select good quaulity tracks
       if(esdTrack->Pt()<fCutPt) continue;
       if(esdTrack->GetNcls(1)<fCutNTPC)continue;
 
-      //Reject kink daughters
-      isKink = kFALSE;
-      for (Int_t j = 0; j < 3; j++)
-       {
-         if (esdTrack->GetKinkIndex(j) != 0) isKink = kTRUE;
-       }
-      if (isKink) continue;
-
       //Loose geometric cut
       Double_t phi = esdTrack->Phi()*TMath::RadToDeg();
       if(TMath::Abs(esdTrack->Eta())>0.8 || phi <= 20 || phi >= 240 ) continue;
@@ -401,40 +397,106 @@ Int_t AliEMCALTracker::FindMatchedCluster(AliESDtrack *track)
     trkParam = const_cast<AliExternalTrackParam*>(track->GetInnerParam());
   if(!trkParam) return index;
 
+
+  AliExternalTrackParam trkParamTmp(*trkParam);
+  Double_t eta, phi;
+  if(!ExtrapolateTrackToEMCalSurface(&trkParamTmp, 430., track->GetMass(), fStep, eta, phi)) return index;
+  if(TMath::Abs(eta)>0.75 || (phi) < 70*TMath::DegToRad() || (phi) > 190*TMath::DegToRad()) return index;
+
   //Perform extrapolation
   Double_t trkPos[3];
+  trkParamTmp.GetXYZ(trkPos);
   Int_t nclusters = fClusters->GetEntries();
   for(Int_t ic=0; ic<nclusters; ic++)
     {
-      //AliExternalTrackParam *trkParamTmp = new AliExternalTrackParam(*trkParam);
-      AliExternalTrackParam trkParamTmp(*trkParam);
       AliEMCALMatchCluster *cluster = (AliEMCALMatchCluster*)fClusters->At(ic);
-      TVector3 vec(cluster->X(),cluster->Y(),cluster->Z());
-      Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
-      //Rotate to proper local coordinate
-      vec.RotateZ(-alpha); 
-      //extrapolation is done here
-      if(!AliTrackerBase::PropagateTrackToBxByBz(&trkParamTmp, vec.X(), track->GetMass(), fStep, kTRUE, 0.8, -1)) continue; 
-
-      //Calculate the residuals
-      trkParamTmp.GetXYZ(trkPos);        
-      TVector3 clsPosVec(cluster->X(),cluster->Y(),cluster->Z());
-      TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
-
-      Double_t tmpPhi = clsPosVec.DeltaPhi(trkPosVec);
-      Double_t tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
-  
+      Float_t clsPos[3] = {cluster->X(),cluster->Y(),cluster->Z()};
+      Double_t dR = TMath::Sqrt(TMath::Power(trkPos[0]-clsPos[0],2)+TMath::Power(trkPos[1]-clsPos[1],2)+TMath::Power(trkPos[2]-clsPos[2],2));
+      if(dR > fClusterWindow) continue;
+      
+      AliExternalTrackParam trkParTmp(trkParamTmp);
+
+      Double_t tmpEta, tmpPhi;
+      if(!ExtrapolateTrackToPosition(&trkParTmp, clsPos,track->GetMass(), fStep, tmpEta, tmpPhi)) continue;
       if(TMath::Abs(tmpPhi)<TMath::Abs(maxPhi) && TMath::Abs(tmpEta)<TMath::Abs(maxEta))
         {
           maxPhi=tmpPhi;
           maxEta=tmpEta;
           index=ic;
         }
-      //delete trkParamTmp;
       }
   return index;
 }
 
+
+//
+//------------------------------------------------------------------------------
+//
+Bool_t AliEMCALTracker::ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam, Double_t emcalR, Double_t mass, Double_t step, Double_t &eta, Double_t &phi)
+{
+  eta = -999, phi = -999;
+  if(!trkParam) return kFALSE;
+  if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, emcalR, mass, step, kTRUE, 0.8, -1)) return kFALSE;
+  Double_t trkPos[3] = {0.,0.,0.};
+  if(!trkParam->GetXYZ(trkPos)) return kFALSE;
+  TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+  eta = trkPosVec.Eta();
+  phi = trkPosVec.Phi();
+  if(phi<0)
+    phi += 2*TMath::Pi();
+
+  return kTRUE;
+}
+
+
+//
+//------------------------------------------------------------------------------
+//
+Bool_t AliEMCALTracker::ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam, Float_t *clsPos, Double_t mass, Double_t step, Double_t &tmpEta, Double_t &tmpPhi)
+{
+  //
+  //Return the residual by extrapolating a track param to a global position
+  //
+  tmpEta = -999;
+  tmpPhi = -999;
+  if(!trkParam) return kFALSE;
+  Double_t trkPos[3] = {0.,0.,0.};
+  TVector3 vec(clsPos[0],clsPos[1],clsPos[2]);
+  Double_t alpha =  ((int)(vec.Phi()*TMath::RadToDeg()/20)+0.5)*20*TMath::DegToRad();
+  vec.RotateZ(-alpha); //Rotate the cluster to the local extrapolation coordinate system
+  if(!AliTrackerBase::PropagateTrackToBxByBz(trkParam, vec.X(), mass, step,kTRUE, 0.8, -1)) return kFALSE;
+  if(!trkParam->GetXYZ(trkPos)) return kFALSE; //Get the extrapolated global position
+
+  TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
+  TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+
+  // track cluster matching
+  tmpPhi = clsPosVec.DeltaPhi(trkPosVec);    // tmpPhi is between -pi and pi
+  tmpEta = clsPosVec.Eta()-trkPosVec.Eta();
+
+  return kTRUE;
+}
+
+
+//
+//------------------------------------------------------------------------------
+//
+Bool_t AliEMCALTracker::ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Double_t mass, Double_t step, Double_t &tmpEta, Double_t &tmpPhi)
+{
+  //
+  //Return the residual by extrapolating a track param to a cluster
+  //
+  tmpEta = -999;
+  tmpPhi = -999;
+  if(!cluster) return kFALSE;
+
+  Float_t clsPos[3] = {0.,0.,0.};
+  cluster->GetPosition(clsPos);
+
+  return ExtrapolateTrackToPosition(trkParam, clsPos, mass, step, tmpEta, tmpPhi);
+}
+
+
 //
 //------------------------------------------------------------------------------
 //
index 88321a6..739d71d 100644 (file)
@@ -29,8 +29,10 @@ class TList;
 class TTree;
 class TObjArray;
 class AliESDEvent;
+class AliVCluster;
 class AliESDCaloCluster;
 class AliEMCALTrack;
+class AliExternalTrackParam;
 class AliEMCALRecPoint;
 class AliEMCALGeometry;
 
@@ -63,6 +65,10 @@ public:
        void                SetStepLength(Float_t length) {fStep=length;}
        void                SetTrackCorrectionMode(Option_t *option);
 
+       static Bool_t ExtrapolateTrackToEMCalSurface(AliExternalTrackParam *trkParam, Double_t emcalR, Double_t mass, Double_t step, Double_t &eta, Double_t &phi);
+       static Bool_t ExtrapolateTrackToPosition(AliExternalTrackParam *trkParam, Float_t *clsPos, Double_t mass, Double_t step, Double_t &tmpEta, Double_t &tmpPhi);
+       static Bool_t ExtrapolateTrackToCluster(AliExternalTrackParam *trkParam, AliVCluster *cluster, Double_t mass, Double_t step, Double_t &tmpEta, Double_t &tmpPhi);
+
        enum {  kUnmatched = -99999 };
        
        class  AliEMCALMatchCluster : public TObject
@@ -84,7 +90,7 @@ public:
        };
    
 private:
-       Int_t FindMatchedCluster(AliESDtrack *track);
+       Int_t  FindMatchedCluster(AliESDtrack *track);
        
        enum ETrackCorr { 
                kTrackCorrNone  = 0, // do not correct for energy loss
@@ -96,7 +102,8 @@ private:
        Double_t    fCutNTPC;         // mimimum number of track hits in the TPC
        
        Float_t     fStep;            // Length of each step in propagation
-       ETrackCorr  fTrackCorrMode;   // Material budget correction mode        
+       ETrackCorr  fTrackCorrMode;   // Material budget correction mode
+       Double_t    fClusterWindow;   // Select clusters in the window to be matched to tracks
        Double_t    fCutEta;          // cut on eta difference
        Double_t    fCutPhi;          // cut on phi difference
 
@@ -105,7 +112,7 @@ private:
        
        AliEMCALGeometry *fGeom;      //! EMCAL geometry
        
-       ClassDef(AliEMCALTracker, 4)  // EMCAL "tracker"
+       ClassDef(AliEMCALTracker, 5)  // EMCAL "tracker"
 };
 
 #endif