- // 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));
+ /// Compute the multiple scattering covariance matrix
+ /// (assume that track parameters at each cluster are corrects)
+ AliDebug(1,"Enter ComputeMCSCovariances");
+
+ // Reset the size of the covariance matrix if needed
+ Int_t nClusters = GetNClusters();
+ if (mcsCovariances.GetNrows() != nClusters) mcsCovariances.ResizeTo(nClusters,nClusters);
+
+ // Define variables
+ Int_t nChambers = AliMUONConstants::NTrackingCh();
+ AliMUONTrackParam* trackParamAtCluster;
+ AliMUONTrackParam extrapTrackParam;
+ Int_t currentChamber = 0, expectedChamber = 0, size = 0;
+ Double_t *mcsAngle2 = new Double_t[2*nChambers];
+ Double_t *zMCS = new Double_t[2*nChambers];
+ Int_t *indices = new Int_t[2*nClusters];
+
+ // Compute multiple scattering dispersion angle at each chamber
+ // and save the z position where it is calculated
+ for (Int_t iCluster = 0; iCluster < nClusters; iCluster++) {
+ trackParamAtCluster = (AliMUONTrackParam*) fTrackParamAtCluster->UncheckedAt(iCluster);
+
+ // look for missing chambers if any
+ currentChamber = trackParamAtCluster->GetClusterPtr()->GetChamberId();
+ while (currentChamber > expectedChamber) {
+
+ // Save the z position where MCS dispersion is calculated
+ zMCS[size] = AliMUONConstants::DefaultChamberZ(expectedChamber);
+
+ // Do not take into account MCS in chambers prior the first cluster
+ if (iCluster > 0) {
+
+ // Get track parameters at missing chamber z
+ extrapTrackParam = *trackParamAtCluster;
+ AliMUONTrackExtrap::ExtrapToZ(&extrapTrackParam, zMCS[size]);
+
+ // Save multiple scattering dispersion angle in missing chamber
+ mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(extrapTrackParam,AliMUONConstants::ChamberThicknessInX0(expectedChamber),1.);
+
+ } else mcsAngle2[size] = 0.;
+
+ expectedChamber++;
+ size++;
+ }
+
+ // Save z position where MCS dispersion is calculated
+ zMCS[size] = trackParamAtCluster->GetZ();
+
+ // Save multiple scattering dispersion angle in current chamber
+ mcsAngle2[size] = AliMUONTrackExtrap::GetMCSAngle2(*trackParamAtCluster,AliMUONConstants::ChamberThicknessInX0(currentChamber),1.);
+
+ // Save indice in zMCS array corresponding to the current cluster
+ indices[iCluster] = size;
+
+ expectedChamber = currentChamber + 1;
+ size++;
+ }
+
+ // complete array of z if last cluster is on the last but one chamber
+ if (currentChamber != nChambers-1) zMCS[size++] = AliMUONConstants::DefaultChamberZ(nChambers-1);
+
+ // Compute the covariance matrix
+ for (Int_t iCluster1 = 0; iCluster1 < nClusters; iCluster1++) {
+
+ for (Int_t iCluster2 = iCluster1; iCluster2 < nClusters; iCluster2++) {
+
+ // Initialization to 0 (diagonal plus upper triangular part)
+ mcsCovariances(iCluster1,iCluster2) = 0.;
+
+ // Compute contribution from multiple scattering in upstream chambers
+ for (Int_t k = 0; k < indices[iCluster1]; k++) {
+ mcsCovariances(iCluster1,iCluster2) += (zMCS[indices[iCluster1]] - zMCS[k]) * (zMCS[indices[iCluster2]] - zMCS[k]) * mcsAngle2[k];