+
+ //__________________________________________________________________________
+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->GetTrackReconstructor()->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 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 */
+