Alignment framework (C.Cheshkov). More information is available in http://agenda...
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Nov 2005 15:51:35 +0000 (15:51 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 28 Nov 2005 15:51:35 +0000 (15:51 +0000)
21 files changed:
STEER/AliAlignObj.cxx
STEER/AliAlignObj.h
STEER/AliAlignmentTracks.cxx [new file with mode: 0644]
STEER/AliAlignmentTracks.h [new file with mode: 0644]
STEER/AliESDtrack.cxx
STEER/AliESDtrack.h
STEER/AliReconstruction.cxx
STEER/AliReconstruction.h
STEER/AliTrackFitter.cxx [new file with mode: 0644]
STEER/AliTrackFitter.h [new file with mode: 0644]
STEER/AliTrackFitterRieman.cxx [new file with mode: 0644]
STEER/AliTrackFitterRieman.h [new file with mode: 0644]
STEER/AliTrackPointArray.cxx [new file with mode: 0644]
STEER/AliTrackPointArray.h [new file with mode: 0644]
STEER/AliTrackResiduals.cxx [new file with mode: 0644]
STEER/AliTrackResiduals.h [new file with mode: 0644]
STEER/AliTrackResidualsChi2.cxx [new file with mode: 0644]
STEER/AliTrackResidualsChi2.h [new file with mode: 0644]
STEER/AliTracker.h
STEER/STEERLinkDef.h
STEER/libSTEER.pkg

index 737d467..9212ac6 100644 (file)
 
 #include "AliAlignObj.h"
 //#include "AliLog.h"
-
 ClassImp(AliAlignObj)
 
+Int_t AliAlignObj::fgLayerSize[kLastLayer - kFirstLayer] = {
+  80, 160,  // ITS SPD
+  84, 176,  // ITS SDD
+  748, 950, // ITS SSD
+  36, 36,   // TPC
+  90, 90, 90, 90, 90, 90,  // TRD
+  1,        // TOF ??
+  1, 1,     // PHOS ??
+  1,        // RICH ??
+  1         // MUON ??
+};
+
+const char* AliAlignObj::fgLayerName[kLastLayer - kFirstLayer] = {
+  "ITS inner pixels layer", "ITS outer pixels layer",
+  "ITS inner drifts layer", "ITS outer drifts layer",
+  "ITS inner strips layer", "ITS outer strips layer",
+  "TPC inner chambers layer", "TPC outer chambers layer",
+  "TRD chambers layer 1", "TRD chambers layer 2", "TRD chambers layer 3",
+  "TRD chambers layer 4", "TRD chambers layer 5", "TRD chambers layer 6",
+  "TOF layer",
+  "?","?",
+  "?",
+  "?"
+};
+
 //_____________________________________________________________________________
 AliAlignObj::AliAlignObj():
   fVolUID(0)
index 693c1b7..f4c3d6b 100644 (file)
@@ -19,7 +19,8 @@ class AliAlignObj : public TObject {
   AliAlignObj(const AliAlignObj& theAlignObj);
   AliAlignObj& operator= (const AliAlignObj& theAlignObj);
   virtual ~AliAlignObj();
-  enum ELayerID{kSPD1=1, kSPD2=2,
+  enum ELayerID{kFirstLayer=1,
+               kSPD1=1, kSPD2=2,
                kSDD1=3, kSDD2=4,
                kSSD1=5, kSSD2=6,
                kTPC1=7, kTPC2=8,
@@ -27,7 +28,8 @@ class AliAlignObj : public TObject {
                kTOF=15,
                kPHOS1=16, kPHOS2=17,
                kRICH=18,
-               kMUON=19};
+               kMUON=19,
+               kLastLayer=20};
 
   //Setters
   virtual void SetTranslation(Double_t x, Double_t y, Double_t z) = 0;
@@ -52,6 +54,9 @@ class AliAlignObj : public TObject {
 
   void  Print(Option_t *) const;
 
+  static Int_t       LayerSize(Int_t layer) { return fgLayerSize[layer]; }
+  static const char* LayerName(Int_t layer) { return fgLayerName[layer]; }
+
   static UShort_t LayerToVolUID(ELayerID layerId, Int_t modId);
   static ELayerID VolUIDToLayer(UShort_t voluid, Int_t &modId);
   static ELayerID VolUIDToLayer(UShort_t voluid);
@@ -61,8 +66,11 @@ class AliAlignObj : public TObject {
   Bool_t MatrixToAngles(const Double_t *rot, Double_t *angles) const;
 
   //Volume identifiers
-  TString fVolPath; // Volume path inside TGeo geometry
-  UShort_t fVolUID; // Unique volume ID
+  TString  fVolPath; // Volume path inside TGeo geometry
+  UShort_t fVolUID;  // Unique volume ID
+
+  static Int_t       fgLayerSize[kLastLayer - kFirstLayer];
+  static const char* fgLayerName[kLastLayer - kFirstLayer];
   
   ClassDef(AliAlignObj, 1)
 };
diff --git a/STEER/AliAlignmentTracks.cxx b/STEER/AliAlignmentTracks.cxx
new file mode 100644 (file)
index 0000000..1c8f2b6
--- /dev/null
@@ -0,0 +1,576 @@
+/**************************************************************************
+ * 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 alignment steering class
+//   It provides an access to the track space points
+//   written along the esd tracks. The class enables
+//   the user to plug any track fitter (deriving from
+//   AliTrackFitter class) and minimization fo the
+//   track residual sums (deriving from the AliTrackResiduals).
+//-----------------------------------------------------------------
+
+#include <TChain.h>
+
+#include "AliAlignmentTracks.h"
+#include "AliTrackPointArray.h"
+#include "AliAlignObjAngles.h"
+#include "AliTrackFitterRieman.h"
+#include "AliTrackResidualsChi2.h"
+#include "AliESD.h"
+#include "AliLog.h"
+
+ClassImp(AliAlignmentTracks)
+
+//______________________________________________________________________________
+AliAlignmentTracks::AliAlignmentTracks():
+  fESDChain(0),
+  fPointsFilename("AliTrackPoints.root"),
+  fPointsFile(0),
+  fPointsTree(0),
+  fLastIndex(0),
+  fArrayIndex(0),
+  fIsIndexBuilt(kFALSE),
+  fTrackFitter(0),
+  fMinimizer(0)
+{
+  // Default constructor
+  InitIndex();
+  InitAlignObjs();
+}
+
+//______________________________________________________________________________
+AliAlignmentTracks::AliAlignmentTracks(TChain *esdchain):
+  fESDChain(esdchain),
+  fPointsFilename("AliTrackPoints.root"),
+  fPointsFile(0),
+  fPointsTree(0),
+  fLastIndex(0),
+  fArrayIndex(0),
+  fIsIndexBuilt(kFALSE),
+  fTrackFitter(0),
+  fMinimizer(0)
+{
+  // Constructor in the case
+  // the user provides an already
+  // built TChain with ESD trees
+  InitIndex();
+  InitAlignObjs();
+}
+
+
+//______________________________________________________________________________
+AliAlignmentTracks::AliAlignmentTracks(const char *esdfilename, const char *esdtreename):
+  fPointsFilename("AliTrackPoints.root"),
+  fPointsFile(0),
+  fPointsTree(0),
+  fLastIndex(0),
+  fArrayIndex(0),
+  fIsIndexBuilt(kFALSE),
+  fTrackFitter(0),
+  fMinimizer(0)
+{
+  // Constructor in the case
+  // the user provides a single ESD file
+  // or a directory containing ESD files
+  fESDChain = new TChain(esdtreename);
+  fESDChain->Add(esdfilename);
+
+  InitIndex();
+  InitAlignObjs();
+}
+
+//______________________________________________________________________________
+AliAlignmentTracks::AliAlignmentTracks(const AliAlignmentTracks &alignment):
+  TObject(alignment)
+{
+  // Copy constructor
+  // not implemented
+  AliWarning("Copy constructor not implemented!");
+}
+
+//______________________________________________________________________________
+AliAlignmentTracks& AliAlignmentTracks::operator= (const AliAlignmentTracks& alignment)
+{
+  // Asignment operator
+  // not implemented
+  if(this==&alignment) return *this;
+
+  AliWarning("Asignment operator not implemented!");
+
+  ((TObject *)this)->operator=(alignment);
+
+  return *this;
+}
+
+//______________________________________________________________________________
+AliAlignmentTracks::~AliAlignmentTracks()
+{
+  // Destructor
+  if (fESDChain) delete fESDChain;
+
+  DeleteIndex();
+  DeleteAlignObjs();
+
+  delete fTrackFitter;
+  delete fMinimizer;
+
+  if (fPointsFile) fPointsFile->Close();
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::AddESD(TChain *esdchain)
+{
+  // Add a chain with ESD files
+  if (fESDChain)
+    fESDChain->Add(esdchain);
+  else
+    fESDChain = esdchain;
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::AddESD(const char *esdfilename, const char *esdtreename)
+{
+  // Add a single file or
+  // a directory to the chain
+  // with the ESD files
+  if (fESDChain)
+    fESDChain->AddFile(esdfilename,TChain::kBigNumber,esdtreename);
+  else {
+    fESDChain = new TChain(esdtreename);
+    fESDChain->Add(esdfilename);
+  }
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::ProcessESD()
+{
+  // Analyzes and filters ESD tracks
+  // Stores the selected track space points
+  // into the output file
+
+  if (!fESDChain) return;
+
+  AliESD *esd = 0;
+  fESDChain->SetBranchAddress("ESD",&esd);
+
+  // Open the output file
+  if (fPointsFilename.Data() == "") {
+    AliWarning("Incorrect output filename!");
+    return;
+  }
+
+  TFile *pointsFile = TFile::Open(fPointsFilename,"RECREATE");
+  if (!pointsFile || !pointsFile->IsOpen()) {
+    AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
+    return;
+  }
+
+  TTree *pointsTree = new TTree("spTree", "Tree with track space point arrays");
+  AliTrackPointArray *array = 0;
+  pointsTree->Branch("SP","AliTrackPointArray", &array);
+
+  Int_t ievent = 0;
+  while (fESDChain->GetEntry(ievent++)) {
+    if (!esd) break;
+    Int_t ntracks = esd->GetNumberOfTracks();
+    for (Int_t itrack=0; itrack < ntracks; itrack++) {
+      AliESDtrack * track = esd->GetTrack(itrack);
+      if (!track) continue;
+      UInt_t status = AliESDtrack::kITSpid; 
+      status|=AliESDtrack::kTPCpid; 
+      status|=AliESDtrack::kTRDpid; 
+      if ((track->GetStatus() & status) != status) continue;
+
+      if (track->GetP() < 0.5) continue;
+
+      array = track->GetTrackPointArray();
+      pointsTree->Fill();
+    }
+  }
+
+  if (!pointsTree->Write()) {
+    AliWarning("Can't write the tree with track point arrays!");
+    return;
+  }
+
+  pointsFile->Close();
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::ProcessESD(TSelector *selector)
+{
+  AliWarning(Form("ESD processing based on selector is not yet implemented (%p) !",selector));
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::BuildIndex()
+{
+  // Build index of points tree entries
+  // Used for access based on the volume IDs
+  if (fIsIndexBuilt)
+    ResetIndex();
+  else
+    fIsIndexBuilt = kTRUE;
+
+  TFile *fPointsFile = TFile::Open(fPointsFilename);
+  if (!fPointsFile || !fPointsFile->IsOpen()) {
+    AliWarning(Form("Can't open %s !",fPointsFilename.Data()));
+    return;
+  }
+  
+  //  AliTrackPointArray* array = new AliTrackPointArray;
+  AliTrackPointArray* array = 0;
+  TTree* pointsTree = (TTree*) fPointsFile->Get("spTree");
+  if (!pointsTree) {
+    AliWarning("No pointsTree found!");
+    return;
+  }
+  pointsTree->SetBranchAddress("SP", &array);
+
+  Int_t nArrays = pointsTree->GetEntries();
+  for (Int_t iArray = 0; iArray < nArrays; iArray++)
+    {
+      pointsTree->GetEvent(iArray);
+      if (!array) continue;
+      for (Int_t ipoint = 0; ipoint < array->GetNPoints(); ipoint++) {
+       UShort_t volId = array->GetVolumeID()[ipoint];
+       Int_t modId;
+       Int_t layerId = AliAlignObj::VolUIDToLayer(volId,modId)
+                     - AliAlignObj::kFirstLayer;
+       if (!fArrayIndex[layerId][modId]) {
+         //first entry for this volume
+         fArrayIndex[layerId][modId] = new TArrayI(1000);
+       }
+       else {
+         Int_t size = fArrayIndex[layerId][modId]->GetSize();
+         // If needed allocate new size
+         if (fLastIndex[layerId][modId] >= size)
+           fArrayIndex[layerId][modId]->Set(size + 1000);
+       }
+
+       // Check if the index is already filled
+       Bool_t fillIndex = kTRUE;
+       if (fLastIndex[layerId][modId] != 0) {
+         if ((*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]-1] == iArray)
+           fillIndex = kFALSE;
+       }
+       // Fill the index array and store last filled index
+       if (fillIndex) {
+         (*fArrayIndex[layerId][modId])[fLastIndex[layerId][modId]] = iArray;
+         fLastIndex[layerId][modId]++;
+       }
+      }
+    }
+
+}
+
+// //______________________________________________________________________________
+// void AliAlignmentTracks::BuildIndexLayer(AliAlignObj::ELayerID layer)
+// {
+// }
+
+// //______________________________________________________________________________
+// void AliAlignmentTracks::BuildIndexVolume(UShort_t volid)
+// {
+// }
+
+//______________________________________________________________________________
+void AliAlignmentTracks::InitIndex()
+{
+  // Initialize the index arrays
+  Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
+  fLastIndex = new Int_t*[nLayers];
+  fArrayIndex = new TArrayI**[nLayers];
+  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
+    fLastIndex[iLayer] = new Int_t[AliAlignObj::LayerSize(iLayer)];
+    fArrayIndex[iLayer] = new TArrayI*[AliAlignObj::LayerSize(iLayer)];
+    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
+      fLastIndex[iLayer][iModule] = 0;
+      fArrayIndex[iLayer][iModule] = 0;
+    }
+  }
+}
+
+//______________________________________________________________________________
+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++) {
+    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
+      fLastIndex[iLayer][iModule] = 0;
+    }
+  }
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::DeleteIndex()
+{
+  // Delete the index arrays
+  // Called by the destructor
+  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
+    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++) {
+      if (fArrayIndex[iLayer][iModule]) {
+       delete fArrayIndex[iLayer][iModule];
+       fArrayIndex[iLayer][iModule] = 0;
+      }
+    }
+    delete [] fLastIndex[iLayer];
+    delete [] fArrayIndex[iLayer];
+  }
+  delete [] fLastIndex;
+  delete [] fArrayIndex;
+}
+
+//______________________________________________________________________________
+Bool_t AliAlignmentTracks::ReadAlignObjs(const char *alignobjfilename)
+{
+  // Read alignment object from a file
+  // To be replaced by a call to CDB
+  AliWarning(Form("Method not yet implemented (%s) !",alignobjfilename));
+
+  return kFALSE;
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::InitAlignObjs()
+{
+  // Initialize the alignment objects array
+  Int_t nLayers = AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer;
+  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;
+  }
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::ResetAlignObjs()
+{
+  // Reset the alignment objects array
+  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
+    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
+      fAlignObjs[iLayer][iModule]->SetPars(0,0,0,0,0,0);
+  }
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::DeleteAlignObjs()
+{
+  // Delete the alignment objects array
+  for (Int_t iLayer = 0; iLayer < (AliAlignObj::kLastLayer - AliAlignObj::kFirstLayer); iLayer++) {
+    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(iLayer); iModule++)
+      if (fAlignObjs[iLayer][iModule])
+       delete fAlignObjs[iLayer][iModule];
+    delete [] fAlignObjs[iLayer];
+  }
+  delete [] fAlignObjs;
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::Align(Int_t iterations)
+{
+  // This method is just an example
+  // how one can user AlignLayer and
+  // AlignVolume methods to construct
+  // a custom alignment procedure
+  while (iterations > 0) {
+    // First inward pass
+    AlignLayer(AliAlignObj::kTPC1);
+    AlignLayer(AliAlignObj::kSSD2);
+    AlignLayer(AliAlignObj::kSSD1);
+    AlignLayer(AliAlignObj::kSDD2);
+    AlignLayer(AliAlignObj::kSDD1);
+    AlignLayer(AliAlignObj::kSPD2);
+    AlignLayer(AliAlignObj::kSPD1);
+    // Outward pass
+    AlignLayer(AliAlignObj::kSPD2);
+    AlignLayer(AliAlignObj::kSDD1);
+    AlignLayer(AliAlignObj::kSDD2);
+    AlignLayer(AliAlignObj::kSSD1);
+    AlignLayer(AliAlignObj::kSSD2);
+    AlignLayer(AliAlignObj::kTPC1);
+    AlignLayer(AliAlignObj::kTPC2);
+    AlignLayer(AliAlignObj::kTRD1);
+    AlignLayer(AliAlignObj::kTRD2);
+    AlignLayer(AliAlignObj::kTRD3);
+    AlignLayer(AliAlignObj::kTRD4);
+    AlignLayer(AliAlignObj::kTRD5);
+    AlignLayer(AliAlignObj::kTRD6);
+    AlignLayer(AliAlignObj::kTOF);
+    // Again inward
+    AlignLayer(AliAlignObj::kTRD6);
+    AlignLayer(AliAlignObj::kTRD5);
+    AlignLayer(AliAlignObj::kTRD4);
+    AlignLayer(AliAlignObj::kTRD3);
+    AlignLayer(AliAlignObj::kTRD2);
+    AlignLayer(AliAlignObj::kTRD1);
+    AlignLayer(AliAlignObj::kTPC2);
+  }
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::AlignLayer(AliAlignObj::ELayerID layer,
+                                   AliAlignObj::ELayerID layerRangeMin,
+                                   AliAlignObj::ELayerID layerRangeMax,
+                                   Int_t iterations)
+{
+  // Align detector volumes within
+  // a given layer.
+  // Tracks are fitted only within
+  // the range defined by the user.
+  while (iterations > 0) {
+    for (Int_t iModule = 0; iModule < AliAlignObj::LayerSize(layer); iModule++) {
+      UShort_t volId = AliAlignObj::LayerToVolUID(layer,iModule);
+      AlignVolume(volId,layerRangeMin,layerRangeMax);
+    }
+    iterations--;
+  }
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::AlignVolume(UShort_t volid,
+                                    AliAlignObj::ELayerID layerRangeMin,
+                                    AliAlignObj::ELayerID layerRangeMax)
+{
+  // Align a single detector volume.
+  // Tracks are fitted only within
+  // the range defined by the user.
+
+  // First take the alignment object to be updated
+  Int_t iModule;
+  AliAlignObj::ELayerID iLayer = AliAlignObj::VolUIDToLayer(volid,iModule);
+  AliAlignObj *alignObj = fAlignObjs[iLayer][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;
+
+  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);
+  }
+  minimizer->Minimize();
+
+  UnloadPoints(nArrays, points);
+}
+  
+//______________________________________________________________________________
+Int_t AliAlignmentTracks::LoadPoints(UShort_t volid, AliTrackPointArray** &points)
+{
+  // Load track point arrays with at least
+  // one space point in a given detector
+  // 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) {
+    points = 0;
+    return 0;
+  }
+
+  if (!fPointsTree) {
+    AliWarning("Tree with the space point arrays not initialized!");
+    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];
+  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) {
+      AliWarning("Wrong space point array index!");
+      continue;
+    }
+    Int_t nPoints = array->GetNPoints();
+    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);
+      points[iArray]->AddPoint(iPoint,&p);
+    }
+  }
+
+  return nArrays;
+}
+
+//______________________________________________________________________________
+void AliAlignmentTracks::UnloadPoints(Int_t n, AliTrackPointArray **points)
+{
+  // Unload track point arrays for a given
+  // detector volume
+  for (Int_t iArray = 0; iArray < n; iArray++)
+    delete points[iArray];
+  delete [] points;
+}
+
+//______________________________________________________________________________
+AliTrackFitter *AliAlignmentTracks::CreateFitter()
+{
+  // Check if the user has already supplied
+  // a track fitter object.
+  // If not, create a default one.
+  if (!fTrackFitter)
+    fTrackFitter = new AliTrackFitterRieman;
+
+  return fTrackFitter;
+}
+
+//______________________________________________________________________________
+AliTrackResiduals *AliAlignmentTracks::CreateMinimizer()
+{
+  // Check if the user has already supplied
+  // a track residuals minimizer object.
+  // If not, create a default one.
+  if (!fMinimizer)
+    fMinimizer = new AliTrackResidualsChi2;
+
+  return fMinimizer;
+}
diff --git a/STEER/AliAlignmentTracks.h b/STEER/AliAlignmentTracks.h
new file mode 100644 (file)
index 0000000..5cf6a8a
--- /dev/null
@@ -0,0 +1,89 @@
+#ifndef ALIALIGNMENTTRACKS_H
+#define ALIALIGNMENTTRACKS_H
+
+//*************************************************************************
+// AliAlignmentTracks: main steering class which deals with the alignment *
+// procedures based on reconstructed tracks.                              *
+// More comments will come with the development of the interfaces and     *
+// functionalities of the class.                                          *
+//*************************************************************************
+
+#include <TObject.h>
+
+#include "AliAlignObj.h"
+
+class TChain;
+class AliTrackPointArray;
+class AliAlignObj;
+class AliTrackFitter;
+class AliTrackResiduals;
+
+class AliAlignmentTracks : public TObject {
+
+ public:
+
+  AliAlignmentTracks();
+  AliAlignmentTracks(TChain *esdchain);
+  AliAlignmentTracks(const char *esdfilename, const char *esdtreename = "esdTree");
+  AliAlignmentTracks(const AliAlignmentTracks & alignment);
+  AliAlignmentTracks& operator= (const AliAlignmentTracks& alignment);
+  virtual ~AliAlignmentTracks();
+
+  void AddESD(TChain *esdchain);
+  void AddESD(const char *esdfilename, const char *esdtreename = "esdTree");
+
+  void SetPointsFilename(const char *pointsfilename = "AliTrackPoints.root") { fPointsFilename = pointsfilename; }
+
+  void ProcessESD(TSelector *selector);
+  void ProcessESD();
+
+  void BuildIndex();
+/*   void BuildIndexLayer(AliAlignObj::ELayerID layer); */
+/*   void BuildIndexVolume(UShort_t volid); */
+
+  Bool_t ReadAlignObjs(const char *alignobjfilename = "AlignObjs.root");
+
+  void SetTrackFitter(AliTrackFitter *fitter) { fTrackFitter = fitter; }
+  void SetMinimizer(AliTrackResiduals *minimizer) { fMinimizer = minimizer; }
+
+  void Align(Int_t iterations = 100);
+  void AlignLayer(AliAlignObj::ELayerID layer,
+                 AliAlignObj::ELayerID layerRangeMin = AliAlignObj::kFirstLayer,
+                 AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer,
+                 Int_t iterations = 1);
+  void AlignVolume(UShort_t volid,
+                  AliAlignObj::ELayerID layerRangeMin = AliAlignObj::kFirstLayer,
+                  AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer);
+
+ protected:
+
+  void InitIndex();
+  void ResetIndex();
+  void DeleteIndex();
+
+  void InitAlignObjs();
+  void ResetAlignObjs();
+  void DeleteAlignObjs();
+
+  Int_t LoadPoints(UShort_t volid, AliTrackPointArray** &points);
+  void  UnloadPoints(Int_t n, AliTrackPointArray **points);
+
+  AliTrackFitter *CreateFitter();
+  AliTrackResiduals *CreateMinimizer();
+
+  TChain           *fESDChain;       //! Chain with ESDs
+  TString           fPointsFilename; //  Name of the file containing the track point arrays
+  TFile            *fPointsFile;     //  File containing the track point arrays
+  TTree            *fPointsTree;     //  Tree with the track point arrays
+  Int_t           **fLastIndex;      //! Last filled index in volume arrays
+  TArrayI        ***fArrayIndex;     //! Volume arrays which contains the tree index
+  Bool_t            fIsIndexBuilt;   //  Is points tree index built
+  AliAlignObj    ***fAlignObjs;      //  Array with alignment objects
+  AliTrackFitter   *fTrackFitter;    //  Pointer to the track fitter
+  AliTrackResiduals*fMinimizer;      //  Pointer to track residuals minimizer
+
+  ClassDef(AliAlignmentTracks,1)
+
+};
+
+#endif
index 45ed565..0270c60 100644 (file)
@@ -24,6 +24,7 @@
 
 #include "AliESDtrack.h"
 #include "AliKalmanTrack.h"
+#include "AliTrackPointArray.h"
 #include "AliLog.h"
 
 ClassImp(AliESDtrack)
@@ -78,7 +79,8 @@ AliESDtrack::AliESDtrack() :
   fRICHtheta(0),
   fRICHphi(0),
   fRICHdx(0),
-  fRICHdy(0)
+  fRICHdy(0),
+  fPoints(0)
 {
   //
   // The default ESD constructor 
@@ -176,7 +178,8 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   fRICHtheta(track.fRICHtheta),
   fRICHphi(track.fRICHphi),
   fRICHdx(track.fRICHdx),
-  fRICHdy(track.fRICHdy)
+  fRICHdy(track.fRICHdy),
+  fPoints(track.fPoints)
 {
   //
   //copy constructor
@@ -232,7 +235,8 @@ AliESDtrack::~AliESDtrack(){
   //
   //printf("Delete track\n");
   delete fITStrack;
-  delete fTRDtrack;  
+  delete fTRDtrack; 
+  delete fPoints;
 }
 
 //_______________________________________________________________________
@@ -341,6 +345,7 @@ void AliESDtrack::MakeMiniESDtrack(){
   fRICHdx = 0;     
   fRICHdy = 0;      
 
+  fPoints = 0;
 } 
 //_______________________________________________________________________
 Double_t AliESDtrack::GetMass() const {
@@ -776,6 +781,58 @@ void AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const
  
 }
 
+Int_t AliESDtrack::GetNcls(Int_t idet) const
+{
+  // Get number of clusters by subdetector index
+  //
+  Int_t ncls = 0;
+  switch(idet){
+  case 0:
+    ncls = fITSncls;
+    break;
+  case 1:
+    ncls = fTPCncls;
+    break;
+  case 2:
+    ncls = fTRDncls;
+    break;
+  case 3:
+    if (fTOFindex != 0)
+      ncls = 1;
+    break;
+  default:
+    break;
+  }
+  return ncls;
+}
+
+Int_t AliESDtrack::GetClusters(Int_t idet, UInt_t *idx) const
+{
+  // Get cluster index array by subdetector index
+  //
+  Int_t ncls = 0;
+  switch(idet){
+  case 0:
+    ncls = GetITSclusters(idx);
+    break;
+  case 1:
+    ncls = GetTPCclusters((Int_t *)idx);
+    break;
+  case 2:
+    ncls = GetTRDclusters(idx);
+    break;
+  case 3:
+    if (fTOFindex != 0) {
+      idx[0] = GetTOFcluster();
+      ncls = 1;
+    }
+    break;
+  default:
+    break;
+  }
+  return ncls;
+}
+
 void  AliESDtrack::GetTRDExternalParameters(Double_t &x, Double_t&alpha, Double_t p[5], Double_t cov[15]) const
 {
   //
index 6650076..3fba245 100644 (file)
@@ -28,6 +28,7 @@
 #include <TVector3.h>
 
 class AliKalmanTrack;
+class AliTrackPointArray;
 
 const Int_t kNPlane = 6;
 
@@ -87,7 +88,10 @@ public:
   void GetInnerExternalParameters(Double_t &x, Double_t p[5]) const;//skowron
   void GetInnerExternalCovariance(Double_t cov[15]) const;//skowron
   Double_t GetInnerAlpha() const {return fIalpha;}
-  
+  Int_t GetNcls(Int_t idet) const;
+  Int_t GetClusters(Int_t idet, UInt_t *idx) const;
   void SetITSpid(const Double_t *p);
   void SetITSChi2MIP(const Float_t *chi2mip);
   void SetITStrack(AliKalmanTrack * track){fITStrack=track;}
@@ -102,6 +106,7 @@ public:
   void SetTPCpid(const Double_t *p);
   void GetTPCpid(Double_t *p) const;
   void SetTPCPoints(Float_t points[4]){for (Int_t i=0;i<4;i++) fTPCPoints[i]=points[i];}
+  Float_t GetTPCPoints(Int_t i){return fTPCPoints[i];}
   void SetKinkIndexes(Int_t points[3]) {for (Int_t i=0;i<3;i++) fKinkIndexes[i] = points[i];}
   void SetV0Indexes(Int_t points[3]) {for (Int_t i=0;i<3;i++) fV0Indexes[i] = points[i];}
   Float_t GetTPCsignal() const {return fTPCsignal;}
@@ -196,6 +201,9 @@ public:
   Bool_t IsPHOS()  const {return fFlags&kPHOSpid;}
   Bool_t IsEMCAL() const {return fFlags&kEMCALpid;}
 
+  void   SetTrackPointArray(AliTrackPointArray *points) { fPoints = points; }
+  AliTrackPointArray *GetTrackPointArray() const { return fPoints; }
+
   virtual void Print(Option_t * opt) const ; 
 
   enum {
@@ -314,8 +322,10 @@ protected:
   Float_t fRICHphi;        // phi of the track extrapolated to the RICH
   Float_t fRICHdx;         // x of the track impact minus x of the MIP
   Float_t fRICHdy;         // y of the track impact minus y of the MIP
-       
-  ClassDef(AliESDtrack,17)  //ESDtrack 
+
+  AliTrackPointArray *fPoints; // Array which contains the track space points in the global frame
+
+  ClassDef(AliESDtrack,18)  //ESDtrack 
 };
 
 #endif 
index 93d9bad..6333f97 100644 (file)
 #include "AliDetectorTag.h"
 #include "AliEventTag.h"
 
-
+#include "AliTrackPointArray.h"
 
 ClassImp(AliReconstruction)
 
@@ -167,7 +167,9 @@ AliReconstruction::AliReconstruction(const char* gAliceFilename,
   fRunLoader(NULL),
   fRawReader(NULL),
 
-  fVertexer(NULL)
+  fVertexer(NULL),
+
+  fWriteAlignmentData(kFALSE)
 {
 // create reconstruction object with default parameters
   
@@ -200,7 +202,9 @@ AliReconstruction::AliReconstruction(const AliReconstruction& rec) :
   fRunLoader(NULL),
   fRawReader(NULL),
 
-  fVertexer(NULL)
+  fVertexer(NULL),
+
+  fWriteAlignmentData(rec.fWriteAlignmentData)
 {
 // copy constructor
 
@@ -750,6 +754,11 @@ Bool_t AliReconstruction::RunTracking(AliESD*& esd)
     }
   }
 
+  // write space-points to the ESD in case alignment data output
+  // is switched on
+  if (fWriteAlignmentData)
+    WriteAlignmentData(esd);
+
   // pass 3: TRD + TPC + ITS refit inwards
   for (Int_t iDet = 2; iDet >= 0; iDet--) {
     if (!fTracker[iDet]) continue;
@@ -1419,3 +1428,56 @@ void AliReconstruction::CreateTag(TFile* file)
   delete evTag;
 }
 
+void AliReconstruction::WriteAlignmentData(AliESD* esd)
+{
+  // Write space-points which are then used in the alignment procedures
+  // For the moment only ITS, TRD and TPC
+
+  // Load TOF clusters
+  fLoader[3]->LoadRecPoints("read");
+  TTree* tree = fLoader[3]->TreeR();
+  if (!tree) {
+    AliError(Form("Can't get the %s cluster tree", fgkDetectorName[3]));
+    return;
+  }
+  fTracker[3]->LoadClusters(tree);
+
+  Int_t ntracks = esd->GetNumberOfTracks();
+  for (Int_t itrack = 0; itrack < ntracks; itrack++)
+    {
+      AliESDtrack *track = esd->GetTrack(itrack);
+      Int_t nsp = 0;
+      UInt_t idx[200];
+      for (Int_t iDet = 3; iDet >= 0; iDet--)
+       nsp += track->GetNcls(iDet);
+      if (nsp) {
+       AliTrackPointArray *sp = new AliTrackPointArray(nsp);
+       track->SetTrackPointArray(sp);
+       Int_t isptrack = 0;
+       for (Int_t iDet = 3; iDet >= 0; iDet--) {
+         AliTracker *tracker = fTracker[iDet];
+         if (!tracker) continue;
+         Int_t nspdet = track->GetNcls(iDet);
+         cout<<iDet<<" "<<nspdet<<endl;
+         if (nspdet <= 0) continue;
+         track->GetClusters(iDet,idx);
+         AliTrackPoint p;
+         Int_t isp = 0;
+         Int_t isp2 = 0;
+         while (isp < nspdet) {
+           Bool_t isvalid = tracker->GetTrackPoint(idx[isp2],p); isp2++;
+           if (!isvalid) continue;
+           sp->AddPoint(isptrack,&p); isptrack++; isp++;
+         }
+         //      for (Int_t isp = 0; isp < nspdet; isp++) {
+           //      AliCluster *cl = tracker->GetCluster(idx[isp]);
+           //      UShort_t volid = tracker->GetVolumeID(idx[isp]);
+         //        tracker->GetTrackPoint(idx[isp],p);
+         //        sp->AddPoint(isptrack,&p); isptrack++;
+         //      }
+       }       
+      }
+    }
+  fTracker[3]->UnloadClusters();
+  fLoader[3]->UnloadRecPoints();
+}
index ba2f126..157b1b8 100644 (file)
@@ -77,6 +77,8 @@ public:
   Bool_t         Run(Int_t firstEvent, Int_t lastEvent = -1)
     {return Run(NULL, firstEvent, lastEvent);};
 
+  void SetWriteAlignmentData(){fWriteAlignmentData=kTRUE;}
+
 private:
   Bool_t         RunLocalReconstruction(const TString& detectors);
   Bool_t         RunLocalEventReconstruction(const TString& detectors);
@@ -100,6 +102,7 @@ private:
   void           CreateTag(TFile* file);
   //==========================================//
 
+  void           WriteAlignmentData(AliESD* esd);
 
 
 
@@ -127,7 +130,9 @@ private:
   AliVertexer*   fVertexer;                //! vertexer for ITS
   AliTracker*    fTracker[fgkNDetectors];  //! trackers
 
-  ClassDef(AliReconstruction, 5)      // class for running the reconstruction
+  Bool_t         fWriteAlignmentData; // write track space-points flag
+
+  ClassDef(AliReconstruction, 6)      // class for running the reconstruction
 };
 
 #endif
diff --git a/STEER/AliTrackFitter.cxx b/STEER/AliTrackFitter.cxx
new file mode 100644 (file)
index 0000000..2385322
--- /dev/null
@@ -0,0 +1,111 @@
+/**************************************************************************
+ * 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 base class for fast track fitters
+//
+//
+//-----------------------------------------------------------------
+
+#include <TMatrixDSym.h>
+
+#include "AliTrackFitter.h"
+#include "AliTrackPointArray.h"
+
+ClassImp(AliTrackFitter)
+
+//_____________________________________________________________________________
+AliTrackFitter::AliTrackFitter()
+{
+  // default constructor
+  //
+  for (Int_t i=0;i<6;i++) fParams[i] = 0;
+  fCov = 0;
+  fPoints = 0;
+  fIsOwner = kFALSE;
+}
+
+//_____________________________________________________________________________
+AliTrackFitter::AliTrackFitter(AliTrackPointArray *array, Bool_t owner)
+{
+  // constructor from space points array
+  //
+  for (Int_t i=0;i<6;i++) fParams[i] = 0;
+  fCov = new TMatrixDSym(6);
+  fIsOwner = kFALSE;
+  SetTrackPointArray(array,owner);
+}
+
+//_____________________________________________________________________________
+AliTrackFitter::AliTrackFitter(const AliTrackFitter &fitter):
+  TObject(fitter)
+{
+  // Copy constructor
+  //
+  for (Int_t i=0;i<6;i++) fParams[i] = fitter.fParams[i];
+  fCov = new TMatrixDSym(*fitter.fCov);
+  fIsOwner = kFALSE;
+  SetTrackPointArray(fitter.fPoints,fitter.fIsOwner);
+}
+
+//_____________________________________________________________________________
+AliTrackFitter &AliTrackFitter::operator =(const AliTrackFitter& fitter)
+{
+  // assignment operator
+  //
+  if(this==&fitter) return *this;
+
+  for (Int_t i=0;i<6;i++) fParams[i] = fitter.fParams[i];
+  fCov = new TMatrixDSym(*fitter.fCov);
+  fIsOwner = kFALSE;
+  SetTrackPointArray(fitter.fPoints);
+  
+  return *this;
+}
+
+//_____________________________________________________________________________
+AliTrackFitter::~AliTrackFitter()
+{
+  if (fIsOwner)
+    delete fPoints;
+  delete fCov;
+}
+
+//_____________________________________________________________________________
+void AliTrackFitter::Reset()
+{
+  for (Int_t i=0;i<6;i++) fParams[i] = 0;
+  delete fCov;
+  fCov = new TMatrixDSym(6);
+}
+
+void AliTrackFitter::SetTrackPointArray(AliTrackPointArray *array, Bool_t owner)
+{
+  // Load space points from array
+  // By default we don't copy them but
+  // just put the pointers to them
+  Reset();
+
+  if (fIsOwner) delete fPoints;
+
+  if (owner) {
+    fPoints = new AliTrackPointArray(*array);
+    fIsOwner = kTRUE;
+  }
+  else {
+    fPoints = array;
+    fIsOwner = kFALSE;
+  }
+}
diff --git a/STEER/AliTrackFitter.h b/STEER/AliTrackFitter.h
new file mode 100644 (file)
index 0000000..73caed1
--- /dev/null
@@ -0,0 +1,55 @@
+#ifndef ALITRACKFITTER_H
+#define ALITRACKFITTER_H
+
+/*************************************************************************
+ * AliTrackFitter: base class for the fast track fitters                 *
+ *                                                                       *
+ *                                                                       *
+ *                                                                       *
+ *************************************************************************/
+
+#include "TObject.h"
+
+#include "AliTrackPointArray.h"
+#include "AliAlignObj.h"
+
+class TMatrixDSym;
+
+class AliTrackFitter : public TObject {
+
+ public:
+
+  AliTrackFitter();
+  AliTrackFitter(AliTrackPointArray *array, Bool_t owner = kTRUE);
+  AliTrackFitter(const AliTrackFitter &fitter);
+  AliTrackFitter& operator= (const AliTrackFitter& fitter);
+  virtual ~AliTrackFitter();
+
+  virtual void   Reset();
+  virtual void   SetTrackPointArray(AliTrackPointArray *array, Bool_t owner = kTRUE);
+  virtual Bool_t Fit(UShort_t volId,
+                    AliTrackPointArray *pVolId, AliTrackPointArray *pTrack,
+                    AliAlignObj::ELayerID layerRangeMin = AliAlignObj::kFirstLayer,
+                    AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer) = 0;
+  virtual Bool_t GetPCA(const AliTrackPoint &pIn, AliTrackPoint &pOut) const = 0;
+
+  const Float_t* GetX() const {return fPoints->GetX();}
+  const Float_t* GetY() const {return fPoints->GetY();}
+  const Float_t* GetZ() const {return fPoints->GetZ();}
+  const Double_t* GetParam() const {return &fParams[0];}
+  const TMatrixDSym &  GetCovariance() const {return *fCov;}
+
+ protected:
+
+  Double_t      fParams[6];    // Track parameters
+  TMatrixDSym   *fCov;         // Track cov matrix
+  AliTrackPointArray *fPoints; // Pointer to the array with track space points
+  Bool_t  fIsOwner;            // Is the object owner of the space points array
+
+ private:
+
+  ClassDef(AliTrackFitter,1) // Abstract class of fast track fitters
+
+};
+
+#endif
diff --git a/STEER/AliTrackFitterRieman.cxx b/STEER/AliTrackFitterRieman.cxx
new file mode 100644 (file)
index 0000000..793d77b
--- /dev/null
@@ -0,0 +1,402 @@
+#include "TMatrixDSym.h"
+#include "TMatrixD.h"
+#include "AliTrackFitterRieman.h"
+
+ClassImp(AliTrackFitterRieman)
+
+AliTrackFitterRieman::AliTrackFitterRieman():
+  AliTrackFitter()
+{
+  //
+  // default constructor
+  //
+  fAlpha = 0.;
+  for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
+  for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
+  fConv = kFALSE;
+}
+
+
+AliTrackFitterRieman::AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner):
+  AliTrackFitter(array,owner)
+{
+  //
+  // Constructor
+  //
+  fAlpha = 0.;
+  for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
+  for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
+  fConv = kFALSE;
+}
+
+AliTrackFitterRieman::AliTrackFitterRieman(const AliTrackFitterRieman &rieman):
+  AliTrackFitter(rieman)
+{
+  //
+  // copy constructor
+  //
+  fAlpha = rieman.fAlpha;
+  for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
+  for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
+  fConv = rieman.fConv;
+}
+
+//_____________________________________________________________________________
+AliTrackFitterRieman &AliTrackFitterRieman::operator =(const AliTrackFitterRieman& rieman)
+{
+  // assignment operator
+  //
+  if(this==&rieman) return *this;
+  ((AliTrackFitter *)this)->operator=(rieman);
+
+  fAlpha = rieman.fAlpha;
+  for (Int_t i=0;i<9;i++) fSumXY[i]  = rieman.fSumXY[i];
+  for (Int_t i=0;i<9;i++) fSumXZ[i]  = rieman.fSumXZ[i];
+  fConv = rieman.fConv;
+
+  return *this;
+}
+
+AliTrackFitterRieman::~AliTrackFitterRieman()
+{
+  // destructor
+  //
+}
+
+void AliTrackFitterRieman::Reset()
+{
+  // Reset the track parameters and
+  // rieman sums
+  AliTrackFitter::Reset();
+  fAlpha = 0.;
+  for (Int_t i=0;i<9;i++) fSumXY[i] = 0;
+  for (Int_t i=0;i<9;i++) fSumXZ[i] = 0;
+  fConv =kFALSE;
+}
+
+Bool_t AliTrackFitterRieman::Fit(UShort_t volId,
+                                AliTrackPointArray *pVolId, AliTrackPointArray *pTrack,
+                                AliAlignObj::ELayerID layerRangeMin,
+                                AliAlignObj::ELayerID layerRangeMax)
+{
+  // Fit the track points. The method takes as an input
+  // the id (volid) of the volume to be skipped from fitting.
+  // The following two parameters are used to define the
+  // range of volumes to be used in the fitting
+  // As a result two AliTrackPointArray's obects are filled.
+  // The first one contains the space points with
+  // volume id = volid. The second array of points represents
+  // the track extrapolation corresponding to the space points
+  // in the first array. The two arrays can be used to find
+  // the residuals in the volid and consequently construct a
+  // chi2 function to be minimized during the alignment
+  // procedures. For the moment the track extrapolation is taken
+  // as follows: in XY plane - at the CDA between track circle
+  // and the space point; in Z - the track extrapolation on the Z
+  // plane defined by the space point.
+
+  pVolId = pTrack = 0x0;
+  fConv = kFALSE;
+
+  Int_t npoints = fPoints->GetNPoints();
+  if (npoints < 3) return kFALSE;
+
+  AliTrackPoint p;
+  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;
+  Int_t *pindex = new Int_t[npoints];
+  for (Int_t ipoint = 0; ipoint < npoints; ipoint++)
+    {
+      fPoints->GetPoint(p,ipoint);
+      UShort_t iVolId = p.GetVolumeID();
+      if (iVolId == 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 (npused < 3) {
+    delete [] pindex;
+    return kFALSE;
+  }
+
+  Update();
+
+  if (!fConv) {
+    delete [] pindex;
+    return kFALSE;
+  }
+
+  if ((fParams[0] == 0) ||
+      ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0)) {
+    delete [] pindex;
+    return kFALSE;
+  }
+
+  pVolId = new AliTrackPointArray(npVolId);
+  pTrack = 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);
+      }
+    }  
+
+  delete [] pindex;
+
+  return kTRUE;
+}
+
+void AliTrackFitterRieman::AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz)
+{
+  //
+  //  Rieman update
+  //
+  //------------------------------------------------------
+  // XY direction
+  //
+  //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
+  //
+  //  (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
+  //
+  //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
+  //     
+  //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
+  //
+  //  final linear equation :   a + u*b +t*c - v =0;
+  //
+  // Minimization :
+  //
+  // sum( (a + ui*b +ti*c - vi)^2)/(sigmai)^2 = min;
+  //
+  // where sigmai is the error of  maesurement  (a + ui*b +ti*c - vi)
+  //
+  // neglecting error of xi, and supposing  xi>>yi    sigmai ~ sigmaVi ~ 2*sigmay*t  
+  //
+  //
+  // XY part
+  //
+  Double_t  t  =  x*x+y*y;
+  if (t<2) return;
+  t            = 1./t;
+  Double_t  u  =  2.*x*t;
+  Double_t  v  =  2.*y*t;
+  Double_t  error = 2.*sy*t;
+  error *=error;
+  Double_t weight = 1./error;
+  fSumXY[0] +=weight;
+  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;
+  //
+  // 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;
+}
+
+void AliTrackFitterRieman::Update(){
+  //
+  //  Rieman update
+  //
+  //
+  for (Int_t i=0;i<6;i++)fParams[i]=0;
+  Int_t conv=0;
+  //
+  // XY part
+  //
+  TMatrixDSym     smatrix(3);
+  TMatrixD        sums(1,3);
+  //
+  //   smatrix(0,0) = s0; smatrix(1,1)=su2; smatrix(2,2)=st2;
+  //   smatrix(0,1) = su; smatrix(0,2)=st; smatrix(1,2)=sut;
+  //   sums(0,0)    = sv; sums(0,1)=suv; sums(0,2)=svt;
+
+  smatrix(0,0) = fSumXY[0]; smatrix(1,1)=fSumXY[4]; smatrix(2,2)=fSumXY[5];
+  smatrix(0,1) = fSumXY[1]; smatrix(0,2)=fSumXY[3]; smatrix(1,2)=fSumXY[7];
+  sums(0,0)    = fSumXY[2]; sums(0,1)   =fSumXY[6]; sums(0,2)   =fSumXY[8];
+  smatrix.Invert();
+  if (smatrix.IsValid()){
+    for (Int_t i=0;i<3;i++)
+      for (Int_t j=0;j<=i;j++){
+       (*fCov)(i,j)=smatrix(i,j);
+      }
+    TMatrixD  res = sums*smatrix;
+    fParams[0] = res(0,0);
+    fParams[1] = res(0,1);
+    fParams[2] = res(0,2);
+    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);
+    }
+    conv++;
+  }
+
+  //  (x-x0)^2+(y-y0)^2-R^2=0 ===>
+  //
+  //  (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
+  //
+  //  1 - u*x0 - v*y0 - t *D0 =0 ;  - linear equation
+  //     
+  //  next substition   a = 1/y0    b = -x0/y0   c = -D0/y0
+  //  final linear equation :   a + u*b +t*c - v =0;
+  //
+  //  fParam[0]  = 1/y0
+  //  fParam[1]  = -x0/y0
+  //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
+  if (conv>1) fConv =kTRUE;
+  else
+    fConv=kFALSE;
+}
+
+Double_t AliTrackFitterRieman::GetYat(Double_t x){
+  if (!fConv) return 0.;
+  Double_t res2 = (x*fParams[0]+fParams[1]);
+  res2*=res2;
+  res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
+  if (res2<0) return 0;
+  res2 = TMath::Sqrt(res2);
+  res2 = (1-res2)/fParams[0];
+  return res2;
+}
+
+Double_t AliTrackFitterRieman::GetDYat(Double_t x){
+  if (!fConv) return 0.;
+  Double_t x0 = -fParams[1]/fParams[0];
+  if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<0) return 0;
+  Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
+  if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)<=0) return 0;
+  Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
+  if (fParams[0]<0) res*=-1.;
+  return res;
+}
+
+
+
+Double_t AliTrackFitterRieman::GetZat(Double_t x){
+  if (!fConv) return 0.;
+  Double_t x0 = -fParams[1]/fParams[0];
+  if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
+  Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
+  Double_t phi  = TMath::ASin((x-x0)*Rm1);
+  Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
+  Double_t dphi = (phi-phi0);
+  Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
+  return res;
+}
+
+Double_t AliTrackFitterRieman::GetDZat(Double_t x){
+  if (!fConv) return 0.;
+  Double_t x0 = -fParams[1]/fParams[0]; 
+  if  (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
+  Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
+  Double_t res = fParams[4]/TMath::Cos(TMath::ASin((x-x0)*Rm1));
+  return res;
+}
+
+
+Double_t AliTrackFitterRieman::GetC(){
+  return fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1);
+}
+
+Bool_t AliTrackFitterRieman::GetXYZat(Double_t r, Float_t *xyz){
+  if (!fConv) return kFALSE;
+  Double_t res2 = (r*fParams[0]+fParams[1]);
+  res2*=res2;
+  res2 = 1.-fParams[2]*fParams[0]+fParams[1]*fParams[1]-res2;
+  if (res2<0) return kFALSE;
+  res2 = TMath::Sqrt(res2);
+  res2 = (1-res2)/fParams[0];
+
+  if (!fConv) return kFALSE;
+  Double_t x0 = -fParams[1]/fParams[0];
+  if (-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1<=0) return 0;
+  Double_t Rm1  = fParams[0]/TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1); 
+  Double_t phi  = TMath::ASin((r-x0)*Rm1);
+  Double_t phi0 = TMath::ASin((0.-x0)*Rm1);
+  Double_t dphi = (phi-phi0);
+  Double_t res = fParams[3]+fParams[4]*dphi/Rm1;
+
+  Double_t sin = TMath::Sin(fAlpha);
+  Double_t cos = TMath::Cos(fAlpha);
+  xyz[0] = r*cos - res2*sin;
+  xyz[1] = res2*cos + r*sin;
+  xyz[2] = res;
+
+  return kTRUE;
+}
+
+Bool_t AliTrackFitterRieman::GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const
+{
+  // Get the closest to a given spacepoint track trajectory point
+  // Look for details in the description of the Fit() method
+
+  if (!fConv) return kFALSE;
+
+  // First X and Y coordinates
+  Double_t sin = TMath::Sin(fAlpha);
+  Double_t cos = TMath::Cos(fAlpha);
+  //  fParam[0]  = 1/y0
+  //  fParam[1]  = -x0/y0
+  //  fParam[2]  = - (R^2 - x0^2 - y0^2)/y0
+  if (fParams[0] == 0) return kFALSE;
+  Double_t x0 = -fParams[1]/fParams[0]*cos -         1./fParams[0]*sin;
+  Double_t y0 =          1./fParams[0]*cos - fParams[1]/fParams[0]*sin;
+  if ((-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1) <= 0) return kFALSE;
+  Double_t R  = TMath::Sqrt(-fParams[2]*fParams[0]+fParams[1]*fParams[1]+1)/
+                fParams[0];
+
+  Double_t x  = p.GetX();
+  Double_t y  = p.GetY();
+  Double_t dR = TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));
+  Double_t xprime = TMath::Abs(R)*(x-x0)/dR + x0;
+  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 dphi = (phi-phi0);
+  Double_t zprime = fParams[3]+fParams[4]*dphi*R;
+
+  p2.SetXYZ(xprime,yprime,zprime);
+
+  return kTRUE;
+}
diff --git a/STEER/AliTrackFitterRieman.h b/STEER/AliTrackFitterRieman.h
new file mode 100644 (file)
index 0000000..1c9c1fb
--- /dev/null
@@ -0,0 +1,46 @@
+#ifndef ALITRACKFITTERRIEMAN_H
+#define ALITRACKFITTERRIEMAN_H
+
+#include "TMatrixDSym.h"
+
+#include "AliTrackFitter.h"
+
+class AliTrackFitterRieman : public AliTrackFitter{
+ public:
+  AliTrackFitterRieman();
+  AliTrackFitterRieman(AliTrackPointArray *array, Bool_t owner = kTRUE);
+  AliTrackFitterRieman(const AliTrackFitterRieman &rieman);
+  AliTrackFitterRieman &operator =(const AliTrackFitterRieman& rieman);
+  virtual ~AliTrackFitterRieman();
+
+  Bool_t Fit(UShort_t volId,
+            AliTrackPointArray *pVolId, AliTrackPointArray *pTrack,
+            AliAlignObj::ELayerID layerRangeMin = AliAlignObj::kFirstLayer,
+            AliAlignObj::ELayerID layerRangeMax = AliAlignObj::kLastLayer);
+  Bool_t GetPCA(const AliTrackPoint &p, AliTrackPoint &p2) const;
+
+  void Reset();
+  void AddPoint(Float_t x, Float_t y, Float_t z, Float_t sy, Float_t sz);
+  void Update();
+
+  Double_t GetC(); 
+  Double_t GetYat(Double_t x);
+  Double_t GetZat(Double_t x);
+  Double_t GetDYat(Double_t x);
+  Double_t GetDZat(Double_t x);
+  Bool_t   GetXYZat(Double_t r, Float_t *xyz);
+
+ protected:
+
+  Double_t      fAlpha;     //angle to transform to the fitting coordinate system
+  Double_t      fSumXY[9];  //sums for XY part
+  Double_t      fSumXZ[9];  //sums for XZ part
+  Bool_t        fConv;      // indicates convergation
+
+ private:
+
+  ClassDef(AliTrackFitterRieman,1)  // Fast fit of helices on ITS RecPoints
+
+};
+
+#endif
diff --git a/STEER/AliTrackPointArray.cxx b/STEER/AliTrackPointArray.cxx
new file mode 100644 (file)
index 0000000..10b0380
--- /dev/null
@@ -0,0 +1,230 @@
+/**************************************************************************
+ * 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.                  *
+ **************************************************************************/
+
+//////////////////////////////////////////////////////////////////////////////
+//                          Class AliTrackPointArray                        //
+//   This class contains the ESD track space-points which are used during   //
+//   the alignment procedures. Each space-point consist of 3 coordinates    //
+//   (and their errors) and the index of the sub-detector which contains    //
+//   the space-point.                                                       //
+//   cvetan.cheshkov@cern.ch 3/11/2005                                      //
+//////////////////////////////////////////////////////////////////////////////
+
+#include "AliTrackPointArray.h"
+
+ClassImp(AliTrackPointArray)
+
+//______________________________________________________________________________
+AliTrackPointArray::AliTrackPointArray()
+{
+  fNPoints = fSize = 0;
+  fX = fY = fZ = 0;
+  fVolumeID = 0;
+  fCov = 0;
+}
+
+//______________________________________________________________________________
+AliTrackPointArray::AliTrackPointArray(Int_t npoints):
+  fNPoints(npoints)
+{
+  // Constructor
+  //
+  fSize = 6*npoints;
+  fX = new Float_t[npoints];
+  fY = new Float_t[npoints];
+  fZ = new Float_t[npoints];
+  fVolumeID = new UShort_t[npoints];
+  fCov = new Float_t[fSize];
+}
+
+//______________________________________________________________________________
+AliTrackPointArray::AliTrackPointArray(const AliTrackPointArray &array):
+  TObject(array)
+{
+  // Copy constructor
+  //
+  fNPoints = array.fNPoints;
+  fSize = array.fSize;
+  fX = new Float_t[fNPoints];
+  fY = new Float_t[fNPoints];
+  fZ = new Float_t[fNPoints];
+  fVolumeID = new UShort_t[fNPoints];
+  fCov = new Float_t[fSize];
+  memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
+  memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
+  memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
+  memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
+  memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
+}
+
+//_____________________________________________________________________________
+AliTrackPointArray &AliTrackPointArray::operator =(const AliTrackPointArray& array)
+{
+  // assignment operator
+  //
+  if(this==&array) return *this;
+  ((TObject *)this)->operator=(array);
+
+  fNPoints = array.fNPoints;
+  fSize = array.fSize;
+  fX = new Float_t[fNPoints];
+  fY = new Float_t[fNPoints];
+  fZ = new Float_t[fNPoints];
+  fVolumeID = new UShort_t[fNPoints];
+  fCov = new Float_t[fSize];
+  memcpy(fX,array.fX,fNPoints*sizeof(Float_t));
+  memcpy(fY,array.fY,fNPoints*sizeof(Float_t));
+  memcpy(fZ,array.fZ,fNPoints*sizeof(Float_t));
+  memcpy(fVolumeID,array.fVolumeID,fNPoints*sizeof(UShort_t));
+  memcpy(fCov,array.fCov,fSize*sizeof(Float_t));
+
+  return *this;
+}
+
+//______________________________________________________________________________
+AliTrackPointArray::~AliTrackPointArray()
+{
+  // Destructor
+  //
+  delete [] fX;
+  delete [] fY;
+  delete [] fZ;
+  delete [] fVolumeID;
+  delete [] fCov;
+}
+
+
+//______________________________________________________________________________
+Bool_t AliTrackPointArray::AddPoint(Int_t i, const AliTrackPoint *p)
+{
+  // Add a point to the array at position i
+  //
+  if (i >= fNPoints) return kFALSE;
+  fX[i] = p->GetX();
+  fY[i] = p->GetY();
+  fZ[i] = p->GetZ();
+  fVolumeID[i] = p->GetVolumeID();
+  memcpy(&fCov[6*i],p->GetCov(),6*sizeof(Float_t));
+  return kTRUE;
+}
+
+//______________________________________________________________________________
+Bool_t AliTrackPointArray::GetPoint(AliTrackPoint &p, Int_t i) const
+{
+  // Get the point at position i
+  //
+  if (i >= fNPoints) return kFALSE;
+  p.SetXYZ(fX[i],fY[i],fZ[i],&fCov[6*i]);
+  p.SetVolumeID(fVolumeID[i]);
+  return kTRUE;
+}
+
+//______________________________________________________________________________
+Bool_t AliTrackPointArray::HasVolumeID(UShort_t volid) const
+{
+  // This method checks if the array
+  // has at least one hit in the detector
+  // volume defined by volid
+  Bool_t check = kFALSE;
+  for (Int_t ipoint = 0; ipoint < fNPoints; ipoint++)
+    if (fVolumeID[ipoint] == volid) check = kTRUE;
+
+  return check;
+}
+
+ClassImp(AliTrackPoint)
+
+//______________________________________________________________________________
+AliTrackPoint::AliTrackPoint()
+{
+  // Default constructor
+  //
+  fX = fY = fZ = 0;
+  fVolumeID = 0;
+  memset(fCov,0,6*sizeof(Float_t));
+}
+
+
+//______________________________________________________________________________
+AliTrackPoint::AliTrackPoint(Float_t x, Float_t y, Float_t z, const Float_t *cov, UShort_t volid)
+{
+  // Constructor
+  //
+  SetXYZ(x,y,z,cov);
+  SetVolumeID(volid);
+}
+
+//______________________________________________________________________________
+AliTrackPoint::AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid)
+{
+  // Constructor
+  //
+  SetXYZ(xyz[0],xyz[1],xyz[2],cov);
+  SetVolumeID(volid);
+}
+
+//______________________________________________________________________________
+AliTrackPoint::AliTrackPoint(const AliTrackPoint &p):
+  TObject(p)
+{
+  // Copy constructor
+  //
+  SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
+  SetVolumeID(p.fVolumeID);
+}
+
+//_____________________________________________________________________________
+AliTrackPoint &AliTrackPoint::operator =(const AliTrackPoint& p)
+{
+  // assignment operator
+  //
+  if(this==&p) return *this;
+  ((TObject *)this)->operator=(p);
+
+  SetXYZ(p.fX,p.fY,p.fZ,&(p.fCov[0]));
+  SetVolumeID(p.fVolumeID);
+
+  return *this;
+}
+
+//______________________________________________________________________________
+void AliTrackPoint::SetXYZ(Float_t x, Float_t y, Float_t z, const Float_t *cov)
+{
+  // Set XYZ coordinates and their cov matrix
+  //
+  fX = x;
+  fY = y;
+  fZ = z;
+  if (cov)
+    memcpy(fCov,cov,6*sizeof(Float_t));
+}
+
+//______________________________________________________________________________
+void AliTrackPoint::SetXYZ(const Float_t *xyz, const Float_t *cov)
+{
+  // Set XYZ coordinates and their cov matrix
+  //
+  SetXYZ(xyz[0],xyz[1],xyz[2],cov);
+}
+
+//______________________________________________________________________________
+void AliTrackPoint::GetXYZ(Float_t *xyz, Float_t *cov) const
+{
+  xyz[0] = fX;
+  xyz[1] = fY;
+  xyz[2] = fZ;
+  if (cov)
+    memcpy(cov,fCov,6*sizeof(Float_t));
+}
diff --git a/STEER/AliTrackPointArray.h b/STEER/AliTrackPointArray.h
new file mode 100644 (file)
index 0000000..7493625
--- /dev/null
@@ -0,0 +1,98 @@
+#ifndef ALITRACKPOINTARRAY_H
+#define ALITRACKPOINTARRAY_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//////////////////////////////////////////////////////////////////////////////
+//                          Class AliTrackPoint                             //
+//   This class represent a single track space-point.                       //
+//   It is used to access the points array defined in AliTrackPointArray.   //
+//   Note that the space point coordinates are given in the global frame.   //
+//                                                                          //
+//   cvetan.cheshkov@cern.ch 3/11/2005                                      //
+//////////////////////////////////////////////////////////////////////////////
+
+#include <TObject.h>
+
+class AliTrackPoint : public TObject {
+
+ public:
+  AliTrackPoint();
+  AliTrackPoint(Float_t x, Float_t y, Float_t z, const Float_t *cov, UShort_t volid);
+  AliTrackPoint(const Float_t *xyz, const Float_t *cov, UShort_t volid);
+  AliTrackPoint(const AliTrackPoint &p);
+  AliTrackPoint& operator= (const AliTrackPoint& p);
+  virtual ~AliTrackPoint() {}
+
+  // Multiplication with TGeoMatrix and distance between points (chi2) to be implemented
+
+  void     SetXYZ(Float_t x, Float_t y, Float_t z, const Float_t *cov = 0);
+  void     SetXYZ(const Float_t *xyz, const Float_t *cov = 0);
+  void     SetVolumeID(UShort_t volid) { fVolumeID = volid; }
+
+  Float_t  GetX() const { return fX; }
+  Float_t  GetY() const { return fY; }
+  Float_t  GetZ() const { return fZ; }
+  void     GetXYZ(Float_t *xyz, Float_t *cov = 0) const;
+  const Float_t *GetCov() const { return &fCov[0]; }
+  UShort_t GetVolumeID() const { return fVolumeID; }
+
+ private:
+
+  Float_t  fX;        // X coordinate
+  Float_t  fY;        // Y coordinate
+  Float_t  fZ;        // Z coordinate
+  Float_t  fCov[6];   // Cov matrix
+  UShort_t fVolumeID; // Volume ID
+
+  ClassDef(AliTrackPoint,1)
+};
+
+//////////////////////////////////////////////////////////////////////////////
+//                          Class AliTrackPointArray                        //
+//   This class contains the ESD track space-points which are used during   //
+//   the alignment procedures. Each space-point consist of 3 coordinates    //
+//   (and their errors) and the index of the sub-detector which contains    //
+//   the space-point.                                                       //
+//   cvetan.cheshkov@cern.ch 3/11/2005                                      //
+//////////////////////////////////////////////////////////////////////////////
+
+class AliTrackPointArray : public TObject {
+
+ public:
+
+  AliTrackPointArray();
+  AliTrackPointArray(Int_t npoints);
+  AliTrackPointArray(const AliTrackPointArray &array);
+  AliTrackPointArray& operator= (const AliTrackPointArray& array);
+  virtual ~AliTrackPointArray();
+
+  //  Bool_t    AddPoint(Int_t i, AliCluster *cl, UShort_t volid);
+  Bool_t    AddPoint(Int_t i, const AliTrackPoint *p);
+
+  Int_t     GetNPoints() const { return fNPoints; }
+  Int_t     GetCovSize() const { return fSize; }
+  Bool_t    GetPoint(AliTrackPoint &p, Int_t i) const;
+  // Getters for fast access to the coordinate arrays
+  const Float_t*  GetX() const { return &fX[0]; }
+  const Float_t*  GetY() const { return &fY[0]; }
+  const Float_t*  GetZ() const { return &fZ[0]; }
+  const Float_t*  GetCov() const { return &fCov[0]; }
+  const UShort_t* GetVolumeID() const { return &fVolumeID[0]; }
+
+  Bool_t    HasVolumeID(UShort_t volid) const;
+  
+ private:
+
+  Int_t     fNPoints;    // Number of space points
+  Float_t   *fX;         //[fNPoints] Array with space point X coordinates
+  Float_t   *fY;         //[fNPoints] Array with space point Y coordinates
+  Float_t   *fZ;         //[fNPoints] Array with space point Z coordinates
+  Int_t     fSize;       // Size of array with cov matrices = 6*N of points
+  Float_t   *fCov;       //[fSize] Array with space point coordinates cov matrix
+  UShort_t  *fVolumeID;  //[fNPoints] Array of space point volume IDs
+
+  ClassDef(AliTrackPointArray,1)
+};
+
+#endif
diff --git a/STEER/AliTrackResiduals.cxx b/STEER/AliTrackResiduals.cxx
new file mode 100644 (file)
index 0000000..4604fd6
--- /dev/null
@@ -0,0 +1,176 @@
+/**************************************************************************
+ * 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 base class for track residuals
+//
+//
+//-----------------------------------------------------------------
+
+#include "AliTrackResiduals.h"
+
+#include "AliAlignObj.h"
+#include "AliTrackPointArray.h"
+
+ClassImp(AliTrackResiduals)
+
+//_____________________________________________________________________________
+AliTrackResiduals::AliTrackResiduals():
+  fN(0),
+  fLast(0),
+  fIsOwner(kTRUE)
+{
+  // Default constructor
+  fAlignObj = 0x0;
+  fVolArray = fTrackArray = 0x0;
+}
+
+//_____________________________________________________________________________
+AliTrackResiduals::AliTrackResiduals(Int_t ntracks, AliAlignObj *alignobj):
+  fN(ntracks),
+  fLast(0),
+  fAlignObj(alignobj),
+  fIsOwner(kTRUE)
+{
+  // Constructor
+  if (ntracks > 0) {
+    fVolArray = new AliTrackPointArray*[ntracks];
+    fTrackArray = new AliTrackPointArray*[ntracks];
+    for (Int_t itrack = 0; itrack < ntracks; itrack++)
+      fVolArray[itrack] = fTrackArray[itrack] = 0x0;
+  }
+}
+
+//_____________________________________________________________________________
+AliTrackResiduals::AliTrackResiduals(const AliTrackResiduals &res):
+  TObject(res),
+  fN(res.fN),
+  fLast(res.fLast),
+  fAlignObj(res.fAlignObj),
+  fIsOwner(kTRUE)
+{
+  // Copy constructor
+  // By default the created copy owns the track point arrays
+  if (fN > 0) {
+    fVolArray = new AliTrackPointArray*[fN];
+    fTrackArray = new AliTrackPointArray*[fN];
+    for (Int_t itrack = 0; itrack < fN; itrack++)
+      {
+       if (res.fVolArray[itrack])
+         fVolArray[itrack] = new AliTrackPointArray(*res.fVolArray[itrack]);
+       else
+         fVolArray = 0x0;
+       if (res.fTrackArray[itrack])
+         fTrackArray[itrack] = new AliTrackPointArray(*res.fTrackArray[itrack]);
+       else
+         fTrackArray = 0x0;
+      }
+  }
+}
+
+//_____________________________________________________________________________
+AliTrackResiduals &AliTrackResiduals::operator =(const AliTrackResiduals& res)
+{
+  // assignment operator
+  // Does not copy the track point arrays
+  if(this==&res) return *this;
+  ((TObject *)this)->operator=(res);
+
+  fN = res.fN;
+  fLast = res.fLast;
+  fIsOwner = kFALSE;
+  fAlignObj = res.fAlignObj;
+
+  fVolArray = res.fVolArray;
+  fTrackArray = res.fTrackArray;
+
+  return *this;
+}
+
+//_____________________________________________________________________________
+AliTrackResiduals::~AliTrackResiduals()
+{
+  // Destructor
+  DeleteTrackPointArrays();
+}
+
+//_____________________________________________________________________________
+void AliTrackResiduals::SetNTracks(Int_t ntracks)
+{
+  // Set new size for the track point arrays.
+  // Delete the old arrays and allocate the
+  // new ones.
+  DeleteTrackPointArrays();
+
+  fN = ntracks;
+  fLast = 0;
+  fIsOwner = kTRUE;
+
+  if (ntracks > 0) {
+    fVolArray = new AliTrackPointArray*[ntracks];
+    fTrackArray = new AliTrackPointArray*[ntracks];
+    for (Int_t itrack = 0; itrack < ntracks; itrack++)
+      fVolArray[itrack] = fTrackArray[itrack] = 0x0;
+  }
+}
+
+//_____________________________________________________________________________
+Bool_t AliTrackResiduals::AddTrackPointArrays(AliTrackPointArray *volarray, AliTrackPointArray *trackarray)
+{
+  // Adds pair of track space point and
+  // track extrapolation point arrays
+  if (!fVolArray || !fTrackArray) return kFALSE;
+
+  if (fLast >= fN) return kFALSE;
+
+  fVolArray[fLast] = volarray;
+  fTrackArray[fLast] = trackarray;
+  fLast++;
+
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliTrackResiduals::GetTrackPointArrays(Int_t i, AliTrackPointArray* &volarray, AliTrackPointArray* &trackarray) const
+{
+  // Provide an access to a pair of track point arrays
+  // with given index
+  if (i >= fLast) {
+    volarray = trackarray = 0x0;
+    return kFALSE;
+  }
+  else {
+    volarray = fVolArray[i];
+    trackarray = fTrackArray[i];
+    return kTRUE;
+  }
+}
+
+//_____________________________________________________________________________
+void AliTrackResiduals::DeleteTrackPointArrays()
+{
+  // Deletes the track point arrays only in case
+  // the object is their owner.
+  // Called by the destructor and SetNTracks methods.
+  if (fIsOwner) {
+    for (Int_t itrack = 0; itrack < fLast; itrack++)
+      {
+       delete fVolArray[itrack];
+       delete fTrackArray[itrack];
+      }
+    delete [] fVolArray;
+    delete [] fTrackArray;
+  }
+}
diff --git a/STEER/AliTrackResiduals.h b/STEER/AliTrackResiduals.h
new file mode 100644 (file)
index 0000000..34d230c
--- /dev/null
@@ -0,0 +1,55 @@
+#ifndef ALITRACKRESIDUALS_H
+#define ALITRACKRESIDUALS_H
+
+//************************************************************************
+// AliTrackResiduals: base class for collecting the track space point    *
+// residuals produced by the fast track fitters (AliTrackFitter class).  *
+// It provides an interface to the arrays which contain the space points *
+// and track extrapolation points within the detector volume to be       *
+// aligned. The derived classes should implement method to analyze the   *
+// track residuals and minimize their sum in order to get the            *
+// AliAlignObj for the given detector volume.                            *
+//************************************************************************
+
+#include "TObject.h"
+
+class AliAlignObj;
+class AliTrackPointArray;
+
+class AliTrackResiduals : public TObject {
+
+ public:
+
+  AliTrackResiduals();
+  AliTrackResiduals(Int_t ntracks, AliAlignObj *alignobj);
+  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);
+
+  virtual Bool_t Minimize() = 0;
+
+  Int_t  GetNTracks() const { return fN; }
+  Int_t  GetNFilledTracks() const { return fLast; }
+  Bool_t GetTrackPointArrays(Int_t i, AliTrackPointArray* &volarray, AliTrackPointArray* &trackarray) const;
+  AliAlignObj *GetAlignObj() const { return fAlignObj; }
+
+ protected:
+
+  void DeleteTrackPointArrays();
+
+  Int_t              fN;            // Number of tracks
+  Int_t              fLast;         // Index of the last filled track arrays
+  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
+  Bool_t             fIsOwner;      // Track point arrays owned by the object
+
+  ClassDef(AliTrackResiduals,1)
+
+};
+
+#endif
diff --git a/STEER/AliTrackResidualsChi2.cxx b/STEER/AliTrackResidualsChi2.cxx
new file mode 100644 (file)
index 0000000..38bb693
--- /dev/null
@@ -0,0 +1,117 @@
+/**************************************************************************
+ * 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 Chi2 minimization
+//
+//-----------------------------------------------------------------
+
+#include <TMinuit.h>
+#include <TGeoMatrix.h>
+
+#include "AliAlignObj.h"
+#include "AliTrackPointArray.h"
+#include "AliTrackResidualsChi2.h"
+
+ClassImp(AliTrackResidualsChi2)
+
+//______________________________________________________________________________
+void Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+{
+  // This function is called by minuit
+  // The corresponding member method is called
+  // using SetObjectFit/GetObjectFit methods of TMinuit
+  AliTrackResidualsChi2* dummy = (AliTrackResidualsChi2 *)gMinuit->GetObjectFit();
+  dummy->Chi2(npar, gin, f, par, iflag);
+}
+
+//______________________________________________________________________________
+Bool_t AliTrackResidualsChi2::Minimize()
+{
+  // Implementation of Chi2 based minimization
+  // of track residuala sum
+  TMinuit *gMinuit = new TMinuit(6);  //initialize TMinuit
+  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);
+
+  // Now ready for minimization step
+  arglist[0] = 500;
+  arglist[1] = 1.;
+  gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
+
+  // Print results
+  Double_t amin,edm,errdef;
+  Int_t nvpar,nparx,icstat;
+  gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
+  gMinuit->mnprin(3,amin);
+
+  return kTRUE;
+}
+
+//______________________________________________________________________________
+void AliTrackResidualsChi2::Chi2(Int_t & /* npar */, Double_t * /* gin */, Double_t &f, Double_t *par, Int_t /* iflag */)
+{
+  // Chi2 function to be minimized
+  // Sums all the track residuals
+  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];
+
+  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);
+      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];
+    }
+  }
+  f = chi2;
+}
diff --git a/STEER/AliTrackResidualsChi2.h b/STEER/AliTrackResidualsChi2.h
new file mode 100644 (file)
index 0000000..7bb8aa0
--- /dev/null
@@ -0,0 +1,33 @@
+#ifndef ALITRACKRESIDUALSCHI2_H
+#define ALITRACKRESIDUALSCHI2_H
+
+//************************************************************************
+// AliTrackResidualsChi2: derived class (from AliTrackResiduals) which   *
+// implements a minimization of the track residuals based on chi2        *
+// approach.                                                             *
+//                                                                       *
+//************************************************************************
+
+#include "AliAlignObj.h"
+#include "AliTrackResiduals.h"
+
+class AliTrackResidualsChi2 : public AliTrackResiduals {
+
+ public:
+  AliTrackResidualsChi2():AliTrackResiduals() { }
+  AliTrackResidualsChi2(Int_t ntracks, AliAlignObj *alignobj):AliTrackResiduals(ntracks,alignobj) { }
+  AliTrackResidualsChi2(const AliTrackResidualsChi2 &res):AliTrackResiduals(res) { }
+  AliTrackResidualsChi2& operator= (const AliTrackResidualsChi2& res) { ((AliTrackResiduals *)this)->operator=(res); return *this; }
+  virtual ~AliTrackResidualsChi2() { }
+
+  Bool_t Minimize();
+
+  void   Chi2(Int_t & /* npar */, Double_t * /* gin */, Double_t &f, Double_t *par, Int_t /* iflag */);
+
+ protected:
+
+  ClassDef(AliTrackResidualsChi2,1)
+
+};
+
+#endif
index 57300f7..dfa53d5 100644 (file)
@@ -17,6 +17,7 @@ class TTree;
 class AliKalmanTrack;
 class AliESD;
 class AliMagF;
+class AliTrackPoint;
 
 class AliTracker : public TObject {
 public:
@@ -41,6 +42,8 @@ public:
   virtual Int_t LoadClusters(TTree *)=0;
   virtual void UnloadClusters()=0;
   virtual AliCluster *GetCluster(Int_t index) const=0;
+  //  virtual UShort_t GetVolumeID(Int_t index) {return 0;}
+  virtual Bool_t GetTrackPoint(Int_t /* index */ , AliTrackPoint& /* p */) { return kFALSE;}
   virtual void  UseClusters(const AliKalmanTrack *t, Int_t from=0) const;
   virtual void  CookLabel(AliKalmanTrack *t,Float_t wrong) const; 
   Double_t GetX() const {return fX;}
index 1b683e6..4a3cabe 100644 (file)
 #pragma link C++ class AliAlignObj+;
 #pragma link C++ class AliAlignObjAngles+;
 #pragma link C++ class AliAlignObjMatrix+;
+#pragma link C++ class AliTrackPointArray+;
+#pragma link C++ class AliTrackPoint+;
+#pragma link C++ class AliTrackFitter+;
+#pragma link C++ class AliTrackFitterRieman+;
+#pragma link C++ class AliTrackResiduals+;
+#pragma link C++ class AliTrackResidualsChi2+;
+#pragma link C++ class AliAlignmentTracks+;
 
 #pragma link C++ class  TTreeDataElement+;
 #pragma link C++ class  TTreeStream+;
index b8e45cc..0ae7178 100644 (file)
@@ -26,7 +26,11 @@ AliCDBPath.cxx AliCDBRunRange.cxx AliCDBManager.cxx\
 AliCDBStorage.cxx AliCDBLocal.cxx AliCDBDump.cxx AliCDBGrid.cxx\
 AliReconstructor.cxx AliDetectorEventHeader.cxx TTreeStream.cxx\
 AliAlignObj.cxx AliAlignObjAngles.cxx AliAlignObjMatrix.cxx \
-AliRieman.cxx
+AliRieman.cxx\
+AliTrackPointArray.cxx\
+AliTrackFitter.cxx AliTrackFitterRieman.cxx\
+AliTrackResiduals.cxx AliTrackResidualsChi2.cxx\
+AliAlignmentTracks.cxx
 
 HDRS:= $(SRCS:.cxx=.h)