New version of alignment framework.
authorcvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Feb 2006 10:21:57 +0000 (10:21 +0000)
committercvetan <cvetan@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 22 Feb 2006 10:21:57 +0000 (10:21 +0000)
18 files changed:
STEER/AliAlignObj.cxx
STEER/AliAlignObj.h
STEER/AliAlignmentTracks.cxx
STEER/AliAlignmentTracks.h
STEER/AliTrackFitter.cxx
STEER/AliTrackFitter.h
STEER/AliTrackFitterRieman.cxx
STEER/AliTrackFitterRieman.h
STEER/AliTrackPointArray.cxx
STEER/AliTrackPointArray.h
STEER/AliTrackResiduals.cxx
STEER/AliTrackResiduals.h
STEER/AliTrackResidualsChi2.cxx
STEER/AliTrackResidualsChi2.h
STEER/AliTrackResidualsFast.cxx [new file with mode: 0644]
STEER/AliTrackResidualsFast.h [new file with mode: 0644]
STEER/STEERLinkDef.h
STEER/libSTEER.pkg

index f414766..7b6d87e 100644 (file)
@@ -257,6 +257,17 @@ UShort_t AliAlignObj::LayerToVolUID(ELayerID layerId, Int_t modId)
 }
 
 //_____________________________________________________________________________
+UShort_t AliAlignObj::LayerToVolUID(Int_t   layerId, Int_t modId)
+{
+  // From detector (layer) index and module number (according to detector numbering)
+  // build fVolUID, unique numerical identity of that volume inside ALICE
+  // fVolUID is 16 bits, first 5 reserved for layerID (32 possible values),
+  // remaining 11 for module ID inside det (2048 possible values).
+  //
+  return ((UShort_t(layerId) << 11) | UShort_t(modId));
+}
+
+//_____________________________________________________________________________
 AliAlignObj::ELayerID AliAlignObj::VolUIDToLayer(UShort_t voluid, Int_t &modId)
 {
   // From detector (layer) name and module number (according to detector numbering)
index e4cb6a4..17db925 100644 (file)
@@ -68,6 +68,7 @@ class AliAlignObj : public TObject {
   static const char* LayerName(Int_t layer) { return fgLayerName[layer]; }
 
   static UShort_t LayerToVolUID(ELayerID layerId, Int_t modId);
+  static UShort_t LayerToVolUID(Int_t    layerId, Int_t modId);
   static ELayerID VolUIDToLayer(UShort_t voluid, Int_t &modId);
   static ELayerID VolUIDToLayer(UShort_t voluid);
 
index 089dece..f4f106e 100644 (file)
@@ -44,6 +44,7 @@ AliAlignmentTracks::AliAlignmentTracks():
   fLastIndex(0),
   fArrayIndex(0),
   fIsIndexBuilt(kFALSE),
+  fMisalignObjs(0),
   fTrackFitter(0),
   fMinimizer(0)
 {
@@ -61,6 +62,7 @@ AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
   fLastIndex(0),
   fArrayIndex(0),
   fIsIndexBuilt(kFALSE),
+  fMisalignObjs(0),
   fTrackFitter(0),
   fMinimizer(0)
 {
@@ -80,6 +82,7 @@ AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdt
   fLastIndex(0),
   fArrayIndex(0),
   fIsIndexBuilt(kFALSE),
+  fMisalignObjs(0),
   fTrackFitter(0),
   fMinimizer(0)
 {
@@ -222,10 +225,9 @@ void AliAlignmentTracks::BuildIndex()
 {
   // Build index of points tree entries
   // Used for access based on the volume IDs
-  if (fIsIndexBuilt)
-    ResetIndex();
-  else
-    fIsIndexBuilt = kTRUE;
+  if (fIsIndexBuilt) return;
+
+  fIsIndexBuilt = kTRUE;
 
   TFile *fPointsFile = TFile::Open(fPointsFilename);
   if (!fPointsFile || !fPointsFile->IsOpen()) {
@@ -235,17 +237,17 @@ void AliAlignmentTracks::BuildIndex()
   
   //  AliTrackPointArray* array = new AliTrackPointArray;
   AliTrackPointArray* array = 0;
-  TTree* pointsTree = (TTree*) fPointsFile->Get("spTree");
-  if (!pointsTree) {
+  fPointsTree = (TTree*) fPointsFile->Get("spTree");
+  if (!fPointsTree) {
     AliWarning("No pointsTree found!");
     return;
   }
-  pointsTree->SetBranchAddress("SP", &array);
+  fPointsTree->SetBranchAddress("SP", &array);
 
-  Int_t nArrays = pointsTree->GetEntries();
+  Int_t nArrays = fPointsTree->GetEntries();
   for (Int_t iArray = 0; iArray < nArrays; iArray++)
     {
-      pointsTree->GetEvent(iArray);
+      fPointsTree->GetEvent(iArray);
       if (!array) continue;
       for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
        UShort_t volId = array->GetVolumeID()[ipoint];
@@ -311,7 +313,10 @@ void AliAlignmentTracks::ResetIndex()
 {
   // Reset the value of the last filled index
   // Do not realocate memory
-  for (Int_t iLayer = AliAlignObj::kFirstLayer; iLayer < AliAlignObj::kLastLayer; iLayer++) {
+
+  fIsIndexBuilt = kFALSE;
+  
+  for (Int_t iLayer = 0; iLayer < AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer; iLayer++) {
     for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
       fLastIndex[iLayer][iModule] = 0;
     }
@@ -338,11 +343,11 @@ void AliAlignmentTracks::DeleteIndex()
 }
 
 //______________________________________________________________________________
-Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignobjfilename)
+Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignObjFileName, const char* arrayName)
 {
   // Read alignment object from a file
   // To be replaced by a call to CDB
-  AliWarning(Form("Method not yet implemented (%s) !",alignobjfilename));
+  AliWarning(Form("Method not yet implemented (%s in %s) !",arrayName,alignObjFileName));
 
   return kFALSE;
 }
@@ -355,8 +360,10 @@ void AliAlignmentTracks::InitAlignObjs()
   fAlignObjs = new AliAlignObj**[nLayers];
   for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
     fAlignObjs[iLayer] = new AliAlignObj*[AliAlignObj::LayerSize(iLayer)];
-    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
-      fAlignObjs[iLayer][iModule] = new AliAlignObjAngles;
+    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
+      UShort_t volid = AliAlignObj::LayerToVolUID(iLayer+ AliAlignObj::kFirstLayer,iModule);
+      fAlignObjs[iLayer][iModule] = new AliAlignObjAngles("",volid,0,0,0,0,0,0);
+    }
   }
 }
 
@@ -381,6 +388,7 @@ void AliAlignmentTracks::DeleteAlignObjs()
     delete [] fAlignObjs[iLayer];
   }
   delete [] fAlignObjs;
+  fAlignObjs = 0;
 }
 
 //______________________________________________________________________________
@@ -435,8 +443,9 @@ void AliAlignmentTracks::AlignLayer(AliAlignObj::ELayerID layer,
   // a given layer.
   // Tracks are fitted only within
   // the range defined by the user.
+  Int_t nModules = AliAlignObj::LayerSize(layer - AliAlignObj::kFirstLayer);
   while (iterations > 0) {
-    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(layer); iModule++) {
+    for (Int_t iModule = 0; iModule < nModules; iModule++) {
       UShort_t volId = AliAlignObj::LayerToVolUID(layer,iModule);
       AlignVolume(volId,layerRangeMin,layerRangeMax);
     }
@@ -445,39 +454,50 @@ void AliAlignmentTracks::AlignLayer(AliAlignObj::ELayerID layer,
 }
 
 //______________________________________________________________________________
-void AliAlignmentTracks::AlignVolume(UShort_t volid,
+void AliAlignmentTracks::AlignVolume(UShort_t volid, UShort_t volidfit,
                                     AliAlignObj::ELayerID layerRangeMin,
-                                    AliAlignObj::ELayerID layerRangeMax)
+                                    AliAlignObj::ELayerID layerRangeMax,
+                                    Int_t iterations)
 {
   // Align a single detector volume.
   // Tracks are fitted only within
-  // the range defined by the user.
+  // the range defined by the user
+  // (by layerRangeMin and layerRangeMax)
+  // or within the volid2
+  // Repeat the procedure 'iterations' times
 
   // First take the alignment object to be updated
   Int_t iModule;
   AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
-  AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
+  AliAlignObj *alignObj = fAlignObjs[iLayer-AliAlignObj::kFirstLayer][iModule];
 
   // Then load only the tracks with at least one
   // space point in the volume (volid)
   BuildIndex();
   AliTrackPointArray **points;
-  Int_t nArrays = LoadPoints(volid, points);
-  if (nArrays == 0) return;
+  // Start the iterations
+  while (iterations > 0) {
+    Int_t nArrays = LoadPoints(volid, points);
+    if (nArrays == 0) return;
+
+    AliTrackResiduals *minimizer = CreateMinimizer();
+    minimizer->SetNTracks(nArrays);
+    minimizer->SetAlignObj(alignObj);
+    AliTrackFitter *fitter = CreateFitter();
+    for (Int_t iArray = 0; iArray < nArrays; iArray++) {
+      fitter->SetTrackPointArray(points[iArray], kFALSE);
+      fitter->Fit(volid,volidfit,layerRangeMin,layerRangeMax);
+      AliTrackPointArray *pVolId,*pTrack;
+      fitter->GetTrackResiduals(pVolId,pTrack);
+      minimizer->AddTrackPointArrays(pVolId,pTrack);
+    }
+    minimizer->Minimize();
+    *alignObj *= *minimizer->GetAlignObj();
 
-  AliTrackResiduals *minimizer = CreateMinimizer();
-  minimizer->SetNTracks(nArrays);
-  minimizer->SetAlignObj(alignObj);
-  AliTrackFitter *fitter = CreateFitter();
-  for (Int_t iArray = 0; iArray < nArrays; iArray++) {
-    fitter->SetTrackPointArray(points[iArray], kFALSE);
-    AliTrackPointArray *pVolId = 0, *pTrack = 0;
-    fitter->Fit(volid,pVolId,pTrack,layerRangeMin,layerRangeMax);
-    minimizer->AddTrackPointArrays(pVolId,pTrack);
+    UnloadPoints(nArrays, points);
+    
+    iterations--;
   }
-  minimizer->Minimize();
-
-  UnloadPoints(nArrays, points);
 }
   
 //______________________________________________________________________________
@@ -488,36 +508,29 @@ Int_t AliAlignmentTracks::LoadPoints(UShort_t volid, AliTrackPointArray** &point
   // volume (volid).
   // Use the already created tree index for
   // fast access.
-  Int_t iModule;
-  AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
-  Int_t nArrays = fLastIndex[iLayer][iModule];
 
-  // In case of empty index
-  if (nArrays == 0) {
+  if (!fPointsTree) {
+    AliWarning("Tree with the space point arrays not initialized!");
     points = 0;
     return 0;
   }
 
-  if (!fPointsTree) {
-    AliWarning("Tree with the space point arrays not initialized!");
+  Int_t iModule;
+  AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
+  Int_t nArrays = fLastIndex[iLayer-AliAlignObj::kFirstLayer][iModule];
+
+  // In case of empty index
+  if (nArrays == 0) {
     points = 0;
     return 0;
   }
 
-  AliAlignObj *alignObj = fAlignObjs[iLayer][iModule];
-  TGeoHMatrix m;
-  alignObj->GetMatrix(m);
-  Double_t *rot = m.GetRotationMatrix();
-  Double_t *tr  = m.GetTranslation();
-
   AliTrackPointArray* array = 0;
   fPointsTree->SetBranchAddress("SP", &array);
 
   points = new AliTrackPointArray*[nArrays];
-  TArrayI *index = fArrayIndex[iLayer][iModule];
+  TArrayI *index = fArrayIndex[iLayer-AliAlignObj::kFirstLayer][iModule];
   AliTrackPoint p;
-  Float_t xyz[3],cov[6];
-  Float_t newxyz[3];
   for (Int_t iArray = 0; iArray < nArrays; iArray++) {
     fPointsTree->GetEvent((*index)[iArray]);
     if (!array) {
@@ -528,13 +541,21 @@ Int_t AliAlignmentTracks::LoadPoints(UShort_t volid, AliTrackPointArray** &point
     points[iArray] = new AliTrackPointArray(nPoints);
     for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
       array->GetPoint(p,iPoint);
-      p.GetXYZ(xyz,cov);
-      for (Int_t i = 0; i < 3; i++)
-       newxyz[i] = tr[i] 
-                 + xyz[0]*rot[3*i]
-                  + xyz[1]*rot[3*i+1]
-                 + xyz[2]*rot[3*i+2];
-      p.SetXYZ(newxyz,cov);
+      Int_t modnum;
+      AliAlignObj::ELayerID layer = AliAlignObj::VolUIDToLayer(p.GetVolumeID(),modnum);
+
+      // Misalignment is introduced here
+      // Switch it off in case of real
+      // alignment job!
+      if (fMisalignObjs) {
+       AliAlignObj *misalignObj = fMisalignObjs[layer-AliAlignObj::kFirstLayer][modnum];
+       if (misalignObj)
+         misalignObj->Transform(p);
+      }
+      // End of misalignment
+
+      AliAlignObj *alignObj = fAlignObjs[layer-AliAlignObj::kFirstLayer][modnum];
+      alignObj->Transform(p);
       points[iArray]->AddPoint(iPoint,&p);
     }
   }
@@ -575,3 +596,51 @@ AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
 
   return fMinimizer;
 }
+
+//______________________________________________________________________________
+Bool_t AliAlignmentTracks::Misalign(const char *misalignObjFileName, const char* arrayName)
+{
+  // The method reads from a file a set of AliAlignObj which are
+  // then used to apply misalignments directly on the track
+  // space-points. The method is supposed to be used only for
+  // fast development and debugging of the alignment algorithms.
+  // Be careful not to use it in the case of 'real' alignment
+  // scenario since it will bias the results.
+
+  // Initialize the misalignment objects array
+  Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
+  fMisalignObjs = new AliAlignObj**[nLayers];
+  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
+    fMisalignObjs[iLayer] = new AliAlignObj*[AliAlignObj::LayerSize(iLayer)];
+    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
+      fMisalignObjs[iLayer][iModule] = 0x0;
+  }
+
+  // Open the misliagnment file and load the array with
+  // misalignment objects
+  TFile* inFile = TFile::Open(misalignObjFileName,"READ");
+  if (!inFile || !inFile->IsOpen()) {
+    AliError(Form("Could not open misalignment file %s !",misalignObjFileName));
+    return kFALSE;
+  }
+
+  TClonesArray* array = ((TClonesArray*) inFile->Get(arrayName));
+  if (!array) {
+    AliError(Form("Could not find misalignment array %s in the file %s !",arrayName,misalignObjFileName));
+    inFile->Close();
+    return kFALSE;
+  }
+  inFile->Close();
+
+  // Store the misalignment objects for further usage  
+  Int_t nObjs = array->GetEntriesFast();
+  AliAlignObj::ELayerID layerId; // volume layer
+  Int_t modId; // volume ID inside the layer
+  for(Int_t i=0; i<nObjs; i++)
+    {
+      AliAlignObj* alObj = (AliAlignObj*)array->UncheckedAt(i);
+      alObj->GetVolUID(layerId,modId);
+      fMisalignObjs[layerId-AliAlignObj::kFirstLayer][modId] = alObj;
+    }
+  return kTRUE;
+}
index 5cf6a8a..44ad42c 100644 (file)
@@ -41,7 +41,7 @@ class AliAlignmentTracks : public TObject {
 /*   void BuildIndexLayer(AliAlignObj::ELayerID layer); */
 /*   void BuildIndexVolume(UShort_t volid); */
 
-  Bool_t ReadAlignObjs(const char *alignobjfilename = "AlignObjs.root");
+  Bool_t ReadAlignObjs(const char *alignObjFileName = "AlignObjs.root", const char* arrayName = "Alignment");
 
   void SetTrackFitter(AliTrackFitter *fitter) { fTrackFitter = fitter; }
   void SetMinimizer(AliTrackResiduals *minimizer) { fMinimizer = minimizer; }
@@ -51,9 +51,10 @@ class AliAlignmentTracks : public TObject {
                  AliAlignObj::ELayerID layerRangeMin = AliAlignObj::kFirstLayer,
                  AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer,
                  Int_t iterations = 1);
-  void AlignVolume(UShort_t volid,
+  void AlignVolume(UShort_t volid, UShort_t volidfit = 0,
                   AliAlignObj::ELayerID layerRangeMin = AliAlignObj::kFirstLayer,
-                  AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer);
+                  AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer,
+                  Int_t iterations = 1);
 
  protected:
 
@@ -71,6 +72,8 @@ class AliAlignmentTracks : public TObject {
   AliTrackFitter *CreateFitter();
   AliTrackResiduals *CreateMinimizer();
 
+  Bool_t Misalign(const char *misalignObjFileName, const char* arrayName);
+
   TChain           *fESDChain;       //! Chain with ESDs
   TString           fPointsFilename; //  Name of the file containing the track point arrays
   TFile            *fPointsFile;     //  File containing the track point arrays
@@ -79,6 +82,7 @@ class AliAlignmentTracks : public TObject {
   TArrayI        ***fArrayIndex;     //! Volume arrays which contains the tree index
   Bool_t            fIsIndexBuilt;   //  Is points tree index built
   AliAlignObj    ***fAlignObjs;      //  Array with alignment objects
+  AliAlignObj    ***fMisalignObjs;   //  Array with alignment objects used to introduce misalignment of the space-points
   AliTrackFitter   *fTrackFitter;    //  Pointer to the track fitter
   AliTrackResiduals*fMinimizer;      //  Pointer to track residuals minimizer
 
index 2385322..01b21e3 100644 (file)
@@ -34,6 +34,9 @@ AliTrackFitter::AliTrackFitter()
   for (Int_t i=0;i<6;i++) fParams[i] = 0;
   fCov = 0;
   fPoints = 0;
+  fPVolId = fPTrack = 0;
+  fChi2 = 0;
+  fNdf = 0;
   fIsOwner = kFALSE;
 }
 
@@ -44,6 +47,9 @@ AliTrackFitter::AliTrackFitter(AliTrackPointArray *array, Bool_t owner)
   //
   for (Int_t i=0;i<6;i++) fParams[i] = 0;
   fCov = new TMatrixDSym(6);
+  fPVolId = fPTrack = 0;
+  fChi2 = 0;
+  fNdf = 0;
   fIsOwner = kFALSE;
   SetTrackPointArray(array,owner);
 }
@@ -54,10 +60,12 @@ AliTrackFitter::AliTrackFitter(const AliTrackFitter &fitter):
 {
   // Copy constructor
   //
+  SetTrackPointArray(fitter.fPoints,fitter.fIsOwner);
   for (Int_t i=0;i<6;i++) fParams[i] = fitter.fParams[i];
   fCov = new TMatrixDSym(*fitter.fCov);
+  fChi2 = fitter.fChi2;
+  fNdf = fitter.fNdf;
   fIsOwner = kFALSE;
-  SetTrackPointArray(fitter.fPoints,fitter.fIsOwner);
 }
 
 //_____________________________________________________________________________
@@ -67,10 +75,12 @@ AliTrackFitter &AliTrackFitter::operator =(const AliTrackFitter& fitter)
   //
   if(this==&fitter) return *this;
 
+  SetTrackPointArray(fitter.fPoints);
   for (Int_t i=0;i<6;i++) fParams[i] = fitter.fParams[i];
   fCov = new TMatrixDSym(*fitter.fCov);
+  fChi2 = fitter.fChi2;
+  fNdf = fitter.fNdf;
   fIsOwner = kFALSE;
-  SetTrackPointArray(fitter.fPoints);
   
   return *this;
 }
@@ -89,6 +99,9 @@ void AliTrackFitter::Reset()
   for (Int_t i=0;i<6;i++) fParams[i] = 0;
   delete fCov;
   fCov = new TMatrixDSym(6);
+  fPVolId = fPTrack = 0;
+  fChi2 = 0;
+  fNdf = 0;
 }
 
 void AliTrackFitter::SetTrackPointArray(AliTrackPointArray *array, Bool_t owner)
index 5c700d6..f75bdf7 100644 (file)
@@ -27,8 +27,7 @@ class AliTrackFitter : public TObject {
 
   virtual void   Reset();
   virtual void   SetTrackPointArray(AliTrackPointArray *array, Bool_t owner = kTRUE);
-  virtual Bool_t Fit(UShort_t volId,
-                    AliTrackPointArray *pVolId, AliTrackPointArray *pTrack,
+  virtual Bool_t Fit(UShort_t volId,UShort_t volIdFit = 0,
                     AliAlignObj::ELayerID layerRangeMin = AliAlignObj::kFirstLayer,
                     AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer) = 0;
   virtual Bool_t GetPCA(const AliTrackPoint &pIn, AliTrackPoint &pOut) const = 0;
@@ -38,12 +37,21 @@ class AliTrackFitter : public TObject {
   const Float_t* GetZ() const {return fPoints->GetZ();}
   const Double_t* GetParam() const {return &fParams[0];}
   const TMatrixDSym &  GetCovariance() const {return *fCov;}
+  Float_t        GetChi2() const {return fChi2;}
+  Int_t          GetNdf()  const {return fNdf;}
+  Float_t        GetNormChi2() const { return (fNdf != 0) ? fChi2/fNdf : 0; }
+  void           GetTrackResiduals(AliTrackPointArray*& pVolId, AliTrackPointArray*& pTrack) const
+    { pVolId = fPVolId; pTrack = fPTrack; }
 
  protected:
 
   Double_t      fParams[6];    // Track parameters
-  TMatrixDSym   *fCov;         // Track cov matrix
+  TMatrixDSym  *fCov;          // Track cov matrix
   AliTrackPointArray *fPoints; // Pointer to the array with track space points
+  AliTrackPointArray *fPVolId; // Pointer to the array with space-points in volId
+  AliTrackPointArray *fPTrack; // Pointer to the array with track extrapolation points in volId
+  Float_t       fChi2;         // Chi squared of the fit
+  Int_t         fNdf;          // Number of degrees of freedom
   Bool_t  fIsOwner;            // Is the object owner of the space points array
 
  private:
index 793d77b..7fdc505 100644 (file)
@@ -12,7 +12,10 @@ AliTrackFitterRieman::AliTrackFitterRieman():
   //
   fAlpha = 0.;
   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
+  fSumYY = 0;
   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
+  fSumZZ = 0;
+  fNUsed = 0;
   fConv = kFALSE;
 }
 
@@ -25,7 +28,10 @@ AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t own
   //
   fAlpha = 0.;
   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
+  fSumYY = 0;
   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
+  fSumZZ = 0;
+  fNUsed = 0;
   fConv = kFALSE;
 }
 
@@ -37,7 +43,10 @@ AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
   //
   fAlpha = rieman.fAlpha;
   for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
+  fSumYY = rieman.fSumYY;
   for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
+  fSumZZ = rieman.fSumZZ;
+  fNUsed = rieman.fNUsed;
   fConv = rieman.fConv;
 }
 
@@ -51,7 +60,10 @@ AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRiema
 
   fAlpha = rieman.fAlpha;
   for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
+  fSumYY = rieman.fSumYY;
   for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
+  fSumZZ = rieman.fSumZZ;
+  fNUsed = rieman.fNUsed;
   fConv = rieman.fConv;
 
   return *this;
@@ -70,12 +82,14 @@ void AliTrackFitterRieman::Reset()
   AliTrackFitter::Reset();
   fAlpha = 0.;
   for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
+  fSumYY = 0;
   for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
+  fSumZZ = 0;
+  fNUsed = 0;
   fConv =kFALSE;
 }
 
-Bool_t AliTrackFitterRieman::Fit(UShort_t volId,
-                                AliTrackPointArray *pVolId, AliTrackPointArray *pTrack,
+Bool_t AliTrackFitterRieman::Fit(UShort_t volId,UShort_t volIdFit,
                                 AliAlignObj::ELayerID layerRangeMin,
                                 AliAlignObj::ELayerID layerRangeMax)
 {
@@ -95,20 +109,17 @@ Bool_t AliTrackFitterRieman::Fit(UShort_t volId,
   // and the space point; in Z - the track extrapolation on the Z
   // plane defined by the space point.
 
-  pVolId = pTrack = 0x0;
-  fConv = kFALSE;
+  Reset();
 
   Int_t npoints = fPoints->GetNPoints();
   if (npoints < 3) return kFALSE;
 
-  AliTrackPoint p;
+  AliTrackPoint p,plocal;
   fPoints->GetPoint(p,0);
   fAlpha = TMath::ATan2(p.GetY(),p.GetX());
-  Double_t sin = TMath::Sin(fAlpha);
-  Double_t cos = TMath::Cos(fAlpha);
 
   Int_t npVolId = 0;
-  Int_t npused = 0;
+  fNUsed = 0;
   Int_t *pindex = new Int_t[npoints];
   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
     {
@@ -118,15 +129,22 @@ Bool_t AliTrackFitterRieman::Fit(UShort_t volId,
        pindex[npVolId] = ipoint;
        npVolId++;
       }
-      if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
-         iVolId >= AliAlignObj::LayerToVolUID(layerRangeMax,0)) continue;
-      Float_t x = p.GetX()*cos + p.GetY()*sin;
-      Float_t y = p.GetY()*cos - p.GetX()*sin;
-      AddPoint(x,y,p.GetZ(),1,1);
-      npused++;
+      if (volIdFit != 0) {
+       if (iVolId != volIdFit) continue;
+      }
+      else {
+       if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
+           iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,
+                                               AliAlignObj::LayerSize(layerRangeMax-
+                                                                      AliAlignObj::kFirstLayer))) continue;
+      }
+      plocal = p.Rotate(fAlpha);
+      AddPoint(plocal.GetX(),plocal.GetY(),plocal.GetZ(),
+              TMath::Sqrt(plocal.GetCov()[3]),TMath::Sqrt(plocal.GetCov()[5]));
+      fNUsed++;
     }
 
-  if (npused < 3) {
+  if (fNUsed < 3) {
     delete [] pindex;
     return kFALSE;
   }
@@ -144,21 +162,59 @@ Bool_t AliTrackFitterRieman::Fit(UShort_t volId,
     return kFALSE;
   }
 
-  pVolId = new AliTrackPointArray(npVolId);
-  pTrack = new AliTrackPointArray(npVolId);
+  fPVolId = new AliTrackPointArray(npVolId);
+  fPTrack = new AliTrackPointArray(npVolId);
   AliTrackPoint p2;
   for (Int_t ipoint = 0; ipoint < npVolId; ipoint++)
     {
       Int_t index = pindex[ipoint];
       fPoints->GetPoint(p,index);
       if (GetPCA(p,p2)) {
-       pVolId->AddPoint(ipoint,&p);
-       pTrack->AddPoint(ipoint,&p2);
+       fPVolId->AddPoint(ipoint,&p);
+       fPTrack->AddPoint(ipoint,&p2);
       }
     }  
 
   delete [] pindex;
 
+  // debug info
+//   Float_t chi2 = 0, chi22 = 0;
+//   for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
+//     {
+//       fPoints->GetPoint(p,ipoint);
+//       UShort_t iVolId = p.GetVolumeID();
+//       if (volIdFit != 0) {
+//     if (iVolId != volIdFit) continue;
+//       }
+//       else {
+//     if (iVolId < AliAlignObj::LayerToVolUID(layerRangeMin,0) ||
+//         iVolId > AliAlignObj::LayerToVolUID(layerRangeMax,AliAlignObj::LayerSize(layerRangeMax-
+//                                                                                  AliAlignObj::kFirstLayer))) continue;
+//       }
+//       plocal = p.Rotate(fAlpha);
+//       Float_t delta = (fParams[0]*(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY())+
+//                    2.*plocal.GetX()*fParams[1]+
+//                    fParams[2]-
+//                    2.*plocal.GetY())/
+//                   (2.*TMath::Sqrt(plocal.GetCov()[3]));
+// //       Float_t delta2 = (fParams[3]+
+// //                 plocal.GetX()*fParams[4]+
+// //                 plocal.GetX()*plocal.GetX()*fParams[5]-
+// //                 plocal.GetZ())/
+// //                (TMath::Sqrt(plocal.GetCov()[5]));
+//       Double_t r = TMath::Sqrt(plocal.GetX()*plocal.GetX()+plocal.GetY()*plocal.GetY());
+//       Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
+//       Float_t delta2 = (fParams[3]+
+//                     r*fParams[4]+r*r*r*fParams[4]*Rm1*Rm1/24-
+//                     plocal.GetZ())/
+//                    (TMath::Sqrt(plocal.GetCov()[5]));
+//       chi2 += delta*delta;
+//       chi22 += delta2*delta2;
+//       //      printf("pulls %d %d %f %f\n",ipoint,iVolId,delta,delta2);
+      
+//     }
+//   printf("My chi2 = %f + %f = %f\n",chi2,chi22,chi2+chi22);
+
   return kTRUE;
 }
 
@@ -174,7 +230,7 @@ void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy,
   //
   //  (x^2+y^2 -2*x*x0 - 2*y*y0+ x0^2 -y0^2 -R^2 =0;  ==>
   //
-  //   substitution t = 1/(x^2+y^2),   u = 2*x*t, y = 2*y*t,  D0 = R^2 - x0^2- y0^2
+  //   substitution t = 1/(x^2+y^2),   u = 2*x*t, v = 2*y*t,  D0 = R^2 - x0^2- y0^2
   //
   //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
   //     
@@ -205,13 +261,27 @@ void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy,
   fSumXY[1] +=u*weight;      fSumXY[2]+=v*weight;  fSumXY[3]+=t*weight;
   fSumXY[4] +=u*u*weight;    fSumXY[5]+=t*t*weight;
   fSumXY[6] +=u*v*weight;    fSumXY[7]+=u*t*weight; fSumXY[8]+=v*t*weight;
+  fSumYY += v*v*weight;
   //
   // XZ part
   //
-  weight = 1./sz;
-  fSumXZ[0] +=weight;
-  fSumXZ[1] +=weight*x;   fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
-  fSumXZ[5] +=weight*z;   fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
+  if (0) {
+    weight = 1./(sz*sz);
+    fSumXZ[0] +=weight;
+    fSumXZ[1] +=weight*x;   fSumXZ[2] +=weight*x*x; fSumXZ[3] +=weight*x*x*x; fSumXZ[4] += weight*x*x*x*x;
+    fSumXZ[5] +=weight*z;   fSumXZ[6] +=weight*x*z; fSumXZ[7] +=weight*x*x*z;
+    fSumZZ += z*z*weight;
+  }
+  else {
+    weight = 1./(sz*sz);
+    fSumXZ[0] +=weight;
+    Double_t r = TMath::Sqrt(x*x+y*y);
+    fSumXZ[1] +=weight*r;   fSumXZ[2] +=weight*r*r; fSumXZ[3] +=weight*z; fSumXZ[4] += weight*r*z;
+    // Now the auxulary sums
+    fSumXZ[5] +=weight*r*r*r/24; fSumXZ[6] +=weight*r*r*r*r/12; fSumXZ[7] +=weight*r*r*r*z/24;
+    fSumXZ[8] +=weight*r*r*r*r*r*r/(24*24);
+    fSumZZ += z*z*weight;
+  }
 }
 
 void AliTrackFitterRieman::Update(){
@@ -220,6 +290,8 @@ void AliTrackFitterRieman::Update(){
   //
   //
   for (Int_t i=0;i<6;i++)fParams[i]=0;
+  fChi2 = 0;
+  fNdf = 0;
   Int_t conv=0;
   //
   // XY part
@@ -244,29 +316,65 @@ void AliTrackFitterRieman::Update(){
     fParams[0] = res(0,0);
     fParams[1] = res(0,1);
     fParams[2] = res(0,2);
+    TMatrixD  tmp = res*sums.T();
+    fChi2 += fSumYY - tmp(0,0);
+    fNdf  += fNUsed - 3;
     conv++;
   }
   //
   // XZ part
   //
-  TMatrixDSym     smatrixz(3);
-  smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2];
-  smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3];
-  smatrixz(2,2) = fSumXZ[4];
-  smatrixz.Invert();
-  if (smatrixz.IsValid()){
-    sums(0,0)    = fSumXZ[5];
-    sums(0,1)    = fSumXZ[6];
-    sums(0,2)    = fSumXZ[7];
-    TMatrixD res = sums*smatrixz;
-    fParams[3] = res(0,0);
-    fParams[4] = res(0,1);
-    fParams[5] = res(0,2);
-    for (Int_t i=0;i<3;i++)
-      for (Int_t j=0;j<=i;j++){
-       (*fCov)(i+2,j+2)=smatrixz(i,j);
+  if (0) {
+    TMatrixDSym     smatrixz(3);
+    smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(0,2) = fSumXZ[2];
+    smatrixz(1,1) = fSumXZ[2]; smatrixz(1,2) = fSumXZ[3];
+    smatrixz(2,2) = fSumXZ[4];
+    smatrixz.Invert();
+    TMatrixD        sumsxz(1,3);
+    if (smatrixz.IsValid()){
+      sumsxz(0,0)    = fSumXZ[5];
+      sumsxz(0,1)    = fSumXZ[6];
+      sumsxz(0,2)    = fSumXZ[7];
+      TMatrixD res = sumsxz*smatrixz;
+      fParams[3] = res(0,0);
+      fParams[4] = res(0,1);
+      fParams[5] = res(0,2);
+      for (Int_t i=0;i<3;i++)
+       for (Int_t j=0;j<=i;j++){
+         (*fCov)(i+2,j+2)=smatrixz(i,j);
+       }
+      TMatrixD  tmp = res*sumsxz.T();
+      fChi2 += fSumZZ - tmp(0,0);
+      fNdf  += fNUsed - 3;
+      conv++;
+    }
+  }
+  else {
+    Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
+    fSumXZ[1] += fSumXZ[5]*Rm1*Rm1;
+    fSumXZ[2] += fSumXZ[6]*Rm1*Rm1 + fSumXZ[8]*Rm1*Rm1*Rm1*Rm1;
+    fSumXZ[4] += fSumXZ[7]*Rm1*Rm1;
+
+    TMatrixDSym     smatrixz(2);
+    smatrixz(0,0) = fSumXZ[0]; smatrixz(0,1) = fSumXZ[1]; smatrixz(1,1) = fSumXZ[2];
+    smatrixz.Invert();
+    TMatrixD        sumsxz(1,2);
+    if (smatrixz.IsValid()){
+      sumsxz(0,0)    = fSumXZ[3];
+      sumsxz(0,1)    = fSumXZ[4];
+      TMatrixD res = sumsxz*smatrixz;
+      fParams[3] = res(0,0);
+      fParams[4] = res(0,1);
+      fParams[5] = 0;
+      for (Int_t i=0;i<2;i++)
+       for (Int_t j=0;j<=i;j++){
+         (*fCov)(i+3,j+3)=smatrixz(i,j);
+       }
+      TMatrixD  tmp = res*sumsxz.T();
+      fChi2 += fSumZZ - tmp(0,0);
+      fNdf  += fNUsed - 2;
+      conv++;
     }
-    conv++;
   }
 
   //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
@@ -391,8 +499,10 @@ Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) c
   Double_t yprime = TMath::Abs(R)*(y-y0)/dR + y0;
 
   // Now Z coordinate
-  Double_t phi  = TMath::ASin((x-x0)/R);
-  Double_t phi0 = TMath::ASin((0.-x0)/R);
+  Double_t x2 = xprime*cos + yprime*sin;
+  Double_t x02 = -fParams[1]/fParams[0];
+  Double_t phi  = TMath::ASin((x2-x02)/R);
+  Double_t phi0 = TMath::ASin((0.-x02)/R);
   Double_t dphi = (phi-phi0);
   Double_t zprime = fParams[3]+fParams[4]*dphi*R;
 
index 1c9c1fb..207756e 100644 (file)
@@ -13,8 +13,7 @@ class AliTrackFitterRieman : public AliTrackFitter{
   AliTrackFitterRieman &operator =(const AliTrackFitterRieman& rieman);
   virtual ~AliTrackFitterRieman();
 
-  Bool_t Fit(UShort_t volId,
-            AliTrackPointArray *pVolId, AliTrackPointArray *pTrack,
+  Bool_t Fit(UShort_t volId,UShort_t volIdFit = 0,
             AliAlignObj::ELayerID layerRangeMin = AliAlignObj::kFirstLayer,
             AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer);
   Bool_t GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const;
@@ -34,7 +33,10 @@ class AliTrackFitterRieman : public AliTrackFitter{
 
   Double_t      fAlpha;     //angle to transform to the fitting coordinate system
   Double_t      fSumXY[9];  //sums for XY part
+  Double_t      fSumYY;     //sum for YY part
   Double_t      fSumXZ[9];  //sums for XZ part
+  Double_t      fSumZZ;     //sum for ZZ part
+  Int_t         fNUsed;     //actual number of space-points used in the fit
   Bool_t        fConv;      // indicates convergation
 
  private:
index 10b0380..26106f1 100644 (file)
@@ -22,6 +22,9 @@
 //   cvetan.cheshkov@cern.ch 3/11/2005                                      //
 //////////////////////////////////////////////////////////////////////////////
 
+#include <TMath.h>
+#include <TMatrixDSym.h>
+
 #include "AliTrackPointArray.h"
 
 ClassImp(AliTrackPointArray)
@@ -228,3 +231,139 @@ void AliTrackPoint::GetXYZ(Float_t *xyz, Float_t *cov) const
   if (cov)
     memcpy(cov,fCov,6*sizeof(Float_t));
 }
+
+//______________________________________________________________________________
+Float_t AliTrackPoint::GetResidual(const AliTrackPoint &p, Bool_t weighted) const
+{
+  // This method calculates the track to space-point residuals. The track
+  // interpolation is also stored as AliTrackPoint. Using the option
+  // 'weighted' one can calculate the residual either with or without
+  // taking into account the covariance matrix of the space-point and
+  // track interpolation. The second case the residual becomes a pull.
+
+  Float_t res = 0;
+
+  if (!weighted) {
+    Float_t xyz[3],xyzp[3];
+    GetXYZ(xyz);
+    p.GetXYZ(xyzp);
+    res = (xyz[0]-xyzp[0])*(xyz[0]-xyzp[0])+
+          (xyz[1]-xyzp[1])*(xyz[1]-xyzp[1])+
+          (xyz[2]-xyzp[2])*(xyz[2]-xyzp[2]);
+  }
+  else {
+    Float_t xyz[3],xyzp[3];
+    Float_t cov[6],covp[6];
+    GetXYZ(xyz,cov);
+    TMatrixDSym mcov(3);
+    mcov(0,0) = cov[0]; mcov(0,1) = cov[1]; mcov(0,2) = cov[2];
+    mcov(1,0) = cov[1]; mcov(1,1) = cov[3]; mcov(1,2) = cov[4];
+    mcov(2,0) = cov[2]; mcov(2,1) = cov[4]; mcov(2,2) = cov[5];
+    p.GetXYZ(xyzp,covp);
+    TMatrixDSym mcovp(3);
+    mcovp(0,0) = covp[0]; mcovp(0,1) = covp[1]; mcovp(0,2) = covp[2];
+    mcovp(1,0) = covp[1]; mcovp(1,1) = covp[3]; mcovp(1,2) = covp[4];
+    mcovp(2,0) = covp[2]; mcovp(2,1) = covp[4]; mcovp(2,2) = covp[5];
+    TMatrixDSym msum = mcov + mcovp;
+    msum.Invert();
+    if (msum.IsValid()) {
+      for (Int_t i = 0; i < 3; i++)
+       for (Int_t j = 0; j < 3; j++)
+         res += (xyz[i]-xyzp[i])*(xyz[j]-xyzp[j])*msum(i,j);
+    }
+  }
+
+  return res;
+}
+
+//______________________________________________________________________________
+Float_t AliTrackPoint::GetAngle() const
+{
+  // The method uses the covariance matrix of
+  // the space-point in order to extract the
+  // orientation of the detector plane.
+  // The rotation in XY plane only is calculated.
+
+  if ((fCov[2] != 0) || (fCov[4] != 0))
+    return TMath::ATan2(-fCov[4],-fCov[2]);
+  else {
+    if (fCov[1] == 0) return 0;
+    Float_t phi= TMath::ATan2(TMath::Abs(fCov[1]),fCov[3]);
+    if (fCov[1] < 0) {
+      phi += TMath::PiOver2();
+      if ((fY-fX) < 0) phi += TMath::Pi();
+    }
+    else {
+      if ((fX+fY) < 0) phi += TMath::Pi();
+    }
+    return phi;
+  } 
+//     return TMath::ATan2(TMath::Sign((Float_t)1.0,fY)*TMath::Abs(fCov[1]),
+//                     TMath::Sign((Float_t)1.0,fX)*fCov[3]);
+}
+
+//_____________________________________________________________________________
+AliTrackPoint& AliTrackPoint::Rotate(Float_t alpha) const
+{
+  // Transform the space-point coordinates
+  // and covariance matrix from global to
+  // local (detector plane) coordinate system
+  // XY plane rotation only
+
+  static AliTrackPoint p;
+  p = *this;
+
+  Float_t xyz[3],cov[6];
+  GetXYZ(xyz,cov);
+
+  Float_t sin = TMath::Sin(alpha), cos = TMath::Cos(alpha);
+
+  Float_t newxyz[3],newcov[6];
+  newxyz[0] = cos*xyz[0] + sin*xyz[1];
+  newxyz[1] = cos*xyz[1] - sin*xyz[0];
+  newxyz[2] = xyz[2];
+
+  newcov[0] = cov[0]*cos*cos+
+            2*cov[1]*sin*cos+
+              cov[3]*sin*sin;
+  newcov[1] = cov[1]*(cos*cos-sin*sin)+
+             (cov[3]-cov[0])*sin*cos;
+  newcov[2] = cov[2]*cos+
+              cov[4]*sin;
+  newcov[3] = cov[0]*sin*sin-
+            2*cov[1]*sin*cos+
+              cov[3]*cos*cos;
+  newcov[4] = cov[4]*cos-
+              cov[2]*sin;
+  newcov[5] = cov[5];
+
+  p.SetXYZ(newxyz,newcov);
+  p.SetVolumeID(GetVolumeID());
+
+  return p;
+}
+
+//_____________________________________________________________________________
+AliTrackPoint& AliTrackPoint::MasterToLocal() const
+{
+  // Transform the space-point coordinates
+  // and the covariance matrix from the
+  // (master) to the local (tracking)
+  // coordinate system
+
+  Float_t alpha = GetAngle();
+  return Rotate(alpha);
+}
+
+//_____________________________________________________________________________
+void AliTrackPoint::Print(Option_t *) const
+{
+  // Print the space-point coordinates and
+  // covariance matrix
+
+  printf("VolumeID=%d\n", GetVolumeID());
+  printf("X = %12.6f    Tx = %12.6f%12.6f%12.6f\n", fX, fCov[0], fCov[1], fCov[2]);
+  printf("Y = %12.6f    Ty = %12.6f%12.6f%12.6f\n", fY, fCov[1], fCov[3], fCov[4]);
+  printf("Z = %12.6f    Tz = %12.6f%12.6f%12.6f\n", fZ, fCov[2], fCov[4], fCov[5]);
+
+}
index 7493625..b7c6aa4 100644 (file)
@@ -37,6 +37,14 @@ class AliTrackPoint : public TObject {
   const Float_t *GetCov() const { return &fCov[0]; }
   UShort_t GetVolumeID() const { return fVolumeID; }
 
+  Float_t  GetResidual(const AliTrackPoint &p, Bool_t weighted = kFALSE) const;
+
+  Float_t  GetAngle() const;
+  AliTrackPoint& Rotate(Float_t alpha) const;
+  AliTrackPoint& MasterToLocal() const;
+
+  void     Print(Option_t *) const;
+
  private:
 
   Float_t  fX;        // X coordinate
@@ -81,7 +89,7 @@ class AliTrackPointArray : public TObject {
   const UShort_t* GetVolumeID() const { return &fVolumeID[0]; }
 
   Bool_t    HasVolumeID(UShort_t volid) const;
-  
+
  private:
 
   Int_t     fNPoints;    // Number of space points
index 4604fd6..157d620 100644 (file)
@@ -30,18 +30,23 @@ ClassImp(AliTrackResiduals)
 AliTrackResiduals::AliTrackResiduals():
   fN(0),
   fLast(0),
+  fAlignObj(0),
+  fVolArray(0),
+  fTrackArray(0),
+  fChi2(0),
+  fNdf(0),
   fIsOwner(kTRUE)
 {
   // Default constructor
-  fAlignObj = 0x0;
-  fVolArray = fTrackArray = 0x0;
 }
 
 //_____________________________________________________________________________
-AliTrackResiduals::AliTrackResiduals(Int_t ntracks, AliAlignObj *alignobj):
+AliTrackResiduals::AliTrackResiduals(Int_t ntracks):
   fN(ntracks),
   fLast(0),
-  fAlignObj(alignobj),
+  fAlignObj(0),
+  fChi2(0),
+  fNdf(0),
   fIsOwner(kTRUE)
 {
   // Constructor
@@ -58,11 +63,15 @@ AliTrackResiduals::AliTrackResiduals(const AliTrackResiduals &res):
   TObject(res),
   fN(res.fN),
   fLast(res.fLast),
-  fAlignObj(res.fAlignObj),
+  fChi2(res.fChi2),
+  fNdf(res.fNdf),
   fIsOwner(kTRUE)
 {
   // Copy constructor
   // By default the created copy owns the track point arrays
+  if (res.fAlignObj)
+    fAlignObj = (AliAlignObj *)res.fAlignObj->Clone();
+
   if (fN > 0) {
     fVolArray = new AliTrackPointArray*[fN];
     fTrackArray = new AliTrackPointArray*[fN];
@@ -90,6 +99,8 @@ AliTrackResiduals &AliTrackResiduals::operator =(const AliTrackResiduals& res)
 
   fN = res.fN;
   fLast = res.fLast;
+  fChi2 = res.fChi2;
+  fNdf  = res.fNdf;
   fIsOwner = kFALSE;
   fAlignObj = res.fAlignObj;
 
@@ -103,6 +114,7 @@ AliTrackResiduals &AliTrackResiduals::operator =(const AliTrackResiduals& res)
 AliTrackResiduals::~AliTrackResiduals()
 {
   // Destructor
+  if (fAlignObj) delete fAlignObj;
   DeleteTrackPointArrays();
 }
 
@@ -116,6 +128,8 @@ void AliTrackResiduals::SetNTracks(Int_t ntracks)
 
   fN = ntracks;
   fLast = 0;
+  fChi2 = 0;
+  fNdf  = 0;
   fIsOwner = kTRUE;
 
   if (ntracks > 0) {
@@ -143,6 +157,17 @@ Bool_t AliTrackResiduals::AddTrackPointArrays(AliTrackPointArray *volarray, AliT
 }
 
 //_____________________________________________________________________________
+void AliTrackResiduals::SetAlignObj(AliAlignObj *alignobj)
+{
+  // Copy the alignment object to be updated in fAlignObj
+  // and flush its parameters
+  if (fAlignObj) delete fAlignObj;
+  fAlignObj = (AliAlignObj *)alignobj->Clone();
+  fAlignObj->SetPars(0,0,0,0,0,0);
+}
+
+
+//_____________________________________________________________________________
 Bool_t AliTrackResiduals::GetTrackPointArrays(Int_t i, AliTrackPointArray* &volarray, AliTrackPointArray* &trackarray) const
 {
   // Provide an access to a pair of track point arrays
index 34d230c..45340a7 100644 (file)
@@ -13,7 +13,8 @@
 
 #include "TObject.h"
 
-class AliAlignObj;
+#include "AliAlignObjAngles.h"
+
 class AliTrackPointArray;
 
 class AliTrackResiduals : public TObject {
@@ -21,14 +22,14 @@ class AliTrackResiduals : public TObject {
  public:
 
   AliTrackResiduals();
-  AliTrackResiduals(Int_t ntracks, AliAlignObj *alignobj);
+  AliTrackResiduals(Int_t ntracks);
   AliTrackResiduals(const AliTrackResiduals &res);
   AliTrackResiduals& operator= (const AliTrackResiduals& res);
   virtual ~AliTrackResiduals();
 
   void   SetNTracks(Int_t ntracks);
-  void   SetAlignObj(AliAlignObj *alignobj) { fAlignObj = alignobj; }
   Bool_t AddTrackPointArrays(AliTrackPointArray *volarray, AliTrackPointArray *trackarray);
+  void   SetAlignObj(AliAlignObj *alignobj);
 
   virtual Bool_t Minimize() = 0;
 
@@ -36,6 +37,8 @@ class AliTrackResiduals : public TObject {
   Int_t  GetNFilledTracks() const { return fLast; }
   Bool_t GetTrackPointArrays(Int_t i, AliTrackPointArray* &volarray, AliTrackPointArray* &trackarray) const;
   AliAlignObj *GetAlignObj() const { return fAlignObj; }
+  Float_t GetChi2() const { return fChi2; }
+  Int_t   GetNdf() const  { return fNdf; }
 
  protected:
 
@@ -46,6 +49,8 @@ class AliTrackResiduals : public TObject {
   AliAlignObj        *fAlignObj;    // Pointer to the volume alignment object to be fitted
   AliTrackPointArray **fVolArray;   //! Pointers to the arrays containing space points
   AliTrackPointArray **fTrackArray; //! Pointers to the arrays containing track extrapolation points
+  Float_t            fChi2;         // Chi2 (or distance) of residuals minimization
+  Int_t              fNdf;          // Number of degrees of freedom
   Bool_t             fIsOwner;      // Track point arrays owned by the object
 
   ClassDef(AliTrackResiduals,1)
index 38bb693..575a818 100644 (file)
@@ -26,6 +26,8 @@
 #include "AliTrackPointArray.h"
 #include "AliTrackResidualsChi2.h"
 
+void Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
+
 ClassImp(AliTrackResidualsChi2)
 
 //______________________________________________________________________________
@@ -43,26 +45,27 @@ Bool_t AliTrackResidualsChi2::Minimize()
 {
   // Implementation of Chi2 based minimization
   // of track residuala sum
+  Double_t arglist[10];
+  Int_t ierflg = 0;
   TMinuit *gMinuit = new TMinuit(6);  //initialize TMinuit
+  arglist[0] = -1;
+  gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg);
+
   gMinuit->SetObjectFit(this);
   gMinuit->SetFCN(Fcn);
 
-  Double_t arglist[10];
-  Int_t ierflg = 0;
-
   arglist[0] = 1;
   gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
 
   // Set starting values and step sizes for parameters
-  Double_t tr[3],rot[3];
-  fAlignObj->GetPars(tr,rot);
-  static Double_t step[6] = {0,0,0,0,0,0};
-  gMinuit->mnparm(0, "x", tr[0], step[0], 0,0,ierflg);
-  gMinuit->mnparm(1, "y", tr[1], step[1], 0,0,ierflg);
-  gMinuit->mnparm(2, "z", tr[2], step[2], 0,0,ierflg);
-  gMinuit->mnparm(3, "psi", rot[0], step[3], 0,0,ierflg);
-  gMinuit->mnparm(4, "theta", rot[1], step[4], 0,0,ierflg);
-  gMinuit->mnparm(5, "phi", rot[2], step[5], 0,0,ierflg);
+  Double_t pars[6] = {0,0,0,0,0,0};
+  Double_t step[6] = {0.0001,0.0001,0.0001,0.0001,0.0001,0.0001};
+  gMinuit->mnparm(0, "dx", pars[0], step[0], 0,0,ierflg);
+  gMinuit->mnparm(1, "dy", pars[1], step[1], 0,0,ierflg);
+  gMinuit->mnparm(2, "dz", pars[2], step[2], 0,0,ierflg);
+  gMinuit->mnparm(3, "psi", pars[3], step[3], 0,0,ierflg);
+  gMinuit->mnparm(4, "theta", pars[4], step[4], 0,0,ierflg);
+  gMinuit->mnparm(5, "phi", pars[5], step[5], 0,0,ierflg);
 
   // Now ready for minimization step
   arglist[0] = 500;
@@ -73,7 +76,7 @@ Bool_t AliTrackResidualsChi2::Minimize()
   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
-  gMinuit->mnprin(3,amin);
+  fChi2 = amin; fNdf -= nvpar;
 
   return kTRUE;
 }
@@ -86,31 +89,21 @@ void AliTrackResidualsChi2::Chi2(Int_t & /* npar */, Double_t * /* gin */, Doubl
   Double_t chi2 = 0;
 
   fAlignObj->SetPars(par[0],par[1],par[2],par[3],par[4],par[5]);
-  TGeoHMatrix m;
-  fAlignObj->GetMatrix(m);
-  Double_t *rot = m.GetRotationMatrix();
-  Double_t *tr  = m.GetTranslation();
 
   AliTrackPoint p1,p2;
-  Double_t dR[3];
-  Float_t xyzvol[3],xyztr[3];
-  Float_t covvol[6],covtr[6];
+
+  Bool_t count = kFALSE;
+  if (fNdf == 0) count = kTRUE;
 
   for (Int_t itrack = 0; itrack < fLast; itrack++) {
     if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
     for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
       fVolArray[itrack]->GetPoint(p1,ipoint);
+      fAlignObj->Transform(p1);
       fTrackArray[itrack]->GetPoint(p2,ipoint);
-      p1.GetXYZ(xyzvol,covvol);
-      p2.GetXYZ(xyztr,covtr);
-
-      for (Int_t i = 0; i < 3; i++)
-       dR[i] = tr[i] 
-               + xyzvol[0]*rot[3*i]
-               + xyzvol[1]*rot[3*i+1]
-               + xyzvol[2]*rot[3*i+2]
-              - xyztr[i];
-      chi2 += dR[0]*dR[0]+dR[1]*dR[1]+dR[2]*dR[2];
+      Float_t residual = p2.GetResidual(p1,kFALSE);
+      chi2 += residual;
+      if (count) fNdf += 3;
     }
   }
   f = chi2;
index 7bb8aa0..17486cc 100644 (file)
@@ -3,8 +3,8 @@
 
 //************************************************************************
 // AliTrackResidualsChi2: derived class (from AliTrackResiduals) which   *
-// implements a minimization of the track residuals based on chi2        *
-// approach.                                                             *
+// implements a MINUIT minimization of the track residuals chi2.         *
+//                                                                       *
 //                                                                       *
 //************************************************************************
 
@@ -15,7 +15,7 @@ class AliTrackResidualsChi2 : public AliTrackResiduals {
 
  public:
   AliTrackResidualsChi2():AliTrackResiduals() { }
-  AliTrackResidualsChi2(Int_t ntracks, AliAlignObj *alignobj):AliTrackResiduals(ntracks,alignobj) { }
+  AliTrackResidualsChi2(Int_t ntracks):AliTrackResiduals(ntracks) { }
   AliTrackResidualsChi2(const AliTrackResidualsChi2 &res):AliTrackResiduals(res) { }
   AliTrackResidualsChi2& operator= (const AliTrackResidualsChi2& res) { ((AliTrackResiduals *)this)->operator=(res); return *this; }
   virtual ~AliTrackResidualsChi2() { }
diff --git a/STEER/AliTrackResidualsFast.cxx b/STEER/AliTrackResidualsFast.cxx
new file mode 100644 (file)
index 0000000..ca7d469
--- /dev/null
@@ -0,0 +1,175 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//   Implementation of the derived class for track residuals
+//   based on linear chi2 minimization (in approximation of
+//   small alignment angles and translations)
+//
+//-----------------------------------------------------------------
+
+#include <TMinuit.h>
+#include <TGeoMatrix.h>
+
+#include "AliAlignObj.h"
+#include "AliTrackPointArray.h"
+#include "AliTrackResidualsFast.h"
+
+ClassImp(AliTrackResidualsFast)
+
+//______________________________________________________________________________
+AliTrackResidualsFast::AliTrackResidualsFast():AliTrackResiduals()
+{
+  // Default constructor
+  for (Int_t i = 0; i < 16; i++) fSum[i] = 0;
+  fSumR = 0;
+}
+
+//______________________________________________________________________________
+AliTrackResidualsFast::AliTrackResidualsFast(Int_t ntracks):
+  AliTrackResiduals(ntracks)
+{
+  // Constructor
+  for (Int_t i = 0; i < 16; i++) fSum[i] = 0;
+  fSumR = 0;
+}
+//______________________________________________________________________________
+AliTrackResidualsFast::AliTrackResidualsFast(const AliTrackResidualsFast &res):
+  AliTrackResiduals(res)
+{
+  // Copy constructor
+  for (Int_t i = 0; i < 16; i++) fSum[i] = res.fSum[i];
+  fSumR = res.fSumR;
+}
+
+//______________________________________________________________________________
+AliTrackResidualsFast &AliTrackResidualsFast::operator= (const AliTrackResidualsFast& res)
+{
+  // Assignment operator
+ ((AliTrackResiduals *)this)->operator=(res);
+ for (Int_t i = 0; i < 16; i++) fSum[i] = res.fSum[i];
+ fSumR = res.fSumR;
+
+ return *this;
+}
+
+//______________________________________________________________________________
+Bool_t AliTrackResidualsFast::Minimize()
+{
+  // Implementation of fast linear Chi2
+  // based minimization of track residuals sum
+  for (Int_t i = 0; i < 16; i++) fSum[i] = 0;
+  fSumR = 0;
+
+  AliTrackPoint p1,p2;
+
+  for (Int_t itrack = 0; itrack < fLast; itrack++) {
+    if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
+    for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
+      fVolArray[itrack]->GetPoint(p1,ipoint);
+      fTrackArray[itrack]->GetPoint(p2,ipoint);
+      AddPoints(p1,p2);
+    }
+  }
+
+  return Update();
+
+  // debug info
+//   Float_t chi2 = 0;
+//   for (Int_t itrack = 0; itrack < fLast; itrack++) {
+//     if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
+//     for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
+//       fVolArray[itrack]->GetPoint(p1,ipoint);
+//       fAlignObj->Transform(p1);
+//       fTrackArray[itrack]->GetPoint(p2,ipoint);
+//       Float_t residual = p2.GetResidual(p1,kFALSE);
+//       chi2 += residual;
+//     }
+//   }
+//   printf("Final chi2 = %f\n",chi2);
+}
+
+//______________________________________________________________________________
+void AliTrackResidualsFast::AddPoints(AliTrackPoint &p, AliTrackPoint &pprime)
+{
+  // Update the sums used for
+  // the linear chi2 minimization
+  Float_t xyz[3],xyzp[3];
+  p.GetXYZ(xyz); pprime.GetXYZ(xyzp);
+  Double_t weight = 1;
+  weight *= weight;
+  fSum[0] += weight;
+  fSum[1] += xyz[0]*weight;
+  fSum[2] += xyz[1]*weight;
+  fSum[3] += xyz[2]*weight;
+  fSum[4] += (xyz[0]*xyz[0]+xyz[1]*xyz[1])*weight;
+  fSum[5] += (xyz[0]*xyz[0]+xyz[2]*xyz[2])*weight;
+  fSum[6] += (xyz[1]*xyz[1]+xyz[2]*xyz[2])*weight;
+  fSum[7] -= xyz[0]*xyz[1]*weight;
+  fSum[8] -= xyz[0]*xyz[2]*weight;
+  fSum[9] -= xyz[1]*xyz[2]*weight;
+  fSum[10] += (xyzp[0]-xyz[0])*weight;
+  fSum[11] += (xyzp[1]-xyz[1])*weight;
+  fSum[12] += (xyzp[2]-xyz[2])*weight;
+  fSum[13] += (xyzp[2]*xyz[1]-xyzp[1]*xyz[2])*weight;
+  fSum[14] += (xyzp[0]*xyz[2]-xyzp[2]*xyz[0])*weight;
+  fSum[15] += (xyzp[1]*xyz[0]-xyzp[0]*xyz[1])*weight;
+
+  fSumR += ((xyzp[0]-xyz[0])*(xyzp[0]-xyz[0])+
+           (xyzp[1]-xyz[1])*(xyzp[1]-xyz[1])+
+           (xyzp[2]-xyz[2])*(xyzp[2]-xyz[2]))*weight;
+
+  fNdf += 3;
+}
+
+//______________________________________________________________________________
+Bool_t AliTrackResidualsFast::Update()
+{
+  // Find the alignment parameters
+  // by using the already accumulated
+  // sums
+  TMatrixDSym     smatrix(6);
+  TMatrixD        sums(1,6);
+
+  smatrix(0,0) = fSum[0]; smatrix(1,1) = fSum[0]; smatrix(2,2)=fSum[0];
+  smatrix(3,3) = fSum[6]; smatrix(4,4) = fSum[5]; smatrix(5,5)=fSum[4];
+
+  smatrix(0,1) = 0; smatrix(0,2) = 0;
+  smatrix(0,3) = 0; smatrix(0,4) = fSum[3]; smatrix(0,5) = -fSum[2];
+  smatrix(1,2) = 0;
+  smatrix(1,3) = -fSum[3]; smatrix(1,4) = 0; smatrix(1,5) = fSum[1];
+  smatrix(2,3) = fSum[2]; smatrix(2,4) = -fSum[1]; smatrix(2,5) = 0;
+
+  smatrix(3,4) = fSum[7]; smatrix(3,5) = fSum[8];
+  smatrix(4,5) = fSum[9];
+
+  sums(0,0)    = fSum[10]; sums(0,1) = fSum[11]; sums(0,2) = fSum[12];
+  sums(0,3)    = fSum[13]; sums(0,4) = fSum[14]; sums(0,5) = fSum[15];
+
+  smatrix.Invert();
+  if (!smatrix.IsValid()) return kFALSE;
+
+  TMatrixD res = sums*smatrix;
+  fAlignObj->SetPars(res(0,0),res(0,1),res(0,2),
+                    TMath::RadToDeg()*res(0,3),
+                    TMath::RadToDeg()*res(0,4),
+                    TMath::RadToDeg()*res(0,5));
+  TMatrixD  tmp = res*sums.T();
+  fChi2 = fSumR - tmp(0,0);
+  fNdf -= 6;
+
+  return kTRUE;
+}
diff --git a/STEER/AliTrackResidualsFast.h b/STEER/AliTrackResidualsFast.h
new file mode 100644 (file)
index 0000000..5d043b9
--- /dev/null
@@ -0,0 +1,42 @@
+#ifndef ALITRACKRESIDUALSFAST_H
+#define ALITRACKRESIDUALSFAST_H
+
+//************************************************************************
+// AliTrackResidualsFast: derived class (from AliTrackResiduals) which   *
+// implements a simple linear minimization of the track residuals chi2.  *
+// The minimization relies on the fact that the alignment parameters     *
+// (angles and translations) are small.                                  *
+//                                                                       *
+//                                                                       *
+//************************************************************************
+
+#include "TMatrixDSym.h"
+#include "TMatrixD.h"
+
+#include "AliAlignObj.h"
+#include "AliTrackResiduals.h"
+
+class AliTrackResidualsFast : public AliTrackResiduals {
+
+ public:
+  AliTrackResidualsFast();
+  AliTrackResidualsFast(Int_t ntracks);
+  AliTrackResidualsFast(const AliTrackResidualsFast &res);
+  AliTrackResidualsFast& operator= (const AliTrackResidualsFast& res);
+  virtual ~AliTrackResidualsFast() { }
+
+  Bool_t Minimize();
+
+ protected:
+
+  void   AddPoints(AliTrackPoint &p, AliTrackPoint &pprime);
+  Bool_t Update();
+
+  Double_t fSum[16]; // Sums used during the chi2 minimization
+  Double_t fSumR;    // Sum of r squared
+
+  ClassDef(AliTrackResidualsFast,1)
+
+};
+
+#endif
index 5b07ce0..3285814 100644 (file)
 #pragma link C++ class AliTrackFitterRieman+;
 #pragma link C++ class AliTrackResiduals+;
 #pragma link C++ class AliTrackResidualsChi2+;
+#pragma link C++ class AliTrackResidualsFast+;
 #pragma link C++ class AliAlignmentTracks+;
 
 #pragma link C++ class  TTreeDataElement+;
index f3847cf..ba20b64 100644 (file)
@@ -30,7 +30,7 @@ AliTriggerCondition.cxx \
 AliTriggerDescriptor.cxx \
 AliCentralTrigger.cxx AliRieman.cxx\
 AliTrackFitter.cxx AliTrackFitterRieman.cxx\
-AliTrackResiduals.cxx AliTrackResidualsChi2.cxx\
+AliTrackResiduals.cxx AliTrackResidualsChi2.cxx AliTrackResidualsFast.cxx\
 AliAlignmentTracks.cxx \
 AliExpression.cxx