- // 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));