From 98937d9305d2bea0af83b8af07963a8b5062caf2 Mon Sep 17 00:00:00 2001 From: hristov Date: Mon, 28 Nov 2005 15:51:35 +0000 Subject: [PATCH] Alignment framework (C.Cheshkov). More information is available in http://agenda.cern.ch/askArchive.php?base=agenda&categ=a057717&id=a057717s1t6/transparencies --- STEER/AliAlignObj.cxx | 27 +- STEER/AliAlignObj.h | 16 +- STEER/AliAlignmentTracks.cxx | 576 ++++++++++++++++++++++++++++++++ STEER/AliAlignmentTracks.h | 89 +++++ STEER/AliESDtrack.cxx | 63 +++- STEER/AliESDtrack.h | 16 +- STEER/AliReconstruction.cxx | 68 +++- STEER/AliReconstruction.h | 7 +- STEER/AliTrackFitter.cxx | 111 ++++++ STEER/AliTrackFitter.h | 55 +++ STEER/AliTrackFitterRieman.cxx | 402 ++++++++++++++++++++++ STEER/AliTrackFitterRieman.h | 46 +++ STEER/AliTrackPointArray.cxx | 230 +++++++++++++ STEER/AliTrackPointArray.h | 98 ++++++ STEER/AliTrackResiduals.cxx | 176 ++++++++++ STEER/AliTrackResiduals.h | 55 +++ STEER/AliTrackResidualsChi2.cxx | 117 +++++++ STEER/AliTrackResidualsChi2.h | 33 ++ STEER/AliTracker.h | 3 + STEER/STEERLinkDef.h | 7 + STEER/libSTEER.pkg | 6 +- 21 files changed, 2185 insertions(+), 16 deletions(-) create mode 100644 STEER/AliAlignmentTracks.cxx create mode 100644 STEER/AliAlignmentTracks.h create mode 100644 STEER/AliTrackFitter.cxx create mode 100644 STEER/AliTrackFitter.h create mode 100644 STEER/AliTrackFitterRieman.cxx create mode 100644 STEER/AliTrackFitterRieman.h create mode 100644 STEER/AliTrackPointArray.cxx create mode 100644 STEER/AliTrackPointArray.h create mode 100644 STEER/AliTrackResiduals.cxx create mode 100644 STEER/AliTrackResiduals.h create mode 100644 STEER/AliTrackResidualsChi2.cxx create mode 100644 STEER/AliTrackResidualsChi2.h diff --git a/STEER/AliAlignObj.cxx b/STEER/AliAlignObj.cxx index 737d467ed08..9212ac6626c 100644 --- a/STEER/AliAlignObj.cxx +++ b/STEER/AliAlignObj.cxx @@ -31,9 +31,34 @@ #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) diff --git a/STEER/AliAlignObj.h b/STEER/AliAlignObj.h index 693c1b77bbb..f4c3d6b179b 100644 --- a/STEER/AliAlignObj.h +++ b/STEER/AliAlignObj.h @@ -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 index 00000000000..1c8f2b64dfb --- /dev/null +++ b/STEER/AliAlignmentTracks.cxx @@ -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 + +#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 index 00000000000..5cf6a8a99ad --- /dev/null +++ b/STEER/AliAlignmentTracks.h @@ -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 + +#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 diff --git a/STEER/AliESDtrack.cxx b/STEER/AliESDtrack.cxx index 45ed5655b4a..0270c607c1e 100644 --- a/STEER/AliESDtrack.cxx +++ b/STEER/AliESDtrack.cxx @@ -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 { // diff --git a/STEER/AliESDtrack.h b/STEER/AliESDtrack.h index 6650076619e..3fba2454498 100644 --- a/STEER/AliESDtrack.h +++ b/STEER/AliESDtrack.h @@ -28,6 +28,7 @@ #include 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 diff --git a/STEER/AliReconstruction.cxx b/STEER/AliReconstruction.cxx index 93d9bad67c0..6333f977d16 100644 --- a/STEER/AliReconstruction.cxx +++ b/STEER/AliReconstruction.cxx @@ -137,7 +137,7 @@ #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<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(); +} diff --git a/STEER/AliReconstruction.h b/STEER/AliReconstruction.h index ba2f126aed8..157b1b83efc 100644 --- a/STEER/AliReconstruction.h +++ b/STEER/AliReconstruction.h @@ -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 index 00000000000..23853223e59 --- /dev/null +++ b/STEER/AliTrackFitter.cxx @@ -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 + +#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 index 00000000000..73caed152ef --- /dev/null +++ b/STEER/AliTrackFitter.h @@ -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 index 00000000000..793d77bbb1a --- /dev/null +++ b/STEER/AliTrackFitterRieman.cxx @@ -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 index 00000000000..1c9c1fb48e1 --- /dev/null +++ b/STEER/AliTrackFitterRieman.h @@ -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 index 00000000000..10b03802ec0 --- /dev/null +++ b/STEER/AliTrackPointArray.cxx @@ -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 index 00000000000..74936259f72 --- /dev/null +++ b/STEER/AliTrackPointArray.h @@ -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 + +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 index 00000000000..4604fd6b796 --- /dev/null +++ b/STEER/AliTrackResiduals.cxx @@ -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 index 00000000000..34d230c6da4 --- /dev/null +++ b/STEER/AliTrackResiduals.h @@ -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 index 00000000000..38bb693d863 --- /dev/null +++ b/STEER/AliTrackResidualsChi2.cxx @@ -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 +#include + +#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 index 00000000000..7bb8aa08966 --- /dev/null +++ b/STEER/AliTrackResidualsChi2.h @@ -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 diff --git a/STEER/AliTracker.h b/STEER/AliTracker.h index 57300f7d522..dfa53d5458e 100644 --- a/STEER/AliTracker.h +++ b/STEER/AliTracker.h @@ -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;} diff --git a/STEER/STEERLinkDef.h b/STEER/STEERLinkDef.h index 1b683e6599b..4a3cabe8ae8 100644 --- a/STEER/STEERLinkDef.h +++ b/STEER/STEERLinkDef.h @@ -107,6 +107,13 @@ #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+; diff --git a/STEER/libSTEER.pkg b/STEER/libSTEER.pkg index b8e45cc7c24..0ae7178436f 100644 --- a/STEER/libSTEER.pkg +++ b/STEER/libSTEER.pkg @@ -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) -- 2.39.3