]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ITS/AliITStrackerMI.cxx
changes in AddTask pPb
[u/mrichter/AliRoot.git] / ITS / AliITStrackerMI.cxx
index f1484fb266600665e29d82a99491014ec368831d..59ddfa12f2a4bfad7ae79c9dda5fddaf71fd7acd 100644 (file)
@@ -707,7 +707,9 @@ Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
   fTrackingPhase="PropagateBack";
   Int_t nentr=event->GetNumberOfTracks();
   //  Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
-
+  double bz0 = GetBz();
+  const double kWatchStep=10.; // for larger steps watch arc vs segment difference
+  //
   Int_t ntrk=0;
   for (Int_t i=0; i<nentr; i++) {
      AliESDtrack *esd=event->GetTrack(i);
@@ -718,7 +720,22 @@ Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
      Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
      t.GetXYZ(xyzTrk); 
      Double_t dst2 = 0.;
-     for (Int_t icoord=0; icoord<3; icoord++) {Double_t di = xyzTrk[icoord] - xyzVtx[icoord];dst2 += di*di; } 
+     {
+       double dxs = xyzTrk[0] - xyzVtx[0];
+       double dys = xyzTrk[1] - xyzVtx[1];
+       double dzs = xyzTrk[2] - xyzVtx[2];
+       // RS: for large segment steps d use approximation of cicrular arc s by
+       // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
+       // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
+       // Hence s^2/d^2 = (1+1/6 p^2)^2
+       dst2 = dxs*dxs + dys*dys;
+       if (dst2 > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
+        double crv = TMath::Abs(esd->GetC(bz0));
+        double fcarc = 1.+crv*crv*dst2/6.;
+        dst2 *= fcarc*fcarc;
+       }
+       dst2 += dzs*dzs;
+     }
      t.StartTimeIntegral();
      t.AddTimeStep(TMath::Sqrt(dst2));
      //
@@ -3074,7 +3091,7 @@ Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID
     else{
       //
       Int_t tracks2[24], cluster[24];
-      for (Int_t i=0;i<trackindex;i++){ tracks2[i]=-1; cluster[i]=0;}
+      for (Int_t i=0;i<24;i++){ tracks2[i]=-1; cluster[i]=0;}
       Int_t index =0;
       //
       for (Int_t i=0;i<trackindex;i++){
@@ -3400,13 +3417,11 @@ void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
   // add track to the list of hypothesys
   //------------------------------------------------------------------
 
-  if (esdindex>=fTrackHypothesys.GetEntriesFast()) 
-    fTrackHypothesys.Expand(TMath::Max(fTrackHypothesys.GetSize(),esdindex*2+10));
   //
   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
   if (!array) {
     array = new TObjArray(10);
-    fTrackHypothesys.AddAt(array,esdindex);
+    fTrackHypothesys.AddAtAndExpand(array,esdindex);
   }
   array->AddLast(track);
 }
@@ -4788,6 +4803,7 @@ Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
     }
   }
 
+  if(zlocmin>zlocmax)return 0;
   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
   if (!nChipsInRoad) return 0;
@@ -4869,7 +4885,7 @@ Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
 }
 //------------------------------------------------------------------------
 //------------------------------------------------------------------------
-Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
+Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer){
 //
 // Method to be optimized further: 
 // Aim: decide whether a track can be used for PlaneEff evaluation
@@ -4905,19 +4921,19 @@ Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t
   AliITStrackMI tmp(*track);
 
 // require a minimal number of cluster in other layers and eventually clusters in closest layers 
-  Int_t ncl_out=0; Int_t ncl_in=0;
+  Int_t nclout=0; Int_t nclin=0;
   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
  AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
-   // if (tmp.GetClIndex(lay)>=0) ncl_out++;
-if(index[lay]>=0)ncl_out++;
+   // if (tmp.GetClIndex(lay)>=0) nclout++;
+if(index[lay]>=0)nclout++;
   }
   for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
-   if (index[lay]>=0) ncl_in++; 
+   if (index[lay]>=0) nclin++; 
   }
-  Int_t ncl=ncl_out+ncl_in;
+  Int_t ncl=nclout+nclin;
   Bool_t nextout = kFALSE;
   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
@@ -4925,7 +4941,7 @@ if(index[lay]>=0)ncl_out++;
   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
   // maximum number of missing clusters allowed in outermost layers
-  if(ncl_out<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff()) 
+  if(nclout<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff()) 
      return kFALSE; 
   // maximum number of missing clusters allowed (both in innermost and in outermost layers)
   if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff()) 
@@ -4960,9 +4976,20 @@ if(index[lay]>=0)ncl_out++;
   //
   Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
   Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
-  Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
-  Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
-                                                    // done for RecPoints
+  Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXFromBoundaryPlaneEff();
+  Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZFromBoundaryPlaneEff();
+  Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
+  // those are precisions in the tracking reference system
+  Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
+  // Use it also for the module reference system, as it is done for RecPoints
+
+  if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigFrmBndPlaneEff()){
+   if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) dx -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
+   else dx -= distx;
+
+   if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) dz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
+   else dz -= distz;
+  }
 
   // exclude tracks at boundary between detectors
   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
@@ -4974,6 +5001,28 @@ if(index[lay]>=0)ncl_out++;
        (locx+dx > blockXmx-boundaryWidth) ||
        (locz-dz < blockZmn+boundaryWidth) ||
        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
+
+  if(ilayer==0){
+       const AliESDEvent *myesd = ((AliESDtrack*)tmp.GetESDtrack())->GetESDEvent();
+       //The beam pipe
+       if (CorrectForPipeMaterial(&tmp,"inward")) {
+          const AliESDVertex* vtx = myesd->GetVertex();
+          if(!vtx) return kFALSE;
+          Double_t ddz[2],cov[3];
+          Double_t maxD=3.;
+          if(!tmp.PropagateToDCA(vtx,AliTracker::GetBz(),maxD,ddz,cov)) return kFALSE;
+          if(TMath::Abs(ddz[0])>=AliITSReconstructor::GetRecoParam()->GetDCACutPlaneEff()) return kFALSE;
+
+          Double_t covar[6]; vtx->GetCovMatrix(covar);
+          Double_t p[2]={tmp.GetParameter()[0]-ddz[0],tmp.GetParameter()[1]-ddz[1]};
+          Double_t c[3]={covar[2],0.,covar[5]};
+          Double_t chi2= ((AliESDtrack*)tmp.GetESDtrack())->GetPredictedChi2(p,c);
+          if (chi2>AliITSReconstructor::GetRecoParam()->GetVertexChi2CutPlaneEff()) return kFALSE;  // Use this to cut on chi^2
+
+       } else return kFALSE;
+    }
+
+
   return kTRUE;
 }
 //------------------------------------------------------------------------
@@ -5052,9 +5101,32 @@ void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilay
     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
   }
+
+  if(AliITSReconstructor::GetRecoParam()->GetSwitchOffStdSearchClusPlaneEff()){
+     Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXSearchClusterPlaneEff();
+     Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZSearchClusterPlaneEff();
+     Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXSearchClusterPlaneEff();
+     Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZSearchClusterPlaneEff();
+     msy = nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
+     msz = nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
+
+  if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigSrhClusPlaneEff()){
+     if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) msy -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
+     else msy -= distx;
+
+     if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) msz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
+     else msz -= distz;
+  }
+
+  msy *= msy;
+  msz *= msz;
+
+  }
+
+  if(msz==0 || msy==0){AliWarning("UseTrackForPlaneEff: null search frame"); return;}
+    
   msz = 1./msz; // 1/RoadZ^2
   msy = 1./msy; // 1/RoadY^2
-//
 
   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
   Int_t idetc=-1;
@@ -5130,7 +5202,7 @@ void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilay
     Int_t cltype[2]={-999,-999};
                                                           // and the module
 
-Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
+    Float_t angleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
 
     tr[0]=locx;
     tr[1]=locz;
@@ -5200,17 +5272,17 @@ Float_t AngleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute
     if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
     anglet *= 180./TMath::Pi();
 
-     AngleModTrack[2]=(Float_t) angle;
-     AngleModTrack[0]=(Float_t) anglet;
+     angleModTrack[2]=(Float_t) angle;
+     angleModTrack[0]=(Float_t) anglet;
      // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
-    AngleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
-    AngleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
-    AngleModTrack[1]*=180./TMath::Pi(); // in degree
+    angleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
+    angleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
+    angleModTrack[1]*=180./TMath::Pi(); // in degree
 
-    fPlaneEff->FillHistos(key,found,tr,clu,cltype,AngleModTrack);
+    fPlaneEff->FillHistos(key,found,tr,clu,cltype,angleModTrack);
 
     // For FO efficiency studies of SPD 
-    if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,AngleModTrack);
+    if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,angleModTrack);
   }
   if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
 return;
@@ -5261,4 +5333,3 @@ Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *stor
   }
   return nfound;
 }
-