#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)
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,
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;
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);
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)
};
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
#include "AliESDtrack.h"
#include "AliKalmanTrack.h"
+#include "AliTrackPointArray.h"
#include "AliLog.h"
ClassImp(AliESDtrack)
fRICHtheta(0),
fRICHphi(0),
fRICHdx(0),
- fRICHdy(0)
+ fRICHdy(0),
+ fPoints(0)
{
//
// The default ESD constructor
fRICHtheta(track.fRICHtheta),
fRICHphi(track.fRICHphi),
fRICHdx(track.fRICHdx),
- fRICHdy(track.fRICHdy)
+ fRICHdy(track.fRICHdy),
+ fPoints(track.fPoints)
{
//
//copy constructor
//
//printf("Delete track\n");
delete fITStrack;
- delete fTRDtrack;
+ delete fTRDtrack;
+ delete fPoints;
}
//_______________________________________________________________________
fRICHdx = 0;
fRICHdy = 0;
+ fPoints = 0;
}
//_______________________________________________________________________
Double_t AliESDtrack::GetMass() 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
{
//
#include <TVector3.h>
class AliKalmanTrack;
+class AliTrackPointArray;
const Int_t kNPlane = 6;
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;}
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;}
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 {
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
#include "AliDetectorTag.h"
#include "AliEventTag.h"
-
+#include "AliTrackPointArray.h"
ClassImp(AliReconstruction)
fRunLoader(NULL),
fRawReader(NULL),
- fVertexer(NULL)
+ fVertexer(NULL),
+
+ fWriteAlignmentData(kFALSE)
{
// create reconstruction object with default parameters
fRunLoader(NULL),
fRawReader(NULL),
- fVertexer(NULL)
+ fVertexer(NULL),
+
+ fWriteAlignmentData(rec.fWriteAlignmentData)
{
// copy constructor
}
}
+ // 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;
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();
+}
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);
void CreateTag(TFile* file);
//==========================================//
+ void WriteAlignmentData(AliESD* esd);
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
--- /dev/null
+/**************************************************************************
+ * 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;
+ }
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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));
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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;
+ }
+}
--- /dev/null
+#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
--- /dev/null
+/**************************************************************************
+ * 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;
+}
--- /dev/null
+#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
class AliKalmanTrack;
class AliESD;
class AliMagF;
+class AliTrackPoint;
class AliTracker : public TObject {
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;}
#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+;
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)
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff