X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=MUON%2FAliMUONTrack.cxx;h=919c11fd5f8c43f407f3d58820985962117e5ab5;hp=8a74b212d6a77438a009dfbf5fd3c7481b4602bc;hb=e804eb46751d1af08efc9eec79c166d52ba04c6b;hpb=a9e2aefa97f1153d6f61e580a32d396156706b7b diff --git a/MUON/AliMUONTrack.cxx b/MUON/AliMUONTrack.cxx index 8a74b212d6a..919c11fd5f8 100644 --- a/MUON/AliMUONTrack.cxx +++ b/MUON/AliMUONTrack.cxx @@ -15,6 +15,64 @@ /* $Log$ +Revision 1.10 2001/04/09 12:25:09 gosset +Inversion of covariance matrices with local copy of TMinuit::mnvert, +for symmetric positive definite matrices, instead of TMatrixD::Invert + +Revision 1.9 2001/01/17 20:59:24 hristov +chPrev initialised + +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 @@ -36,33 +94,47 @@ Addition of files for track reconstruction in C++ #include #include -#include +#include +#include +#include +#include #include "AliMUONEventReconstructor.h" #include "AliMUONHitForRec.h" #include "AliMUONSegment.h" #include "AliMUONTrackHit.h" -static AliMUONTrack *trackBeingFitted; -static AliMUONTrackParam *trackParamBeingFitted; +#include -// Function to be minimized with Minuit +// 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 +TVirtualFitter* AliMUONTrack::fgFitter = NULL; + //__________________________________________________________________________ AliMUONTrack::AliMUONTrack(AliMUONSegment* BegSegment, AliMUONSegment* EndSegment, AliMUONEventReconstructor* EventReconstructor) { // Constructor from two Segment's fEventReconstructor = EventReconstructor; // link back to EventReconstructor - // memory allocation for the TClonesArray of reconstructed TrackHit's - fTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 10); + // memory allocation for the TObjArray of pointers to reconstructed TrackHit's + fTrackHitsPtr = new TObjArray(10); fNTrackHits = 0; AddSegment(BegSegment); // add hits from BegSegment AddSegment(EndSegment); // add hits from EndSegment fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z SetTrackParamAtVertex(); // set track parameters at vertex + // set fit conditions... + fFitMCS = 0; + fFitNParam = 3; + fFitStart = 1; + fFitFMin = -1.0; return; } @@ -71,40 +143,125 @@ AliMUONTrack::AliMUONTrack(AliMUONSegment* Segment, AliMUONHitForRec* HitForRec, { // Constructor from one Segment and one HitForRec fEventReconstructor = EventReconstructor; // link back to EventReconstructor - // memory allocation for the TClonesArray of reconstructed TrackHit's - fTrackHitsPtr = new TClonesArray("AliMUONTrackHit", 10); + // memory allocation for the TObjArray of pointers to reconstructed TrackHit's + fTrackHitsPtr = new TObjArray(10); fNTrackHits = 0; AddSegment(Segment); // add hits from Segment AddHitForRec(HitForRec); // add HitForRec fTrackHitsPtr->Sort(); // sort TrackHits according to increasing Z SetTrackParamAtVertex(); // set track parameters at vertex + // set fit conditions... + fFitMCS = 0; + fFitNParam = 3; + fFitStart = 1; + fFitFMin = -1.0; return; } + //__________________________________________________________________________ +AliMUONTrack::~AliMUONTrack() +{ + // Destructor + if (fTrackHitsPtr) { + delete fTrackHitsPtr; // delete the TObjArray of pointers to TrackHit's + fTrackHitsPtr = NULL; + } +} + + //__________________________________________________________________________ AliMUONTrack::AliMUONTrack (const AliMUONTrack& MUONTrack) { // Dummy copy constructor } + //__________________________________________________________________________ AliMUONTrack & AliMUONTrack::operator=(const AliMUONTrack& MUONTrack) { // Dummy assignment operator return *this; } -// Inline functions for Get and Set: inline removed because it does not work !!!! -AliMUONTrackParam* AliMUONTrack::GetTrackParamAtVertex(void) { - // Get pointer to fTrackParamAtVertex - return &fTrackParamAtVertex;} + //__________________________________________________________________________ +void AliMUONTrack::Remove() +{ + // Remove current track from array of tracks, + // and corresponding track hits from array of track hits. + // Compress the TClonesArray it belongs to. + AliMUONTrackHit *nextTrackHit; + AliMUONEventReconstructor *eventRec = this->fEventReconstructor; + TClonesArray *trackHitsPtr = eventRec->GetRecTrackHitsPtr(); + // Loop over all track hits of track + AliMUONTrackHit *trackHit = (AliMUONTrackHit*) fTrackHitsPtr->First(); + while (trackHit) { + nextTrackHit = (AliMUONTrackHit*) fTrackHitsPtr->After(trackHit); + // Remove TrackHit from event TClonesArray. + // Destructor is called, + // hence links between HitForRec's and TrackHit's are updated + trackHitsPtr->Remove(trackHit); + trackHit = nextTrackHit; + } + // Remove the track from event TClonesArray + // Destructor is called, + // hence space for TObjArray of pointers to TrackHit's is freed + eventRec->GetRecTracksPtr()->Remove(this); + // Number of tracks decreased by 1 + eventRec->SetNRecTracks(eventRec->GetNRecTracks() - 1); + // Compress event TClonesArray of Track's: + // this is essential to retrieve the TClonesArray afterwards + eventRec->GetRecTracksPtr()->Compress(); + // Compress event TClonesArray of TrackHit's: + // this is probably also essential to retrieve the TClonesArray afterwards + trackHitsPtr->Compress(); +} + + //__________________________________________________________________________ +void AliMUONTrack::SetFitMCS(Int_t FitMCS) +{ + // Set multiple Coulomb scattering option for track fit "fFitMCS" + // from "FitMCS" argument: 0 without, 1 with + if ((FitMCS == 0) || (FitMCS == 1)) fFitMCS = FitMCS; + // better implementation with enum(with, without) ???? + else { + cout << "ERROR in AliMUONTrack::SetFitMCS(FitMCS)" << endl; + cout << "FitMCS = " << FitMCS << " is neither 0 nor 1" << endl; + exit(0); + } + return; +} + + //__________________________________________________________________________ +void AliMUONTrack::SetFitNParam(Int_t FitNParam) +{ + // Set number of parameters for track fit "fFitNParam" from "FitNParam": + // 3 for momentum, 5 for momentum and position + if ((FitNParam == 3) || (FitNParam == 5)) fFitNParam = FitNParam; + else { + cout << "ERROR in AliMUONTrack::SetFitNParam(FitNParam)" << endl; + cout << "FitNParam = " << FitNParam << " is neither 3 nor 5" << endl; + exit(0); + } + return; +} + + //__________________________________________________________________________ +void AliMUONTrack::SetFitStart(Int_t FitStart) +{ + // Set multiple Coulomb scattering option for track fit "fFitStart" + // from "FitStart" argument: 0 without, 1 with + if ((FitStart == 0) || (FitStart == 1)) fFitStart = FitStart; + // better implementation with enum(vertex, firstHit) ???? + else { + cout << "ERROR in AliMUONTrack::SetFitStart(FitStart)" << endl; + cout << "FitStart = " << FitStart << " is neither 0 nor 1" << endl; + exit(0); + } + return; +} + + //__________________________________________________________________________ AliMUONTrackParam* AliMUONTrack::GetTrackParamAtFirstHit(void) { // Get pointer to TrackParamAtFirstHit return ((AliMUONTrackHit*) (fTrackHitsPtr->First()))->GetTrackParam();} -TClonesArray* AliMUONTrack::GetTrackHitsPtr(void) { - // Get fTrackHitsPtr - return fTrackHitsPtr;} -Int_t AliMUONTrack::GetNTrackHits(void) { - // Get fNTrackHits - return fNTrackHits;} //__________________________________________________________________________ void AliMUONTrack::RecursiveDump(void) @@ -129,77 +286,112 @@ void AliMUONTrack::RecursiveDump(void) } //__________________________________________________________________________ -void AliMUONTrack::Fit(AliMUONTrackParam *TrackParam, Int_t NParam) +Int_t AliMUONTrack::HitsInCommon(AliMUONTrack* Track) +{ + // Returns the number of hits in common + // between the current track ("this") + // and the track pointed to by "Track". + Int_t hitsInCommon = 0; + AliMUONTrackHit *trackHit1, *trackHit2; + // Loop over hits of first track + trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->First(); + while (trackHit1) { + // Loop over hits of second track + trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->First(); + while (trackHit2) { + // Increment "hitsInCommon" if both TrackHits point to the same HitForRec + if ( (trackHit1->GetHitForRecPtr()) == + (trackHit2->GetHitForRecPtr()) ) hitsInCommon++; + trackHit2 = (AliMUONTrackHit*) Track->GetTrackHitsPtr()->After(trackHit2); + } // trackHit2 + trackHit1 = (AliMUONTrackHit*) this->GetTrackHitsPtr()->After(trackHit1); + } // trackHit1 + return hitsInCommon; +} + + //__________________________________________________________________________ +void AliMUONTrack::Fit() { // Fit the current track ("this"), - // starting with track parameters pointed to by "TrackParam", - // and with 3 or 5 parameters ("NParam"): - // 3 if one keeps X and Y fixed in "TrackParam", - // 5 if one lets them vary. - if ((NParam != 3) && (NParam != 5)) { - cout << "ERROR in AliMUONTrack::Fit, NParam = " << NParam; - cout << " , i.e. neither 3 nor 5 ====> EXIT" << endl; - exit(0); // right instruction for exit ???? - } - Int_t error = 0; + // with or without multiple Coulomb scattering according to "fFitMCS", + // with the number of parameters given by "fFitNParam" + // (3 if one keeps X and Y fixed in "TrackParam", 5 if one lets them vary), + // starting, according to "fFitStart", + // with track parameters at vertex or at the first TrackHit. + // "fFitMCS", "fFitNParam" and "fFitStart" have to be set before + // by calling the corresponding Set methods. Double_t arg[1], benC, errorParam, invBenP, lower, nonBenC, upper, x, y; - TString parName; - TMinuit *minuit = new TMinuit(5); - trackBeingFitted = this; // for the track to be known from the function to minimize - trackParamBeingFitted = TrackParam; // for the track parameters to be known from the function to minimize; possible to use only Minuit parameters ???? - minuit->mninit(5, 10, 7); // sysrd, syswr, syssa: useful ???? + char parName[50]; + AliMUONTrackParam *trackParam; + // Check if Minuit is initialized... + fgFitter = TVirtualFitter::Fitter(this); // add 3 or 5 for the maximum number of parameters ??? + fgFitter->Clear(); // necessary ???? probably yes + // how to reset the printout number at every fit ???? + // is there any risk to leave it like that ???? // how to go faster ???? choice of Minuit parameters like EDM ???? - minuit->SetFCN(TrackChi2); - minuit->SetPrintLevel(1); // More printing !!!! - // set first 3 parameters (try with no limits for the search: min=max=0) - minuit->mnparm(0, "InvBenP", - TrackParam->GetInverseBendingMomentum(), - 0.003, 0.0, 0.0, error); - minuit->mnparm(1, "BenS", - TrackParam->GetBendingSlope(), - 0.001, 0.0, 0.0, error); - minuit->mnparm(2, "NonBenS", - TrackParam->GetNonBendingSlope(), - 0.001, 0.0, 0.0, error); - if (NParam == 5) { - // set last 2 parameters (try with no limits for the search: min=max=0) - minuit->mnparm(3, "X", - TrackParam->GetNonBendingCoor(), - 0.03, 0.0, 0.0, error); - minuit->mnparm(4, "Y", - TrackParam->GetBendingCoor(), - 0.10, 0.0, 0.0, error); + // choice of function to be minimized according to fFitMCS + if (fFitMCS == 0) fgFitter->SetFCN(TrackChi2); + else fgFitter->SetFCN(TrackChi2MCS); + 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 ????) + if (fFitStart == 0) trackParam = &fTrackParamAtVertex; + else trackParam = this->GetTrackParamAtFirstHit(); + // set first 3 Minuit parameters + // could be tried with no limits for the search (min=max=0) ???? + fgFitter->SetParameter(0, "InvBenP", + trackParam->GetInverseBendingMomentum(), + 0.003, -0.4, 0.4); + fgFitter->SetParameter(1, "BenS", + trackParam->GetBendingSlope(), + 0.001, -0.5, 0.5); + fgFitter->SetParameter(2, "NonBenS", + trackParam->GetNonBendingSlope(), + 0.001, -0.5, 0.5); + if (fFitNParam == 5) { + // set last 2 Minuit parameters + // mandatory limits in Bending to avoid NaN values of parameters + fgFitter->SetParameter(3, "X", + trackParam->GetNonBendingCoor(), + 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, -500.0, 500.0); } // search without gradient calculation in the function - minuit->mnexcm("SET NOGRADIENT", arg, 0, error); + fgFitter->ExecuteCommand("SET NOGRADIENT", arg, 0); // minimization - minuit->mnexcm("MINIMIZE", arg, 0, error); + fgFitter->ExecuteCommand("MINIMIZE", arg, 0); // exit from Minuit - minuit->mnexcm("EXIT", arg, 0, error); // necessary ???? - // print results - minuit->mnpout(0, parName, invBenP, errorParam, lower, upper, error); - minuit->mnpout(1, parName, benC, errorParam, lower, upper, error); - minuit->mnpout(2, parName, nonBenC, errorParam, lower, upper, error); - if (NParam == 5) { - minuit->mnpout(3, parName, x, errorParam, lower, upper, error); - minuit->mnpout(4, parName, y, errorParam, lower, upper, error); + fgFitter->ExecuteCommand("EXIT", arg, 0); // necessary ???? + // get results into "invBenP", "benC", "nonBenC" ("x", "y") + fgFitter->GetParameter(0, parName, invBenP, errorParam, lower, upper); + fgFitter->GetParameter(1, parName, benC, errorParam, lower, upper); + fgFitter->GetParameter(2, parName, nonBenC, errorParam, lower, upper); + if (fFitNParam == 5) { + fgFitter->GetParameter(3, parName, x, errorParam, lower, upper); + fgFitter->GetParameter(4, parName, y, errorParam, lower, upper); } // result of the fit into track parameters - TrackParam->SetInverseBendingMomentum(invBenP); - TrackParam->SetBendingSlope(benC); - TrackParam->SetNonBendingSlope(nonBenC); - if (NParam == 5) { - TrackParam->SetNonBendingCoor(x); - TrackParam->SetBendingCoor(y); + trackParam->SetInverseBendingMomentum(invBenP); + trackParam->SetBendingSlope(benC); + trackParam->SetNonBendingSlope(nonBenC); + if (fFitNParam == 5) { + trackParam->SetNonBendingCoor(x); + trackParam->SetBendingCoor(y); } - trackBeingFitted = NULL; - delete minuit; + // global result of the fit + Double_t fedm, errdef; + Int_t npari, nparx; + fgFitter->GetStats(fFitFMin, fedm, errdef, npari, nparx); } //__________________________________________________________________________ void AliMUONTrack::AddSegment(AliMUONSegment* Segment) { - // Add Segment + // Add Segment to the track AddHitForRec(Segment->GetHitForRec1()); // 1st hit AddHitForRec(Segment->GetHitForRec2()); // 2nd hit } @@ -207,8 +399,17 @@ void AliMUONTrack::AddSegment(AliMUONSegment* Segment) //__________________________________________________________________________ void AliMUONTrack::AddHitForRec(AliMUONHitForRec* HitForRec) { - // Add HitForRec - new ((*fTrackHitsPtr)[fNTrackHits]) AliMUONTrackHit(HitForRec); + // Add HitForRec to the track: + // actual TrackHit into TClonesArray of TrackHit's for the event; + // pointer to actual TrackHit in TObjArray of pointers to TrackHit's for the track + TClonesArray *recTrackHitsPtr = this->fEventReconstructor->GetRecTrackHitsPtr(); + Int_t eventTrackHits = this->fEventReconstructor->GetNRecTrackHits(); + // event + AliMUONTrackHit* trackHit = + new ((*recTrackHitsPtr)[eventTrackHits]) AliMUONTrackHit(HitForRec); + this->fEventReconstructor->SetNRecTrackHits(eventTrackHits + 1); + // track + fTrackHitsPtr->Add(trackHit); fNTrackHits++; } @@ -272,13 +473,18 @@ void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Para // and their current values in array pointed to by "Param". // Assumes that the track hits are sorted according to increasing Z. // Track parameters at each TrackHit are updated accordingly. + // Multiple Coulomb scattering is not taken into account + AliMUONTrack *trackBeingFitted; AliMUONTrackHit* hit; AliMUONTrackParam param1; Int_t hitNumber; Double_t zHit; Chi2 = 0.0; // initialize Chi2 // copy of track parameters to be fitted - param1 = *trackParamBeingFitted; + trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); + if (trackBeingFitted->GetFitStart() == 0) + param1 = *(trackBeingFitted->GetTrackParamAtVertex()); + else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit()); // Minuit parameters to be fitted into this copy param1.SetInverseBendingMomentum(Param[0]); param1.SetBendingSlope(Param[1]); @@ -309,3 +515,293 @@ void TrackChi2(Int_t &NParam, Double_t *Gradient, Double_t &Chi2, Double_t *Para dY * dY / hit->GetHitForRecPtr()->GetBendingReso2(); } } + + //__________________________________________________________________________ +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 + // and their current values in array pointed to by "Param". + // Assumes that the track hits are sorted according to increasing Z. + // Track parameters at each TrackHit are updated accordingly. + // Multiple Coulomb scattering is taken into account with covariance matrix. + AliMUONTrack *trackBeingFitted; + AliMUONTrackParam param1; + Chi2 = 0.0; // initialize Chi2 + // copy of track parameters to be fitted + trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); + if (trackBeingFitted->GetFitStart() == 0) + param1 = *(trackBeingFitted->GetTrackParamAtVertex()); + else param1 = *(trackBeingFitted->GetTrackParamAtFirstHit()); + // Minuit parameters to be fitted into this copy + param1.SetInverseBendingMomentum(Param[0]); + param1.SetBendingSlope(Param[1]); + param1.SetNonBendingSlope(Param[2]); + if (NParam == 5) { + param1.SetNonBendingCoor(Param[3]); + param1.SetBendingCoor(Param[4]); + } + + AliMUONTrackHit *hit; + Int_t chCurrent, chPrev = 0, hitNumber, hitNumber1, hitNumber2, hitNumber3; + Double_t z, z1, z2, z3; + AliMUONTrackHit *hit1, *hit2, *hit3; + Double_t hbc1, hbc2, pbc1, pbc2; + Double_t hnbc1, hnbc2, pnbc1, pnbc2; + Int_t numberOfHit = trackBeingFitted->GetNTrackHits(); + 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 + for (hitNumber = 0; hitNumber < numberOfHit; hitNumber++) { + hit = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber]; + z = hit->GetHitForRecPtr()->GetZ(); + // do something special if 2 hits with same Z ???? + // security against infinite loop ???? + (¶m1)->ExtrapToZ(z); // extrapolation + hit->SetTrackParam(¶m1); + // square of multiple scattering angle at current hit, with one chamber + msa2[hitNumber] = MultipleScatteringAngle2(hit); + // correction for eventual missing hits or multiple hits in a chamber, + // according to the number of chambers + // between the current hit and the previous one + chCurrent = hit->GetHitForRecPtr()->GetChamberNumber(); + if (hitNumber > 0) msa2[hitNumber] = msa2[hitNumber] * (chCurrent - chPrev); + chPrev = chCurrent; + } + + // Calculates the covariance matrix + for (hitNumber1 = 0; hitNumber1 < numberOfHit; hitNumber1++) { + hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; + z1 = hit1->GetHitForRecPtr()->GetZ(); + for (hitNumber2 = hitNumber1; hitNumber2 < numberOfHit; hitNumber2++) { + hit2 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2]; + z2 = hit2->GetHitForRecPtr()->GetZ(); + // initialization to 0 (diagonal plus upper triangular part) + (*covBending)(hitNumber2, hitNumber1) = 0.0; + // contribution from multiple scattering in bending plane: + // loop over upstream hits + for (hitNumber3 = 0; hitNumber3 < hitNumber1; hitNumber3++) { + hit3 = (AliMUONTrackHit*) + (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber3]; + z3 = hit3->GetHitForRecPtr()->GetZ(); + (*covBending)(hitNumber2, hitNumber1) = + (*covBending)(hitNumber2, hitNumber1) + + ((z1 - z3) * (z2 - z3) * msa2[hitNumber3]); + } + // equal contribution from multiple scattering in non bending plane + (*covNonBending)(hitNumber2, hitNumber1) = + (*covBending)(hitNumber2, hitNumber1); + if (hitNumber1 == hitNumber2) { + // Diagonal elements: add contribution from position measurements + // in bending plane + (*covBending)(hitNumber2, hitNumber1) = + (*covBending)(hitNumber2, hitNumber1) + + hit1->GetHitForRecPtr()->GetBendingReso2(); + // and in non bending plane + (*covNonBending)(hitNumber2, hitNumber1) = + (*covNonBending)(hitNumber2, hitNumber1) + + hit1->GetHitForRecPtr()->GetNonBendingReso2(); + } + else { + // Non diagonal elements: symmetrization + // for bending plane + (*covBending)(hitNumber1, hitNumber2) = + (*covBending)(hitNumber2, hitNumber1); + // and non bending plane + (*covNonBending)(hitNumber1, hitNumber2) = + (*covNonBending)(hitNumber2, hitNumber1); + } + } // for (hitNumber2 = hitNumber1;... + } // for (hitNumber1 = 0;... + + // 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 ((ifailBending == 0) && (ifailNonBending == 0)) { + // with Multiple Scattering if inversion correct + for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { + hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; + hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor(); + pbc1 = hit1->GetTrackParam()->GetBendingCoor(); + hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor(); + pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor(); + for (hitNumber2 = 0; hitNumber2 < numberOfHit; hitNumber2++) { + hit2 = (AliMUONTrackHit*) + (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber2]; + hbc2 = hit2->GetHitForRecPtr()->GetBendingCoor(); + pbc2 = hit2->GetTrackParam()->GetBendingCoor(); + hnbc2 = hit2->GetHitForRecPtr()->GetNonBendingCoor(); + pnbc2 = hit2->GetTrackParam()->GetNonBendingCoor(); + Chi2 = Chi2 + + ((*covBending)(hitNumber2, hitNumber1) * + (hbc1 - pbc1) * (hbc2 - pbc2)) + + ((*covNonBending)(hitNumber2, hitNumber1) * + (hnbc1 - pnbc1) * (hnbc2 - pnbc2)); + } + } + } else { + // without Multiple Scattering if inversion impossible + for (hitNumber1 = 0; hitNumber1 < numberOfHit ; hitNumber1++) { + hit1 = (AliMUONTrackHit*) (*(trackBeingFitted->GetTrackHitsPtr()))[hitNumber1]; + hbc1 = hit1->GetHitForRecPtr()->GetBendingCoor(); + pbc1 = hit1->GetTrackParam()->GetBendingCoor(); + hnbc1 = hit1->GetHitForRecPtr()->GetNonBendingCoor(); + pnbc1 = hit1->GetTrackParam()->GetNonBendingCoor(); + Chi2 = Chi2 + + ((hbc1 - pbc1) * (hbc1 - pbc1) / + hit1->GetHitForRecPtr()->GetBendingReso2()) + + ((hnbc1 - pnbc1) * (hnbc1 - pnbc1) / + hit1->GetHitForRecPtr()->GetNonBendingReso2()); + } + } + + delete covBending; + delete covNonBending; + delete [] msa2; +} + +Double_t MultipleScatteringAngle2(AliMUONTrackHit *TrackHit) +{ + // Returns square of multiple Coulomb scattering angle + // at TrackHit pointed to by "TrackHit" + Double_t slopeBending, slopeNonBending, radiationLength, inverseBendingMomentum2, inverseTotalMomentum2; + Double_t varMultipleScatteringAngle; + AliMUONTrack *trackBeingFitted = (AliMUONTrack*) AliMUONTrack::Fitter()->GetObjectFit(); + AliMUONTrackParam *param = TrackHit->GetTrackParam(); + // Better implementation in AliMUONTrack class ???? + slopeBending = param->GetBendingSlope(); + slopeNonBending = param->GetNonBendingSlope(); + // thickness in radiation length for the current track, + // taking local angle into account + radiationLength = + trackBeingFitted->GetEventReconstructor()->GetChamberThicknessInX0() * + TMath::Sqrt(1.0 + + slopeBending * slopeBending + slopeNonBending * slopeNonBending); + inverseBendingMomentum2 = + param->GetInverseBendingMomentum() * param->GetInverseBendingMomentum(); + inverseTotalMomentum2 = + inverseBendingMomentum2 * (1.0 + slopeBending * slopeBending) / + (1.0 + slopeBending *slopeBending + slopeNonBending * slopeNonBending); + varMultipleScatteringAngle = 0.0136 * (1.0 + 0.038 * TMath::Log(radiationLength)); + // The velocity is assumed to be 1 !!!! + varMultipleScatteringAngle = inverseTotalMomentum2 * radiationLength * + 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 a_offset; + + /* Local variables */ + Double_t si; + Int_t i, j, k, kp1, km1; + + /* Parameter adjustments */ + a_offset = l + 1; + a -= a_offset; + + /* 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 */ +