X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=MUON%2FAliMUONTrack.cxx;h=2a3e830045eaddd961c8c23a73b457f9e2d52698;hp=a3d421106037f9d80de24d8669022c2f2eb2528b;hb=a866ac60756cfa8bc56bfa641f69a548ca60182e;hpb=bbe35c2310063348d5d2898b5eca7afa83c38933 diff --git a/MUON/AliMUONTrack.cxx b/MUON/AliMUONTrack.cxx index a3d42110603..2a3e830045e 100644 --- a/MUON/AliMUONTrack.cxx +++ b/MUON/AliMUONTrack.cxx @@ -13,78 +13,20 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.8 2001/01/08 11:01:02 gosset -Modifications used for addendum to Dimuon TDR (JP Cussonneau): -*. MaxBendingMomentum to make both a segment and a track (default 500) -*. MaxChi2 per degree of freedom to make a track (default 100) -*. MinBendingMomentum used also to make a track - and not only a segment (default 3) -*. wider roads for track search in stations 1 to 3 -*. extrapolation to actual Z instead of Z(chamber) in FollowTracks -*. in track fit: - - limits on parameters X and Y (+/-500) - - covariance matrices in double precision - - normalization of covariance matrices before inversion - - suppression of Minuit printouts -*. correction against memory leak (delete extrapHit) in FollowTracks -*. RMax to 10 degrees with Z(chamber) instead of fixed values; - RMin and Rmax cuts suppressed in NewHitForRecFromGEANT, - because useless with realistic geometry - -Revision 1.7 2000/09/19 15:50:46 gosset -TrackChi2MCS function: covariance matrix better calculated, -taking into account missing planes... - -Revision 1.6 2000/07/20 12:45:27 gosset -New "EventReconstructor..." structure, - hopefully more adapted to tree/streamer. -"AliMUONEventReconstructor::RemoveDoubleTracks" - to keep only one track among similar ones. - -Revision 1.5 2000/07/18 16:04:06 gosset -AliMUONEventReconstructor package: -* a few minor modifications and more comments -* a few corrections - * right sign for Z of raw clusters - * right loop over chambers inside station - * symmetrized covariance matrix for measurements (TrackChi2MCS) - * right sign of charge in extrapolation (ExtrapToZ) - * right zEndAbsorber for Branson correction below 3 degrees -* use of TVirtualFitter instead of TMinuit for AliMUONTrack::Fit -* no parameter for AliMUONTrack::Fit() but more fit parameters in Track object - -Revision 1.4 2000/06/30 10:15:48 gosset -Changes to EventReconstructor...: -precision fit with multiple Coulomb scattering; -extrapolation to vertex with Branson correction in absorber (JPC) - -Revision 1.3 2000/06/25 13:23:28 hristov -stdlib.h needed for non-Linux compilation - -Revision 1.2 2000/06/15 07:58:48 morsch -Code from MUON-dev joined - -Revision 1.1.2.3 2000/06/12 10:11:34 morsch -Dummy copy constructor and assignment operator added - -Revision 1.1.2.2 2000/06/09 12:58:05 gosset -Removed comment beginnings in Log sections of .cxx files -Suppressed most violations of coding rules - -Revision 1.1.2.1 2000/06/07 14:44:53 gosset -Addition of files for track reconstruction in C++ -*/ - -//__________________________________________________________________________ -// -// Reconstructed track in ALICE dimuon spectrometer -//__________________________________________________________________________ +/* $Id$ */ -#include "AliMUONTrack.h" +/////////////////////////////////////////////////// +// +// Reconstructed track +// in +// ALICE +// dimuon +// spectrometer +// +/////////////////////////////////////////////////// -#include +#include // for cout +#include // for exit() #include #include @@ -95,14 +37,15 @@ Addition of files for track reconstruction in C++ #include "AliMUONEventReconstructor.h" #include "AliMUONHitForRec.h" #include "AliMUONSegment.h" +#include "AliMUONTrack.h" #include "AliMUONTrackHit.h" -#include - // Functions to be minimized with Minuit void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag); +void mnvertLocal(Double_t* a, Int_t l, Int_t m, Int_t n, Int_t& ifail); + Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit); ClassImp(AliMUONTrack) // Class implementation in ROOT context @@ -160,16 +103,33 @@ AliMUONTrack::~AliMUONTrack() } //__________________________________________________________________________ -AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack) +AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack):TObject(MUONTrack) { -// Dummy copy constructor + fEventReconstructor = new AliMUONEventReconstructor(*MUONTrack.fEventReconstructor); + fTrackParamAtVertex = MUONTrack.fTrackParamAtVertex; + fTrackHitsPtr = new TObjArray(*MUONTrack.fTrackHitsPtr); + fNTrackHits = MUONTrack.fNTrackHits; + fFitMCS = MUONTrack.fFitMCS; + fFitNParam = MUONTrack.fFitNParam; + fFitFMin = MUONTrack.fFitFMin; + fFitStart = MUONTrack.fFitStart; } //__________________________________________________________________________ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack) { -// Dummy assignment operator + if (this == &MUONTrack) return *this; + + fEventReconstructor = new AliMUONEventReconstructor(*MUONTrack.fEventReconstructor); + fTrackParamAtVertex = MUONTrack.fTrackParamAtVertex; + fTrackHitsPtr = new TObjArray(*MUONTrack.fTrackHitsPtr); + fNTrackHits = MUONTrack.fNTrackHits; + fFitMCS = MUONTrack.fFitMCS; + fFitNParam = MUONTrack.fFitNParam; + fFitFMin = MUONTrack.fFitFMin; + fFitStart = MUONTrack.fFitStart; + return *this; } //__________________________________________________________________________ @@ -250,12 +210,12 @@ void AliMUONTrack::SetFitStart(Int_t FitStart) } //__________________________________________________________________________ -AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) { +AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) const { // Get pointer to TrackParamAtFirstHit return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();} //__________________________________________________________________________ -void AliMUONTrack::RecursiveDump(void) +void AliMUONTrack::RecursiveDump(void) const { // Recursive dump of AliMUONTrack, i.e. with dump of TrackHit's and HitForRec's AliMUONTrackHit *trackHit; @@ -405,11 +365,12 @@ void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec) } //__________________________________________________________________________ -void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) +void AliMUONTrack::SetTrackParamAtHit(Int_t indexHit, AliMUONTrackParam *TrackParam) const { // Set track parameters at TrackHit with index "indexHit" // from the track parameters pointed to by "TrackParam". - AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]); + //PH AliMUONTrackHit* trackHit = (AliMUONTrackHit*) ((*fTrackHitsPtr)[indexHit]); + AliMUONTrackHit* trackHit = (AliMUONTrackHit*) (fTrackHitsPtr->At(indexHit)); trackHit->SetTrackParam(TrackParam); } @@ -457,7 +418,7 @@ void AliMUONTrack::SetTrackParamAtVertex() } //__________________________________________________________________________ -void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag) +void TrackChi2(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) { // Return the "Chi2" to be minimized with Minuit for track fitting, // with "NParam" parameters @@ -508,7 +469,7 @@ void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Para } //__________________________________________________________________________ -void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Param, Int_t Flag) +void TrackChi2MCS(Int_t &NParam, Double_t * /*Gradient*/, Double_t &Chi2, Double_t *Param, Int_t /*Flag*/) { // Return the "Chi2" to be minimized with Minuit for track fitting, // with "NParam" parameters @@ -534,7 +495,6 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P } AliMUONTrackHit *hit; - Bool_t goodDeterminant; Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3; Double_t z, z1, z2, z3; AliMUONTrackHit *hit1, *hit2, *hit3; @@ -607,50 +567,27 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P } } // for (hitNumber2 = hitNumber1;... } // for (hitNumber1 = 0;... - // Normalization of covariance matrices - Double_t normCovBending2 = covBending->E2Norm(); - Double_t normCovNonBending2 = covNonBending->E2Norm(); - (*covBending) *= 1/normCovBending2; - (*covNonBending) *= 1/normCovNonBending2; -// if (covBending->Determinant() < 1.e-33) { -// printf(" *** covBending *** \n"); -// covBending->Print(); -// printf(" *** covNonBending *** \n"); -// covNonBending->Print(); -// cout << " number of hits " << numberOfHit << endl; -// cout << "Momentum = " << 1/Param[0] <Determinant() != 0) { - covNonBending->Invert(); - } else { - goodDeterminant = kFALSE; - cout << "Warning in ChiMCS Determinant non Bending=0: " << endl; - } + // Inversion of covariance matrices + // with "mnvertLocal", local "mnvert" function of Minuit. + // One cannot use directly "mnvert" since "TVirtualFitter" does not know it. + // One will have to replace this local function by the right inversion function + // from a specialized Root package for symmetric positive definite matrices, + // when available!!!! + Int_t ifailBending; + mnvertLocal(&((*covBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, + ifailBending); + Int_t ifailNonBending; + mnvertLocal(&((*covNonBending)(0,0)), numberOfHit, numberOfHit, numberOfHit, + ifailNonBending); // It would be worth trying to calculate the inverse of the covariance matrix // only once per fit, since it cannot change much in principle, // and it would save a lot of computing time !!!! // Calculates Chi2 - if (goodDeterminant) { + if ((ifailBending == 0) && (ifailNonBending == 0)) { // with Multiple Scattering if inversion correct - // Inverse matrices without normalization - (*covBending) *= 1/normCovBending2; - (*covNonBending) *= 1/normCovNonBending2; for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor(); @@ -720,3 +657,103 @@ Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit) varMultipleScatteringAngle * varMultipleScatteringAngle; return varMultipleScatteringAngle; } + +//______________________________________________________________________________ + void mnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail) +{ +//*-*-*-*-*-*-*-*-*-*-*-*Inverts a symmetric matrix*-*-*-*-*-*-*-*-*-*-*-*-* +//*-* ========================== +//*-* inverts a symmetric matrix. matrix is first scaled to +//*-* have all ones on the diagonal (equivalent to change of units) +//*-* but no pivoting is done since matrix is positive-definite. +//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* + + // taken from TMinuit package of Root (l>=n) + // fVERTs, fVERTq and fVERTpp changed to localVERTs, localVERTq and localVERTpp + // Double_t localVERTs[n], localVERTq[n], localVERTpp[n]; + Double_t * localVERTs = new Double_t[n]; + Double_t * localVERTq = new Double_t[n]; + Double_t * localVERTpp = new Double_t[n]; + // fMaxint changed to localMaxint + Int_t localMaxint = n; + + /* System generated locals */ + Int_t aOffset; + + /* Local variables */ + Double_t si; + Int_t i, j, k, kp1, km1; + + /* Parameter adjustments */ + aOffset = l + 1; + a -= aOffset; + + /* Function Body */ + ifail = 0; + if (n < 1) goto L100; + if (n > localMaxint) goto L100; +//*-*- scale matrix by sqrt of diag elements + for (i = 1; i <= n; ++i) { + si = a[i + i*l]; + if (si <= 0) goto L100; + localVERTs[i-1] = 1 / TMath::Sqrt(si); + } + for (i = 1; i <= n; ++i) { + for (j = 1; j <= n; ++j) { + a[i + j*l] = a[i + j*l]*localVERTs[i-1]*localVERTs[j-1]; + } + } +//*-*- . . . start main loop . . . . + for (i = 1; i <= n; ++i) { + k = i; +//*-*- preparation for elimination step1 + if (a[k + k*l] != 0) localVERTq[k-1] = 1 / a[k + k*l]; + else goto L100; + localVERTpp[k-1] = 1; + a[k + k*l] = 0; + kp1 = k + 1; + km1 = k - 1; + if (km1 < 0) goto L100; + else if (km1 == 0) goto L50; + else goto L40; +L40: + for (j = 1; j <= km1; ++j) { + localVERTpp[j-1] = a[j + k*l]; + localVERTq[j-1] = a[j + k*l]*localVERTq[k-1]; + a[j + k*l] = 0; + } +L50: + if (k - n < 0) goto L51; + else if (k - n == 0) goto L60; + else goto L100; +L51: + for (j = kp1; j <= n; ++j) { + localVERTpp[j-1] = a[k + j*l]; + localVERTq[j-1] = -a[k + j*l]*localVERTq[k-1]; + a[k + j*l] = 0; + } +//*-*- elimination proper +L60: + for (j = 1; j <= n; ++j) { + for (k = j; k <= n; ++k) { a[j + k*l] += localVERTpp[j-1]*localVERTq[k-1]; } + } + } +//*-*- elements of left diagonal and unscaling + for (j = 1; j <= n; ++j) { + for (k = 1; k <= j; ++k) { + a[k + j*l] = a[k + j*l]*localVERTs[k-1]*localVERTs[j-1]; + a[j + k*l] = a[k + j*l]; + } + } + delete localVERTs; + delete localVERTq; + delete localVERTpp; + return; +//*-*- failure return +L100: + delete localVERTs; + delete localVERTq; + delete localVERTpp; + ifail = 1; +} /* mnvertLocal */ +