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