L1phase shift corrected
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTracker.cxx
index 23ddeff..15b34c2 100644 (file)
@@ -40,6 +40,7 @@
 #include <AliLog.h>
 #include "AliPHOSTracker.h"
 #include "AliPHOSEmcRecPoint.h"
+#include "AliPHOSGeometry.h"
 #include "AliESDEvent.h"
 #include "AliESDtrack.h"
 #include "AliPHOSTrackSegmentMakerv1.h"
@@ -55,17 +56,6 @@ ClassImp(AliPHOSTracker)
 
 Bool_t AliPHOSTracker::fgDebug = kFALSE ;  
 
-
-// ***** Some geometrical constants (used in PropagateBack) 
-
-const Double_t kR=460.+ 9;  // Radial coord. of the centre of EMC module (cm)
-
-const Double_t kAlpha=20.*TMath::Pi()/180.;     // Segmentation angle (rad)
-const Double_t kYmax=kR*TMath::Tan(0.5*kAlpha); // Maximal possible y-coord.(cm)
-const Double_t kZmax=65.; // Approximately: the maximal possible z-coord.(cm)
-
-
-
 //____________________________________________________________________________
 AliPHOSTracker::AliPHOSTracker(): 
   AliTracker()
@@ -95,10 +85,7 @@ Int_t AliPHOSTracker::LoadClusters(TTree *cTree) {
   //--------------------------------------------------------------------
   // This function loads the PHOS clusters
   //--------------------------------------------------------------------
-  return 0 ; //At this stage we can not strore result 
-             // the closest track and distance to it
-             //We perform same task later in AliPHOSTrackSegmentMakerv1
-/*
+
   TObjArray *arr=NULL;
   TBranch *branch=cTree->GetBranch("PHOSEmcRP");
   if (branch==0) {
@@ -107,6 +94,8 @@ Int_t AliPHOSTracker::LoadClusters(TTree *cTree) {
   }
   branch->SetAddress(&arr);
 
+  for(Int_t m=0;m<5; m++) fModules[m]->Clear("C") ;
+
   Int_t nclusters=0;
   Int_t nentr=(Int_t)branch->GetEntries();
   for (Int_t i=0; i<nentr; i++) {
@@ -117,8 +106,8 @@ Int_t AliPHOSTracker::LoadClusters(TTree *cTree) {
 
       Int_t m=cl->GetPHOSMod();
       if ((m<1)||(m>5)) {
-         AliError("Wrong module index !");
-         return 1;
+         AliError(Form("Wrong module index: %d !",m));
+         continue ;
       }
 
       // Here is how the alignment is treated
@@ -139,7 +128,7 @@ Int_t AliPHOSTracker::LoadClusters(TTree *cTree) {
   Info("LoadClusters","Number of loaded clusters: %d",nclusters);
 
   return 0;
-*/
+
 }
 
 //____________________________________________________________________________
@@ -150,10 +139,6 @@ Int_t AliPHOSTracker::PropagateBack(AliESDEvent *esd) {
   // Makes the PID
   //--------------------------------------------------------------------
 
-  return 0 ; //At this stage we can not strore result 
-             // the closest track and distance to it
-             //We perform same task later in AliPHOSTrackSegmentMakerv1
-/*
   Int_t nt=esd->GetNumberOfTracks();
 
   // *** Select and sort the ESD track in accordance with their quality
@@ -165,85 +150,88 @@ Int_t AliPHOSTracker::PropagateBack(AliESDEvent *esd) {
   }
   TMath::Sort(nt,quality,index,kFALSE);
 
+  AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
 
   // *** Start the matching
-  Double_t bz = GetGz() ; //For approximate matching
-  Double_t b[3]; GetBxByBz(b); //For final matching
+  TVector3 vecEmc ;   // Local position of EMC recpoint
+  Double_t bz = GetBz() ; //For approximate matching
+  Double_t b[3];  //For final matching
+  Double_t gposTrack[3] ;
   Int_t matched=0;
   for (Int_t i=0; i<nt; i++) {
      AliESDtrack *esdTrack=esd->GetTrack(index[i]);
 
-//     // Skip the tracks having "wrong" status (has to be checked/tuned)
-//     ULong_t status = esdTrack->GetStatus();
+     // Skip the tracks having "wrong" status (has to be checked/tuned)
+     ULong_t status = esdTrack->GetStatus();
+     if ((status & AliESDtrack::kTPCout)   == 0) continue;
 //     if ((status & AliESDtrack::kTRDout)   == 0) continue;
 //     if ((status & AliESDtrack::kTRDrefit) == 1) continue;
 
-     AliExternalTrackParam t(*esdTrack);
-
-     Int_t isec=Int_t(t.GetAlpha()/kAlpha);
-     Int_t imod=-isec-2; // PHOS module
-
-     Double_t y;                       // Some tracks do not reach the PHOS
-     if (!t.GetYAt(kR,bz,y)) continue; //    because of the bending
-
-     Double_t z; t.GetZAt(kR,bz,z);   
-     if (TMath::Abs(z) > kZmax) continue; // Some tracks miss the PHOS in Z
-     Bool_t ok=kTRUE;
-     while (TMath::Abs(y) > kYmax) {   // Find the matching module
-        Double_t alp=t.GetAlpha();
-        if (y > kYmax) {
-         if (!t.Rotate(alp+kAlpha)) {ok=kFALSE; break;}
-          imod--;
-        } else if (y < -kYmax) {
-         if (!t.Rotate(alp-kAlpha)) {ok=kFALSE; break;}
-          imod++;
-        }
-        if (!t.GetYAt(kR,bz,y)) {ok=kFALSE; break;}
-     }
-     if (!ok) continue; // Track rotation failed
+     //Continue extrapolation from TPC outer surface
+     const AliExternalTrackParam *outerParam=esdTrack->GetOuterParam();
+     if (!outerParam) continue;
+     AliExternalTrackParam t(*outerParam);
+
+     t.GetBxByBz(b) ;
+
+     //Loop over PHOS modules
+     Double_t dx=0,dz=0;
+     Double_t minDistance=999. ;
+     Int_t emcIndex=0  ;
+     for(Int_t mod=1; mod<=5; mod++){
+       if(fModules[mod-1]->GetEntriesFast()==0) //no PHOS clusters in this module or module does not exist
+         continue ;
+       
+       //Approximate direction to the current PHOS module
+       Double_t phiMod=(330.-20.*mod)/180.*TMath::Pi() ;
+       if(!t.Rotate(phiMod))
+         continue ;
+       TVector3 globaPos ;
+       geom->Local2Global(mod, 0.,0., globaPos) ;
+       const Double_t rPHOS = globaPos.Pt() ; //Distance to center of  PHOS module
+       const Double_t kYmax = 72.+10. ; //Size of the module (with some reserve) in phi direction
+       const Double_t kZmax = 64.+20. ; //Size of the module (with some reserve) in z direction
  
-  
-     if ((imod<0)||(imod>4)) continue; // Some tracks miss the PHOS in azimuth
-
-     //t.CorrectForMaterial(...); // Correct for the TOF material, if needed
-     t.PropagateToBxByBz(kR,b);        // Propagate to the matching module
-
-
-    // *** Search for the "best" cluster (can be improved)
-     TClonesArray &cArray=*fModules[imod];
-     Int_t ncl=cArray.GetEntriesFast();
-     AliPHOSEmcRecPoint *bestCluster=0;            // The "best" cluster
-     Double_t maxd2=400; // (cm^2)
-     for (Int_t j=0; j<ncl; j++) {
-       AliPHOSEmcRecPoint *c=(AliPHOSEmcRecPoint *)cArray.UncheckedAt(j);
-
-       //we looking at the closest track to the cluster, 
-       //not closest cluster to the track.
-//       if (c->TestBit(14)) continue; // This clusters is "used"
-
-       Double_t dy = t.GetY() - c->GetY(), dz = t.GetZ() - c->GetZ();
-       Double_t d2 = dy*dy + dz*dz;
-       if (d2 < maxd2) {
-         maxd2=d2;
-          bestCluster=c;
+       Double_t y;                       // Some tracks do not reach the PHOS
+       if (!t.GetYAt(rPHOS,bz,y)) continue; //    because of the bending
+
+       Double_t z;
+       if(!t.GetZAt(rPHOS,bz,z))
+         continue ;
+       if (TMath::Abs(z) > kZmax)
+         continue; // Some tracks miss the PHOS in Z
+       if(TMath::Abs(y) < kYmax){
+         t.PropagateToBxByBz(rPHOS,b);        // Propagate to the matching module
+      //t.CorrectForMaterial(...); // Correct for the TOF material, if needed
+         t.GetXYZ(gposTrack) ;
+         TVector3 globalPositionTr(gposTrack) ;
+         TVector3 localPositionTr ;
+         geom->Global2Local(localPositionTr,globalPositionTr,mod) ;
+         for(Int_t icl=0;icl<fModules[mod-1]->GetEntriesFast();icl++){
+           AliPHOSEmcRecPoint * clu =static_cast<AliPHOSEmcRecPoint*>(fModules[mod-1]->At(icl)) ;
+           clu->GetLocalPosition(vecEmc) ;
+           Double_t ddx = vecEmc.X()-localPositionTr.X();
+           Double_t ddz = vecEmc.Z()-localPositionTr.Z();
+           Double_t d2 = ddx*ddx + ddz*ddz;
+           if(d2 < minDistance) {
+             dx = ddx ;
+             dz = ddz ;
+             emcIndex=clu->GetIndexInList() ;
+             minDistance=d2 ;
+           }
+         }
+         break ;
        }
+     } //Loop over modules
+     if(minDistance<999.){
+       //found some match
+       esdTrack->SetStatus(AliESDtrack::kPHOSmatch) ;
+       esdTrack->SetPHOScluster(-emcIndex) ; //Should be ESDCaloCluster index which is not known yet. Will be transformed later in FillESD().
+       esdTrack->SetPHOSdxdz(dx,dz) ;
+       matched++;
      }
 
-     if (!bestCluster) continue;   // No reasonable matching found 
-
-//     bestCluster->SetBit(14,kTRUE); // This clusters is now attached to a track
-
-     matched++;
-
-     // *** Now, do the PID with the "bestCluster"
-     // and add the corresponding info to the ESD track pointed by "esdTrack"  
-
-     
-//     printf("%e %e %e %e\n",t.GetSign(), t.GetX() - bestCluster->GetX(),
-//                                      t.GetY() - bestCluster->GetY(),
-//                                      t.GetZ() - bestCluster->GetZ());
-     
   }
     
   Info("PropagateBack","Number of matched tracks: %d",matched);
@@ -252,7 +240,7 @@ Int_t AliPHOSTracker::PropagateBack(AliESDEvent *esd) {
   delete[] index;
 
   return 0;
-*/
+
 }
 
 //____________________________________________________________________________