From d0bfce8dcb84eaa3e16b4590e56db6ef496a427c Mon Sep 17 00:00:00 2001 From: gosset Date: Mon, 8 Jan 2001 11:01:02 +0000 Subject: [PATCH] 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 --- MUON/AliMUONEventReconstructor.cxx | 110 +++++++++++++++++++++++++---- MUON/AliMUONEventReconstructor.h | 7 ++ MUON/AliMUONSegment.cxx | 48 ++++++++----- MUON/AliMUONSegment.h | 24 +++++++ MUON/AliMUONTrack.cxx | 41 ++++++++--- 5 files changed, 192 insertions(+), 38 deletions(-) diff --git a/MUON/AliMUONEventReconstructor.cxx b/MUON/AliMUONEventReconstructor.cxx index fb97ac19247..67bd12794f0 100644 --- a/MUON/AliMUONEventReconstructor.cxx +++ b/MUON/AliMUONEventReconstructor.cxx @@ -15,6 +15,11 @@ /* $Log$ +Revision 1.19 2000/11/20 11:24:10 gosset +New package for reconstructed tracks (A. Gheata): +* tree output in the file "tree_reco.root" +* display events and make histograms from this file + Revision 1.18 2000/10/26 12:47:03 gosset Real distance between chambers of each station taken into account for the reconstruction parameters "fSegmentMaxDistBending[5]" @@ -139,6 +144,8 @@ Addition of files for track reconstruction in C++ //************* Defaults parameters for reconstruction static const Double_t kDefaultMinBendingMomentum = 3.0; +static const Double_t kDefaultMaxBendingMomentum = 500.0; +static const Double_t kDefaultMaxChi2 = 100.0; static const Double_t kDefaultMaxSigma2Distance = 16.0; static const Double_t kDefaultBendingResolution = 0.01; static const Double_t kDefaultNonBendingResolution = 0.144; @@ -239,6 +246,8 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void) // Set reconstruction parameters to default values // Would be much more convenient with a structure (or class) ???? fMinBendingMomentum = kDefaultMinBendingMomentum; + fMaxBendingMomentum = kDefaultMaxBendingMomentum; + fMaxChi2 = kDefaultMaxChi2; fMaxSigma2Distance = kDefaultMaxSigma2Distance; AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); @@ -251,14 +260,10 @@ void AliMUONEventReconstructor::SetReconstructionParametersToDefaults(void) if (ch < 4) fRMin[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) * 2.0 * TMath::Pi() / 180.0; else fRMin[ch] = 30.0; + // maximum radius at 10 degrees and Z of chamber + fRMax[ch] = TMath::Abs((&(pMUON->Chamber(ch)))->Z()) * + 10.0 * TMath::Pi() / 180.0; } - // maximum radius - // like in TRACKF_STAT (10 degrees ????) - fRMax[0] = fRMax[1] = 91.5; - fRMax[2] = fRMax[3] = 122.5; - fRMax[4] = fRMax[5] = 158.3; - fRMax[6] = fRMax[7] = 260.0; - fRMax[8] = fRMax[9] = 260.0; // ******** Parameters for making segments // should be parametrized ???? @@ -601,7 +606,8 @@ AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* H bendCoor = Hit->Y(); nonBendCoor = Hit->X(); radius = TMath::Sqrt((bendCoor * bendCoor) + (nonBendCoor * nonBendCoor)); - if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL; + // This cut is not needed with a realistic chamber geometry !!!! +// if ((radius < fRMin[chamber]) || (radius > fRMax[chamber])) return NULL; // new AliMUONHitForRec from GEANT hit and increment number of AliMUONHitForRec's hitForRec = new ((*fHitsForRecPtr)[fNHitsForRec]) AliMUONHitForRec(Hit); fNHitsForRec++; @@ -609,6 +615,9 @@ AliMUONHitForRec* AliMUONEventReconstructor::NewHitForRecFromGEANT(AliMUONHit* H hitForRec->SetBendingCoor(bendCoor + gRandom->Gaus(0., fBendingResolution)); hitForRec->SetNonBendingCoor(nonBendCoor + gRandom->Gaus(0., fNonBendingResolution)); +// // !!!! without smearing +// hitForRec->SetBendingCoor(bendCoor); +// hitForRec->SetNonBendingCoor(nonBendCoor); // more information into HitForRec // resolution: angular effect to be added here ???? hitForRec->SetBendingReso2(fBendingResolution * fBendingResolution); @@ -795,7 +804,7 @@ void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station) AliMUONSegment *segment; Bool_t last2st; Double_t bendingSlope, distBend, distNonBend, extBendCoor, extNonBendCoor, - impactParam = 0., maxImpactParam = 0.; // =0 to avoid compilation warnings. + impactParam = 0., maxImpactParam = 0., minImpactParam = 0.; // =0 to avoid compilation warnings. AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? if (fPrintLevel >= 1) cout << "enter MakeSegmentsPerStation (0...) " << Station << endl; @@ -809,6 +818,9 @@ void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station) // maximum impact parameter (cm) according to fMinBendingMomentum... maxImpactParam = TMath::Abs(GetImpactParamFromBendingMomentum(fMinBendingMomentum)); + // minimum impact parameter (cm) according to fMaxBendingMomentum... + minImpactParam = + TMath::Abs(GetImpactParamFromBendingMomentum(fMaxBendingMomentum)); } else last2st = kFALSE; // extrapolation factor from Z of first chamber to Z of second chamber @@ -852,7 +864,8 @@ void AliMUONEventReconstructor::MakeSegmentsPerStation(Int_t Station) // Conditions "distBend" and "impactParam" correlated for these stations ???? if ((distBend < fSegmentMaxDistBending[Station]) && (distNonBend < fSegmentMaxDistNonBending[Station]) && - (!last2st || (impactParam < maxImpactParam))) { + (!last2st || (impactParam < maxImpactParam)) && + (!last2st || (impactParam > minImpactParam))) { // make new segment segment = new ((*fSegmentsPtr[Station])[segmentIndex]) AliMUONSegment(hit1Ptr, hit2Ptr); @@ -1061,13 +1074,14 @@ void AliMUONEventReconstructor::FollowTracks(void) { // Follow tracks in stations(1..) 3, 2 and 1 // too long: should be made more modular !!!! - AliMUONHitForRec *bestHit, *extrapHit, *hit; - AliMUONSegment *bestSegment, *extrapSegment, *segment; + AliMUONHitForRec *bestHit, *extrapHit, *extrapCorrHit, *hit; + AliMUONSegment *bestSegment, *extrapSegment, *extrapCorrSegment, *segment; AliMUONTrack *track, *nextTrack; AliMUONTrackParam *trackParam1, trackParam[2], trackParamVertex; // -1 to avoid compilation warnings Int_t ch = -1, chInStation, chBestHit = -1, iHit, iSegment, station, trackIndex; Double_t bestChi2, chi2, dZ1, dZ2, dZ3, maxSigma2Distance, mcsFactor; + Double_t bendingMomentum, chi2Norm = 0.; AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? // local maxSigma2Distance, for easy increase in testing maxSigma2Distance = fMaxSigma2Distance; @@ -1102,6 +1116,7 @@ void AliMUONEventReconstructor::FollowTracks(void) // extrapolation to station trackParam1->ExtrapToStation(station, trackParam); extrapSegment = new AliMUONSegment(); // empty segment + extrapCorrSegment = new AliMUONSegment(); // empty corrected segment // multiple scattering factor corresponding to one chamber // and momentum in bending plane (not total) mcsFactor = 0.0136 * trackParam1->GetInverseBendingMomentum(); @@ -1121,6 +1136,14 @@ void AliMUONEventReconstructor::FollowTracks(void) extrapSegment->UpdateFromStationTrackParam (trackParam, mcsFactor, dZ1, dZ2, dZ3, station, trackParam1->GetInverseBendingMomentum()); + // same thing for corrected segment + // better to use copy constructor, after checking that it works properly !!!! + extrapCorrSegment->SetBendingCoorReso2(fBendingResolution * fBendingResolution); + extrapCorrSegment-> + SetNonBendingCoorReso2(fNonBendingResolution * fNonBendingResolution); + extrapCorrSegment->UpdateFromStationTrackParam + (trackParam, mcsFactor, dZ1, dZ2, dZ3, station, + trackParam1->GetInverseBendingMomentum()); bestChi2 = 5.0; bestSegment = NULL; if (fPrintLevel >= 10) { @@ -1134,7 +1157,20 @@ void AliMUONEventReconstructor::FollowTracks(void) // multiple scattering ???? // separation in 2 functions: Segment and HitForRec ???? segment = (AliMUONSegment*) ((*fSegmentsPtr[station])[iSegment]); - chi2 = segment->NormalizedChi2WithSegment(extrapSegment, maxSigma2Distance); + // correction of corrected segment (fBendingCoor and fNonBendingCoor) + // according to real Z value of "segment" and slopes of "extrapSegment" + extrapCorrSegment-> + SetBendingCoor(extrapSegment->GetBendingCoor() + + extrapSegment->GetBendingSlope() * + (segment->GetHitForRec1()->GetZ() - + (&(pMUON->Chamber(2 * station)))->Z())); + extrapCorrSegment-> + SetNonBendingCoor(extrapSegment->GetNonBendingCoor() + + extrapSegment->GetNonBendingSlope() * + (segment->GetHitForRec1()->GetZ() - + (&(pMUON->Chamber(2 * station)))->Z())); + chi2 = segment-> + NormalizedChi2WithSegment(extrapCorrSegment, maxSigma2Distance); if (chi2 < bestChi2) { // update best Chi2 and Segment if better found bestSegment = segment; @@ -1159,6 +1195,7 @@ void AliMUONEventReconstructor::FollowTracks(void) // should consider all possibilities ???? // multiple scattering ???? do about like for extrapSegment !!!! extrapHit = new AliMUONHitForRec(); // empty hit + extrapCorrHit = new AliMUONHitForRec(); // empty corrected hit bestChi2 = 3.0; bestHit = NULL; if (fPrintLevel >= 10) { @@ -1176,6 +1213,14 @@ void AliMUONEventReconstructor::FollowTracks(void) // resolutions from "extrapSegment" extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2()); extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2()); + // same things for corrected hit + // better to use copy constructor, after checking that it works properly !!!! + extrapCorrHit-> + SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor()); + extrapCorrHit-> + SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor()); + extrapHit->SetBendingReso2(extrapSegment->GetBendingCoorReso2()); + extrapHit->SetNonBendingReso2(extrapSegment->GetNonBendingCoorReso2()); // Loop over hits in the chamber ch = 2 * station + chInStation; for (iHit = fIndexOfFirstHitForRecPerChamber[ch]; @@ -1183,6 +1228,18 @@ void AliMUONEventReconstructor::FollowTracks(void) fNHitsForRecPerChamber[ch]; iHit++) { hit = (AliMUONHitForRec*) ((*fHitsForRecPtr)[iHit]); + // correction of corrected hit (fBendingCoor and fNonBendingCoor) + // according to real Z value of "hit" and slopes of right "trackParam" + extrapCorrHit-> + SetBendingCoor((&(trackParam[chInStation]))->GetBendingCoor() + + (&(trackParam[chInStation]))->GetBendingSlope() * + (hit->GetZ() - + (&(trackParam[chInStation]))->GetZ())); + extrapCorrHit-> + SetNonBendingCoor((&(trackParam[chInStation]))->GetNonBendingCoor() + + (&(trackParam[chInStation]))->GetNonBendingSlope() * + (hit->GetZ() - + (&(trackParam[chInStation]))->GetZ())); // condition for hit not already in segment ???? chi2 = hit->NormalizedChi2WithHitForRec(extrapHit, maxSigma2Distance); if (chi2 < bestChi2) { @@ -1210,16 +1267,31 @@ void AliMUONEventReconstructor::FollowTracks(void) // and corresponding TrackHit's, ... track->Remove(); delete extrapSegment; + delete extrapCorrSegment; + delete extrapHit; + delete extrapCorrHit; break; // stop the search for this candidate: // exit from the loop over station } + delete extrapHit; + delete extrapCorrHit; } delete extrapSegment; + delete extrapCorrSegment; // Sort track hits according to increasing Z track->GetTrackHitsPtr()->Sort(); // Update track parameters at first track hit (smallest Z) trackParam1 = ((AliMUONTrackHit*) (track->GetTrackHitsPtr()->First()))->GetTrackParam(); + bendingMomentum = 0.; + if (TMath::Abs(trackParam1->GetInverseBendingMomentum()) > 0.) + bendingMomentum = TMath::Abs(1/(trackParam1->GetInverseBendingMomentum())); + // Track removed if bendingMomentum not in window [min, max] + if ((bendingMomentum < fMinBendingMomentum) || (bendingMomentum > fMaxBendingMomentum)) { + track->Remove(); + break; // stop the search for this candidate: + // exit from the loop over station + } // Track fit // with multiple Coulomb scattering if all stations if (station == 0) track->SetFitMCS(1); @@ -1228,6 +1300,18 @@ void AliMUONEventReconstructor::FollowTracks(void) track->SetFitNParam(5); // with 5 parameters (momentum and position) track->SetFitStart(1); // from parameters at first hit track->Fit(); + Double_t numberOfDegFree = (2.0 * track->GetNTrackHits() - 5); + if (numberOfDegFree > 0) { + chi2Norm = track->GetFitFMin() / numberOfDegFree; + } else { + chi2Norm = 1.e10; + } + // Track removed if normalized chi2 too high + if (chi2Norm > fMaxChi2) { + track->Remove(); + break; // stop the search for this candidate: + // exit from the loop over station + } if (fPrintLevel >= 10) { cout << "FollowTracks: track candidate(0..): " << trackIndex << " after fit from station(0..): " << station << " to 4" << endl; diff --git a/MUON/AliMUONEventReconstructor.h b/MUON/AliMUONEventReconstructor.h index a1d8f96ce01..5b9013fa62e 100644 --- a/MUON/AliMUONEventReconstructor.h +++ b/MUON/AliMUONEventReconstructor.h @@ -31,6 +31,10 @@ class AliMUONEventReconstructor : public TObject { // Get and Set, Set to defaults Double_t GetMinBendingMomentum(void) {return fMinBendingMomentum;} void SetMinBendingMomentum(Double_t MinBendingMomentum) {fMinBendingMomentum = MinBendingMomentum;} + Double_t GetMaxBendingMomentum(void) {return fMaxBendingMomentum;} + void SetMaxBendingMomentum(Double_t MaxBendingMomentum) {fMaxBendingMomentum = MaxBendingMomentum;} + Double_t GetMaxChi2(void) {return fMaxChi2;} + void SetMaxChi2(Double_t MaxChi2) {fMaxChi2 = MaxChi2;} Double_t GetMaxSigma2Distance(void) {return fMaxSigma2Distance;} void SetMaxSigma2Distance(Double_t MaxSigma2Distance) {fMaxSigma2Distance = MaxSigma2Distance;} Double_t GetBendingResolution(void) {return fBendingResolution;} @@ -84,6 +88,9 @@ class AliMUONEventReconstructor : public TObject { // Parameters for event reconstruction Double_t fMinBendingMomentum; // minimum value (GeV/c) of momentum in bending plane + // Parameters for event reconstruction + Double_t fMaxBendingMomentum; // maximum value (GeV/c) of momentum in bending plane + Double_t fMaxChi2; // maximum Chi2 per degree of Freedom Double_t fMaxSigma2Distance; // maximum square distance in units of the variance (maximum chi2) Double_t fRMin[kMaxMuonTrackingChambers]; // minimum radius (cm) Double_t fRMax[kMaxMuonTrackingChambers]; // maximum radius (cm) diff --git a/MUON/AliMUONSegment.cxx b/MUON/AliMUONSegment.cxx index 01671570e80..fcf369a5750 100644 --- a/MUON/AliMUONSegment.cxx +++ b/MUON/AliMUONSegment.cxx @@ -15,6 +15,11 @@ /* $Log$ +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:06:39 hristov Inline functions moved from *.cxx to *.h files instead of forward declarations @@ -253,18 +258,27 @@ void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, // The "road" is parametrized from the old reco_muon.F // with 8 cm between stations. AliMUONTrackParam *param0; - Double_t cReso2, sReso2; +// Double_t cReso2, sReso2; // parameters to define the widths of the searching roads in station 0,1,2 // width = p0 + p1/ (momentum)^2 // station number: 0 1 2 - static Double_t p0BendingCoor[3] = { 6.43e-2, 1.64e-2, 0.034 }; - static Double_t p1BendingCoor[3] = { 986., 821., 446. }; - static Double_t p0BendingSlope[3] = { 3.54e-6, 3.63e-6, 3.6e-6 }; - static Double_t p1BendingSlope[3] = { 4.49e-3, 4.8e-3, 0.011 }; - static Double_t p0NonBendingCoor[3] = { 4.66e-2, 4.83e-2, 0.049 }; - static Double_t p1NonBendingCoor[3] = { 1444., 866., 354. }; - static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 }; - static Double_t p1NonBendingSlope[3] = { 0., 0., 0. }; +// static Double_t p0BendingCoor[3] = { 6.43e-2, 1.64e-2, 0.034 }; +// static Double_t p1BendingCoor[3] = { 986., 821., 446. }; +// static Double_t p0BendingSlope[3] = { 3.54e-6, 3.63e-6, 3.6e-6 }; +// static Double_t p1BendingSlope[3] = { 4.49e-3, 4.8e-3, 0.011 }; +// static Double_t p0NonBendingCoor[3] = { 4.66e-2, 4.83e-2, 0.049 }; +// static Double_t p1NonBendingCoor[3] = { 1444., 866., 354. }; +// static Double_t p0NonBendingSlope[3] = { 6.14e-4, 6.49e-4, 6.85e-4 }; +// static Double_t p1NonBendingSlope[3] = { 0., 0., 0. }; + + static Double_t p0BendingCoor[3] = { 6.43e-2, 6.43e-2, 6.43e-2 }; + static Double_t p1BendingCoor[3] = { 986., 986., 986. }; + static Double_t p0BendingSlope[3] = { 3.6e-6, 3.6e-6, 3.6e-6 }; + static Double_t p1BendingSlope[3] = { 1.1e-2, 1.1e-2, 1.1e-2 }; + static Double_t p0NonBendingCoor[3] = { 0.049, 0.049, 0.049 }; + static Double_t p1NonBendingCoor[3] = { 1444., 1444., 1444. }; + static Double_t p0NonBendingSlope[3] = { 6.8e-4, 6.8e-4, 6.8e-4 }; + static Double_t p1NonBendingSlope[3] = { 0., 0., 0. }; param0 = &(TrackParam[0]); // OLD version @@ -296,15 +310,15 @@ void AliMUONSegment::UpdateFromStationTrackParam(AliMUONTrackParam *TrackParam, // because they are added in the functions "NormalizedChi2WithSegment" // and "NormalizedChi2WithHitForRec" // Bending plane - cReso2 = fBendingCoorReso2; - sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ; - fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2; - fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2; +// cReso2 = fBendingCoorReso2; +// sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ; + fBendingCoorReso2 = p0BendingCoor[Station] + p1BendingCoor[Station]*InverseMomentum*InverseMomentum ; // - cReso2 + fBendingSlopeReso2 = p0BendingSlope[Station] + p1BendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2; // Non bending plane - cReso2 = fNonBendingCoorReso2; - sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ; - fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum - cReso2; - fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum - sReso2; +// cReso2 = fNonBendingCoorReso2; +// sReso2 = (2. * cReso2 )/ (Dz3*Dz3) ; + fNonBendingCoorReso2 = p0NonBendingCoor[Station] + p1NonBendingCoor[Station]*InverseMomentum*InverseMomentum; // - cReso2; + fNonBendingSlopeReso2 = p0NonBendingSlope[Station] + p1NonBendingSlope[Station]*InverseMomentum*InverseMomentum; // - sReso2; return; } diff --git a/MUON/AliMUONSegment.h b/MUON/AliMUONSegment.h index fd5f2e6a77a..e7a6fbf537e 100644 --- a/MUON/AliMUONSegment.h +++ b/MUON/AliMUONSegment.h @@ -29,6 +29,30 @@ class AliMUONSegment : public TObject { AliMUONHitForRec* GetHitForRec2(void) { // Get fHitForRecPtr2 return fHitForRecPtr2;} + Double_t GetBendingCoor(void) { + // Get fBendingCoor + return fBendingCoor;} + void SetBendingCoor(Double_t BendingCoor) { + // Set fBendingCoor + fBendingCoor = BendingCoor;} + Double_t GetBendingSlope(void) { + // Get fBendingSlope + return fBendingSlope;} + void SetBendingSlope(Double_t BendingSlope) { + // Set fBendingSlope + fBendingSlope = BendingSlope;} + Double_t GetNonBendingCoor(void) { + // Get fNonBendingCoor + return fNonBendingCoor;} + void SetNonBendingCoor(Double_t NonBendingCoor) { + // Set fNonBendingCoor + fNonBendingCoor = NonBendingCoor;} + Double_t GetNonBendingSlope(void) { + // Get fNonBendingSlope + return fNonBendingSlope;} + void SetNonBendingSlope(Double_t NonBendingSlope) { + // Set fNonBendingSlope + fNonBendingSlope = NonBendingSlope;} Double_t GetBendingCoorReso2(void) { // Get fBendingCoorReso2 return fBendingCoorReso2;} diff --git a/MUON/AliMUONTrack.cxx b/MUON/AliMUONTrack.cxx index 83fc6e228c6..6750b083ad7 100644 --- a/MUON/AliMUONTrack.cxx +++ b/MUON/AliMUONTrack.cxx @@ -15,6 +15,10 @@ /* $Log$ +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. @@ -66,7 +70,7 @@ Addition of files for track reconstruction in C++ #include #include -#include +#include #include #include @@ -301,7 +305,7 @@ void AliMUONTrack::Fit() // choice of function to be minimized according to fFitMCS if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2); else fgFitter->SetFCN(TrackChi2MCS); - arg[0] = 1; + arg[0] = -1; fgFitter->ExecuteCommand("SET PRINT", arg, 1); // More printing !!!! // Parameters according to "fFitStart" // (should be a function to be used at every place where needed ????) @@ -319,13 +323,15 @@ void AliMUONTrack::Fit() trackParam->GetNonBendingSlope(), 0.001, -0.5, 0.5); if (fFitNParam == 5) { - // set last 2 Minuit parameters (no limits for the search: min=max=0) + // set last 2 Minuit parameters + // mandatory limits in Bending to avoid NaN values of parameters fgFitter->SetParameter(3, "X", trackParam->GetNonBendingCoor(), - 0.03, 0.0, 0.0); + 0.03, -500.0, 500.0); + // mandatory limits in non Bending to avoid NaN values of parameters fgFitter->SetParameter(4, "Y", trackParam->GetBendingCoor(), - 0.10, 0.0, 0.0); + 0.10, -500.0, 500.0); } // search without gradient calculation in the function fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0); @@ -517,8 +523,8 @@ void TrackChi2MCS(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *P Double_t hbc1, hbc2, pbc1, pbc2; Double_t hnbc1, hnbc2, pnbc1, pnbc2; Int_t numberOfHit = trackBeingFitted->GetNTrackHits(); - TMatrix *covBending = new TMatrix(numberOfHit, numberOfHit); - TMatrix *covNonBending = new TMatrix(numberOfHit, numberOfHit); + TMatrixD *covBending = new TMatrixD(numberOfHit, numberOfHit); + TMatrixD *covNonBending = new TMatrixD(numberOfHit, numberOfHit); Double_t *msa2 = new Double_t[numberOfHit]; // Predicted coordinates and multiple scattering angles are first calculated @@ -583,7 +589,23 @@ 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] <