Merge branch 'master' into LocalDev
[u/mrichter/AliRoot.git] / PHOS / AliPHOSTrackSegmentMakerv1.cxx
index 94c2a73..2dc7037 100644 (file)
 //
 
 // --- ROOT system ---
-#include "AliMagF.h"
-#include "AliTracker.h"
-#include "AliGeomManager.h"
-#include "AliCluster.h"
-#include "AliKalmanTrack.h"
-#include "AliGlobalQADataMaker.h"
-
 #include "TVector3.h"
 #include "TTree.h"
 #include "TBenchmark.h"
 #include "AliPHOSCpvRecPoint.h"
 #include "AliLog.h"
 #include "AliMagF.h"
+#include "AliMagF.h"
+#include "AliTracker.h"
+#include "AliGeomManager.h"
+#include "AliCluster.h"
+#include "AliKalmanTrack.h"
+#include "AliGlobalQADataMaker.h"
+#include "AliVParticle.h"
+
 
 ClassImp( AliPHOSTrackSegmentMakerv1) 
 
@@ -148,7 +149,7 @@ AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1() :
   fNTrackSegments(0),
   fRcpv(0.f),
   fRtpc(0.f),
-  fVtx(0.f), 
+  fVtx(0.f,0.f,0.f), 
   fLinkUpArray(0),
   fEmcFirst(0),
   fEmcLast(0),
@@ -169,7 +170,7 @@ AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(AliPHOSGeometry *geom) :
   fNTrackSegments(0),
   fRcpv(0.f),
   fRtpc(0.f),
-  fVtx(0.f), 
+  fVtx(0.f,0.f,0.f), 
   fLinkUpArray(0),
   fEmcFirst(0),
   fEmcLast(0),
@@ -192,7 +193,7 @@ AliPHOSTrackSegmentMakerv1::AliPHOSTrackSegmentMakerv1(const AliPHOSTrackSegment
   fNTrackSegments(0),
   fRcpv(0.f),
   fRtpc(0.f),
-  fVtx(0.f), 
+  fVtx(0.f,0.f,0.f), 
   fLinkUpArray(0),
   fEmcFirst(0),
   fEmcLast(0),
@@ -229,14 +230,14 @@ void  AliPHOSTrackSegmentMakerv1::FillOneModule()
   //First EMC clusters
   Int_t totalEmc = fEMCRecPoints->GetEntriesFast() ;
   for(fEmcFirst = fEmcLast; (fEmcLast < totalEmc) &&  
-       ((dynamic_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
+       ((static_cast<AliPHOSRecPoint *>(fEMCRecPoints->At(fEmcLast)))->GetPHOSMod() == fModule ); 
       fEmcLast ++)  ;
   
   //Now CPV clusters
   Int_t totalCpv = fCPVRecPoints->GetEntriesFast() ;
 
     for(fCpvFirst = fCpvLast; (fCpvLast < totalCpv) && 
-         ((dynamic_cast<AliPHOSRecPoint *>(fCPVRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule ); 
+         ((static_cast<AliPHOSRecPoint *>(fCPVRecPoints->At(fCpvLast)))->GetPHOSMod() == fModule ); 
        fCpvLast ++) ;
       
 }
@@ -255,12 +256,12 @@ void  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * em
 //  Float_t delta = 1 ;  // Width of the rows in sorting of RecPoints (in cm)
 //                       // if you change this value, change it as well in xxxRecPoint::Compare()
 
+  trackindex = -1;
+  dx         = 999.;
+  dz         = 999.;
+
   if(!cpvClu) {
     
-    trackindex = -1;
-    dx=999.;
-    dz=999.;
-    
     if(!emcClu) {
       return;
     }
@@ -271,69 +272,65 @@ void  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * em
     //Calculate actual distance to PHOS module
     TVector3 globaPos ;
     fGeom->Local2Global(iPHOSMod, 0.,0., globaPos) ;
-    const Double_t rPHOS = globaPos.Pt() ; //Distance to PHOS module
+    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.+10. ; //Size of the module (with some reserve) in z direction
-    const Double_t kAlpha= 20./180.*TMath::Pi() ; //TPC sector angular size
+    const Double_t kAlpha0=330./180.*TMath::Pi() ; //First PHOS module angular direction
+    const Double_t kAlpha= 20./180.*TMath::Pi() ; //PHOS module angular size
     Double_t minDistance = 1.e6;
 
-    TMatrixF gmat;
-    TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
-    emcClu->GetGlobalPosition(gposRecPoint,gmat);
-    Double_t gposTrack[3] ; 
+    TVector3 vecEmc ;   // Local position of EMC recpoint
+    emcClu->GetLocalPosition(vecEmc) ;
 
-    Double_t bz = GetBz() ; //B-Field for approximate matching
+    Double_t gposTrack[3] ; 
+    Double_t bz = AliTracker::GetBz() ; //B-Field for approximate matching
     Double_t b[3]; 
     for (Int_t i=0; i<nt; i++) {
       AliESDtrack *esdTrack=fESD->GetTrack(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
+      //Continue extrapolation from TPC outer surface
+      const AliExternalTrackParam *outerParam=esdTrack->GetOuterParam();
+      if (!outerParam) continue;
+      AliExternalTrackParam t(*outerParam);
 
+      t.GetBxByBz(b) ;
+      //Direction to the current PHOS module
+      Double_t phiMod=kAlpha0-kAlpha*iPHOSMod ;
+      if(!t.Rotate(phiMod))
+        continue ;
       Double_t y;                       // Some tracks do not reach the PHOS
       if (!t.GetYAt(rPHOS,bz,y)) continue; //    because of the bending
 
-      Double_t z; t.GetZAt(rPHOS,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(rPHOS,bz,y)) {ok=kFALSE; break;}
-      }
-      if (!ok) continue; // Track rotation failed
-
-
-      if(imod!= iPHOSMod-1) 
-       continue; //not even approximate coincidence
-
+      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.GetBxByBz(b) ;
-      t.PropagateToBxByBz(rPHOS,b);        // Propagate to the matching module
-      t.GetXYZ(gposTrack) ;
-      
-      Double_t ddx = gposTrack[0] - gposRecPoint.X(), ddy = gposTrack[1] - gposRecPoint.Y(), ddz = gposTrack[2] - gposRecPoint.Z();
-      Double_t d2 = ddx*ddx + ddy*ddy + ddz*ddz;
-      if(d2 < minDistance) {
-       dx = TMath::Sign(TMath::Sqrt(ddx*ddx + ddy*ddy),ddx) ;
-       dz = ddz ;
-       trackindex=i;
-       minDistance=d2 ;
+        t.GetXYZ(gposTrack) ;
+        TVector3 globalPositionTr(gposTrack) ;
+        TVector3 localPositionTr ;
+        fGeom->Global2Local(localPositionTr,globalPositionTr,iPHOSMod) ;
+        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 ;
+         trackindex=i;
+         minDistance=d2 ;
+        }
       }
-    }
+    } //Scanned all tracks
     return ;
   }
 
@@ -395,11 +392,12 @@ void  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * em
       
   }
     
-  if(trackindex>=0)
+  if(trackindex>=0) {
     AliDebug(1,Form("\t\tBest match for (xClu,zClu,eClu)=(%.3f,%.3f,%.3f): iTrack=%d, dR=%.3f",
                    locClu.X(),locClu.Z(),emcClu->GetEnergy(),
                    trackindex,TMath::Sqrt(dx*dx+dz*dz)));        
-  return;
+    return;
+  }
   
   Float_t distance2Track = fRtpc ; 
   
@@ -417,7 +415,7 @@ void  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * em
   
   emcClu->GetLocalPosition(vecEmc) ;
   
-  Double_t xCPV,zCPV ; //EMC-projected coordinates of CPV cluster 
+  Double_t xCPV=0,zCPV=0 ; //EMC-projected coordinates of CPV cluster 
   TVector3 cpvGlobal; // Global position of the CPV recpoint
   fGeom->GetGlobalPHOS((AliPHOSRecPoint*)cpvClu,cpvGlobal);
   Double_t vtxCPV[3]={cpvGlobal.X(),cpvGlobal.Y(),cpvGlobal.Z()} ;
@@ -437,14 +435,9 @@ void  AliPHOSTrackSegmentMakerv1::GetDistanceInPHOSPlane(AliPHOSEmcRecPoint * em
     Double_t rCPV = cpvGlobal.Pt() ;// Radius from IP to current point 
 
     // Extrapolate the global track direction if any to CPV and find the closest track
-    Int_t nTracks = fESD->GetNumberOfTracks();
     Int_t iClosestTrack = -1;
-    Double_t minDistance = 1.e6;
     TVector3 inPHOS ; 
 
-    AliESDtrack *track;
-    Double_t xyz[3] ;
-    Double_t pxyz[3]; 
     for (Int_t iTrack=0; iTrack<nTracks; iTrack++) {
       track = fESD->GetTrack(iTrack);
       if (!track->GetXYZAt(rCPV, fESD->GetMagneticField(), xyz))
@@ -555,20 +548,23 @@ void  AliPHOSTrackSegmentMakerv1::MakePairs()
   // remove them from the list of "unassigned". 
 
   //Make arrays to mark clusters already chosen
+  Int_t index;
+
   Int_t * emcExist = 0;
-  if(fEmcLast > fEmcFirst)
+  if(fEmcLast > fEmcFirst) {
     emcExist = new Int_t[fEmcLast-fEmcFirst] ;
-  
-  Int_t index;
-  for(index = 0; index <fEmcLast-fEmcFirst; index ++)
-    emcExist[index] = 1 ;
-  
+    for(index = 0; index <fEmcLast-fEmcFirst; index ++)
+      emcExist[index] = 1 ;
+  }
+  else 
+    return;
+
   Bool_t * cpvExist = 0;
-  if(fCpvLast > fCpvFirst)
+  if(fCpvLast > fCpvFirst) {
     cpvExist = new Bool_t[fCpvLast-fCpvFirst] ;
-  for(index = 0; index <fCpvLast-fCpvFirst; index ++)
-    cpvExist[index] = kTRUE ;
-  
+    for(index = 0; index <fCpvLast-fCpvFirst; index ++)
+      cpvExist[index] = kTRUE ;
+  }
   
   // Finds the smallest links and makes pairs of CPV and EMC clusters with smallest distance 
   TIter nextUp(fLinkUpArray) ;
@@ -581,22 +577,24 @@ void  AliPHOSTrackSegmentMakerv1::MakePairs()
 
     if(emcExist[linkUp->GetEmc()-fEmcFirst] != -1){
 
-      if(cpvExist[linkUp->GetCpv()-fCpvFirst]){ //CPV still exist
+      //array cpvExist[] should be non-zero as far as linkUp exists
+      //But Coverity requires to check it
+      if(cpvExist && cpvExist[linkUp->GetCpv()-fCpvFirst]){ //CPV still exist
          Float_t dx,dz ;
          linkUp->GetXZ(dx,dz) ;
         new ((* fTrackSegments)[fNTrackSegments]) 
-          AliPHOSTrackSegment(dynamic_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) , 
-                              dynamic_cast<AliPHOSCpvRecPoint *>(fCPVRecPoints->At(linkUp->GetCpv())) , 
+          AliPHOSTrackSegment(static_cast<AliPHOSEmcRecPoint *>(fEMCRecPoints->At(linkUp->GetEmc())) , 
+                              static_cast<AliPHOSCpvRecPoint *>(fCPVRecPoints->At(linkUp->GetCpv())) , 
                               linkUp->GetTrack(),dx,dz) ;
         
-       (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
+       (static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
        fNTrackSegments++ ;
        emcExist[linkUp->GetEmc()-fEmcFirst] = -1 ; //Mark emc  that Cpv was found 
        //mark CPV recpoint as already used 
        cpvExist[linkUp->GetCpv()-fCpvFirst] = kFALSE ;
       } //if CpvUp still exist
-    } 
-  }        
+    }
+  }
 
   //look through emc recPoints left without CPV
   if(emcExist){ //if there is emc rec point
@@ -611,13 +609,14 @@ void  AliPHOSTrackSegmentMakerv1::MakePairs()
          new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment(emcclu,nullpointer) ;
         else
           new ((*fTrackSegments)[fNTrackSegments]) AliPHOSTrackSegment(emcclu,0,track,dx,dz);
-       (dynamic_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
+       (static_cast<AliPHOSTrackSegment *>(fTrackSegments->At(fNTrackSegments)))->SetIndexInList(fNTrackSegments);
        fNTrackSegments++;    
       } 
     }
   }
   delete [] emcExist ; 
-  delete [] cpvExist ; 
+  if(cpvExist)
+    delete [] cpvExist ; 
 }
 
 //____________________________________________________________________________
@@ -697,36 +696,4 @@ void AliPHOSTrackSegmentMakerv1::PrintTrackSegments(Option_t * option)
     }  
   }
 }
-//__________________________________________________________________________
-Double_t AliPHOSTrackSegmentMakerv1::GetBz()const
-{
-  Double_t kAlmost0Field=1.e-13;
-  AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
-  if (!fld) return 0.5*kAlmost0Field;
-  Double_t bz = fld->SolenoidField();
-  return TMath::Sign(0.5*kAlmost0Field,bz) + bz;
-}
-//__________________________________________________________________________
-void AliPHOSTrackSegmentMakerv1::GetBxByBz(const Double_t r[3], Double_t b[3])const {
-  //------------------------------------------------------------------
-  // Returns Bx, By and Bz (kG) at the point "r" .
-  //------------------------------------------------------------------
-  Double_t kAlmost0Field=1.e-13;
-  AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
-  if (!fld) {
-     b[0] = b[1] = 0.;
-     b[2] = 0.5*kAlmost0Field;
-     return;
-  }
-
-  if (fld->IsUniform()) {
-     b[0] = b[1] = 0.;
-     b[2] = fld->SolenoidField();
-  }  else {
-     fld->Field(r,b);
-  }
-  b[2] = (TMath::Sign(0.5*kAlmost0Field,b[2]) + b[2]);
-  return;
-}
-