STEER/AliTracker.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //               Implementation of the AliTracker class
  20 //  that is the base for AliTPCtracker, AliITStrackerV2 and AliTRDtracker
  21 //        Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  22 //-------------------------------------------------------------------------
  23 #include <TClass.h>
  24 #include <TMath.h>
  25 #include <TGeoManager.h>
  26
  27 #include "AliMagF.h"
  28 #include "AliTracker.h"
  29 #include "AliCluster.h"
  30 #include "AliKalmanTrack.h"
  31
  32 extern TGeoManager *gGeoManager;
  33
  34 Bool_t AliTracker::fgUniformField=kTRUE;
  35 Double_t AliTracker::fgBz=kAlmost0Field;
  36 const AliMagF *AliTracker::fgkFieldMap=0;
  37
  38 ClassImp(AliTracker)
  39
  40 AliTracker::AliTracker():
  41   TObject(),
  42   fX(0),
  43   fY(0),
  44   fZ(0),
  45   fSigmaX(0.005),
  46   fSigmaY(0.005),
  47   fSigmaZ(0.010)
  48 {
  49   //--------------------------------------------------------------------
  50   // The default constructor.
  51   //--------------------------------------------------------------------
  52   if (!fgkFieldMap) AliWarning("Field map is not set. Call AliTracker::SetFieldMap before creating a tracker!");
  53 }
  54
  55 //__________________________________________________________________________
  56 AliTracker::AliTracker(const AliTracker &atr):
  57   TObject(atr),
  58   fX(atr.fX),
  59   fY(atr.fY),
  60   fZ(atr.fZ),
  61   fSigmaX(atr.fSigmaX),
  62   fSigmaY(atr.fSigmaY),
  63   fSigmaZ(atr.fSigmaZ)
  64 {
  65   //--------------------------------------------------------------------
  66   // The default constructor.
  67   //--------------------------------------------------------------------
  68   if (!fgkFieldMap) AliWarning("Field map is not set. Call AliTracker::SetFieldMap before creating a tracker!");
  69 }
  70
  71 //__________________________________________________________________________
  72 void AliTracker::SetFieldMap(const AliMagF* map, Bool_t uni) {
  73   //--------------------------------------------------------------------
  74   //This passes the field map to the reconstruction.
  75   //--------------------------------------------------------------------
  76   if (map==0) AliFatalClass("Can't access the field map !");
  77
  78   if (fgkFieldMap) {
  79      AliWarningClass("The magnetic field map has been already set !");
  80      return;
  81   }
  82
  83   fgUniformField=uni;
  84   fgkFieldMap=map;
  85
  86   //Float_t r[3]={0.,0.,0.},b[3]; map->Field(r,b);
  87   //Double_t bz=-b[2];
  88
  89   Double_t bz=-map->SolenoidField();
  90   fgBz=TMath::Sign(kAlmost0Field,bz) + bz;
  91
  92 }
  93
  94 //__________________________________________________________________________
  95 void AliTracker::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
  96   //--------------------------------------------------------------------
  97   //This function "cooks" a track label. If label<0, this track is fake.
  98   //--------------------------------------------------------------------
  99   Int_t noc=t->GetNumberOfClusters();
 100   Int_t *lb=new Int_t[noc];
 101   Int_t *mx=new Int_t[noc];
 102   AliCluster **clusters=new AliCluster*[noc];
 103
 104   Int_t i;
 105   for (i=0; i<noc; i++) {
 106      lb[i]=mx[i]=0;
 107      Int_t index=t->GetClusterIndex(i);
 108      clusters[i]=GetCluster(index);
 109   }
 110
 111   Int_t lab=123456789;
 112   for (i=0; i<noc; i++) {
 113     AliCluster *c=clusters[i];
 114     lab=TMath::Abs(c->GetLabel(0));
 115     Int_t j;
 116     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
 117     lb[j]=lab;
 118     (mx[j])++;
 119   }
 120
 121   Int_t max=0;
 122   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
 123
 124   for (i=0; i<noc; i++) {
 125     AliCluster *c=clusters[i];
 126     //if (TMath::Abs(c->GetLabel(1)) == lab ||
 127     //    TMath::Abs(c->GetLabel(2)) == lab ) max++;
 128     if (TMath::Abs(c->GetLabel(0)!=lab))
 129         if (TMath::Abs(c->GetLabel(1)) == lab ||
 130             TMath::Abs(c->GetLabel(2)) == lab ) max++;
 131   }
 132
 133   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
 134   t->SetFakeRatio((1.- Float_t(max)/noc));
 135   t->SetLabel(lab);
 136
 137   delete[] lb;
 138   delete[] mx;
 139   delete[] clusters;
 140 }
 141
 142 //____________________________________________________________________________
 143 void AliTracker::UseClusters(const AliKalmanTrack *t, Int_t from) const {
 144   //------------------------------------------------------------------
 145   //This function marks clusters associated with the track.
 146   //------------------------------------------------------------------
 147   Int_t noc=t->GetNumberOfClusters();
 148   for (Int_t i=from; i<noc; i++) {
 149      Int_t index=t->GetClusterIndex(i);
 150      AliCluster *c=GetCluster(index);
 151      c->Use();
 152   }
 153 }
 154
 155 Double_t AliTracker::GetBz(Float_t *r) {
 156   //------------------------------------------------------------------
 157   // Returns Bz (kG) at the point "r" .
 158   //------------------------------------------------------------------
 159     Float_t b[3]; fgkFieldMap->Field(r,b);
 160     Double_t bz=-Double_t(b[2]);
 161     return  (TMath::Sign(kAlmost0Field,bz) + bz);
 162 }
 163
 164 Double_t
 165 AliTracker::MeanMaterialBudget(Double_t *start,Double_t *end,Double_t *mparam)
 166 {
 167   //
 168   // Calculate mean material budget and material properties between
 169   //    the points "start" and "end".
 170   //
 171   // "mparam" - parameters used for the energy and multiple scattering
 172   //  corrections:
 173   //
 174   // mparam[0] - mean density: sum(x_i*rho_i)/sum(x_i) [g/cm3]
 175   // mparam[1] - equivalent rad length fraction: sum(x_i/X0_i) [adimensional]
 176   // mparam[2] - mean A: sum(x_i*A_i)/sum(x_i) [adimensional]
 177   // mparam[3] - mean Z: sum(x_i*Z_i)/sum(x_i) [adimensional]
 178   // mparam[4] - length: sum(x_i) [cm]
 179   // mparam[5] - Z/A mean: sum(x_i*Z_i/A_i)/sum(x_i) [adimensional]
 180   // mparam[6] - number of boundary crosses
 181   //
 182   //  Origin:  Marian Ivanov, Marian.Ivanov@cern.ch
 183   //
 184   //  Corrections and improvements by
 185   //        Andrea Dainese, Andrea.Dainese@lnl.infn.it,
 186   //        Andrei Gheata,  Andrei.Gheata@cern.ch
 187   //
 188
 189   mparam[0]=0; mparam[1]=1; mparam[2] =0; mparam[3] =0;
 190   mparam[4]=0; mparam[5]=0; mparam[6]=0;
 191   //
 192   Double_t bparam[6]; // total parameters
 193   Double_t lparam[6]; // local parameters
 194
 195   for (Int_t i=0;i<6;i++) bparam[i]=0;
 196
 197   if (!gGeoManager) {
 198     printf("ERROR: no TGeo\n");
 199     return 0.;
 200   }
 201   //
 202   Double_t length;
 203   Double_t dir[3];
 204   length = TMath::Sqrt((end[0]-start[0])*(end[0]-start[0])+
 205                        (end[1]-start[1])*(end[1]-start[1])+
 206                        (end[2]-start[2])*(end[2]-start[2]));
 207   mparam[4]=length;
 208   if (length<TGeoShape::Tolerance()) return 0.0;
 209   Double_t invlen = 1./length;
 210   dir[0] = (end[0]-start[0])*invlen;
 211   dir[1] = (end[1]-start[1])*invlen;
 212   dir[2] = (end[2]-start[2])*invlen;
 213
 214   // Initialize start point and direction
 215   TGeoNode *currentnode = 0;
 216   TGeoNode *startnode = gGeoManager->InitTrack(start, dir);
 217   //printf("%s length=%f\n",gGeoManager->GetPath(),length);
 218   if (!startnode) {
 219     AliErrorClass(Form("start point out of geometry: x %f, y %f, z %f",
 220                   start[0],start[1],start[2]));
 221     return 0.0;
 222   }
 223   TGeoMaterial *material = startnode->GetVolume()->GetMedium()->GetMaterial();
 224   lparam[0]   = material->GetDensity();
 225   lparam[1]   = material->GetRadLen();
 226   lparam[2]   = material->GetA();
 227   lparam[3]   = material->GetZ();
 228   lparam[4]   = length;
 229   lparam[5]   = lparam[3]/lparam[2];
 230   if (material->IsMixture()) {
 231     TGeoMixture * mixture = (TGeoMixture*)material;
 232     lparam[5] =0;
 233     Double_t sum =0;
 234     for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 235       sum  += mixture->GetWmixt()[iel];
 236       lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 237     }
 238     lparam[5]/=sum;
 239   }
 240
 241   // Locate next boundary within length without computing safety.
 242   // Propagate either with length (if no boundary found) or just cross boundary
 243   gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
 244   Double_t step = 0.0; // Step made
 245   Double_t snext = gGeoManager->GetStep();
 246   // If no boundary within proposed length, return current density
 247   if (!gGeoManager->IsOnBoundary()) {
 248     mparam[0] = lparam[0];
 249     mparam[1] = lparam[4]/lparam[1];
 250     mparam[2] = lparam[2];
 251     mparam[3] = lparam[3];
 252     mparam[4] = lparam[4];
 253     return lparam[0];
 254   }
 255   // Try to cross the boundary and see what is next
 256   Int_t nzero = 0;
 257   while (length>TGeoShape::Tolerance()) {
 258     currentnode = gGeoManager->GetCurrentNode();
 259     if (snext<2.*TGeoShape::Tolerance()) nzero++;
 260     else nzero = 0;
 261     if (nzero>3) {
 262       // This means navigation has problems on one boundary
 263       // Try to cross by making a small step
 264       printf("ERROR: cannot cross boundary\n");
 265       mparam[0] = bparam[0]/step;
 266       mparam[1] = bparam[1];
 267       mparam[2] = bparam[2]/step;
 268       mparam[3] = bparam[3]/step;
 269       mparam[5] = bparam[5]/step;
 270       mparam[4] = step;
 271       mparam[0] = 0.;             // if crash of navigation take mean density 0
 272       mparam[1] = 1000000;        // and infinite rad length
 273       return bparam[0]/step;
 274     }
 275     mparam[6]+=1.;
 276     step += snext;
 277     bparam[1]    += snext/lparam[1];
 278     bparam[2]    += snext*lparam[2];
 279     bparam[3]    += snext*lparam[3];
 280     bparam[5]    += snext*lparam[5];
 281     bparam[0]    += snext*lparam[0];
 282
 283     if (snext>=length) break;
 284     if (!currentnode) break;
 285     length -= snext;
 286     //printf("%s snext=%f length=%f\n", currentnode->GetName(),snext,length);
 287     material = currentnode->GetVolume()->GetMedium()->GetMaterial();
 288     lparam[0] = material->GetDensity();
 289     lparam[1]  = material->GetRadLen();
 290     lparam[2]  = material->GetA();
 291     lparam[3]  = material->GetZ();
 292     //printf("       %f %f %f %f\n",lparam[0],lparam[1],lparam[2],lparam[3]);
 293     lparam[5]   = lparam[3]/lparam[2];
 294     if (material->IsMixture()) {
 295       TGeoMixture * mixture = (TGeoMixture*)material;
 296       lparam[5]=0;
 297       Double_t sum =0;
 298       for (Int_t iel=0;iel<mixture->GetNelements();iel++){
 299         sum+= mixture->GetWmixt()[iel];
 300         lparam[5]+= mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
 301       }
 302       lparam[5]/=sum;
 303     }
 304     gGeoManager->FindNextBoundaryAndStep(length, kFALSE);
 305     snext = gGeoManager->GetStep();
 306     //printf("snext %f\n",snext);
 307   }
 308   mparam[0] = bparam[0]/step;
 309   mparam[1] = bparam[1];
 310   mparam[2] = bparam[2]/step;
 311   mparam[3] = bparam[3]/step;
 312   mparam[5] = bparam[5]/step;
 313   return bparam[0]/step;
 314 }
 315
 316
 317 Bool_t
 318 AliTracker::PropagateTrackTo(AliExternalTrackParam *track, Double_t xToGo,
 319 Double_t mass, Double_t maxStep, Bool_t rotateTo, Double_t maxSnp){
 320   //----------------------------------------------------------------
 321   //
 322   // Propagates the track to the plane X=xk (cm) using the magnetic field map
 323   // and correcting for the crossed material.
 324   //
 325   // mass     - mass used in propagation - used for energy loss correction
 326   // maxStep  - maximal step for propagation
 327   //
 328   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 329   //
 330   //----------------------------------------------------------------
 331   const Double_t kEpsilon = 0.00001;
 332   Double_t xpos     = track->GetX();
 333   Double_t dir      = (xpos<xToGo) ? 1.:-1.;
 334   //
 335   while ( (xToGo-xpos)*dir > kEpsilon){
 336     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 337     Double_t x    = xpos+step;
 338     Double_t xyz0[3],xyz1[3],param[7];
 339     track->GetXYZ(xyz0);   //starting global position
 340
 341     Double_t bz=GetBz(xyz0); // getting the local Bz
 342
 343     if (!track->GetXYZAt(x,bz,xyz1)) return kFALSE;   // no prolongation
 344     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 345
 346     if (TMath::Abs(track->GetSnpAt(x,bz)) >= maxSnp) return kFALSE;
 347     if (!track->PropagateTo(x,bz))  return kFALSE;
 348
 349     MeanMaterialBudget(xyz0,xyz1,param);
 350     Double_t xrho=param[0]*param[4], xx0=param[1];
 351
 352     if (!track->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 353     if (rotateTo){
 354       if (TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 355       track->GetXYZ(xyz0);   // global position
 356       Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
 357       //
 358       Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
 359                sa=TMath::Sin(alphan-track->GetAlpha());
 360       Double_t sf=track->GetSnp(), cf=TMath::Sqrt(1.- sf*sf);
 361       Double_t sinNew =  sf*ca - cf*sa;
 362       if (TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 363       if (!track->Rotate(alphan)) return kFALSE;
 364     }
 365     xpos = track->GetX();
 366   }
 367   return kTRUE;
 368 }
 369