+//_____________________________________________________________________________
+void AliTRDtrackerV1::AliTRDLeastSquare::Reset(){
+ //
+ // Reset the fitter
+ //
+ memset(fParams, 0, sizeof(Double_t) * 2);
+ memset(fCovarianceMatrix, 0, sizeof(Double_t) * 3);
+ memset(fSums, 0, sizeof(Double_t) * 6);
+}
+
+///////////////////////////////////////////////////////
+// //
+// Resources of class AliTRDtrackFitterRieman //
+// //
+///////////////////////////////////////////////////////
+
+//_____________________________________________________________________________
+AliTRDtrackerV1::AliTRDtrackFitterRieman::AliTRDtrackFitterRieman():
+ fTrackFitter(NULL),
+ fZfitter(NULL),
+ fCovarPolY(NULL),
+ fCovarPolZ(NULL),
+ fXref(0.),
+ fSysClusterError(0.)
+{
+ //
+ // Default constructor
+ //
+ fZfitter = new AliTRDLeastSquare;
+ fCovarPolY = new TMatrixD(3,3);
+ fCovarPolZ = new TMatrixD(2,2);
+ memset(fTracklets, 0, sizeof(AliTRDseedV1 *) * 6);
+ memset(fParameters, 0, sizeof(Double_t) * 5);
+ memset(fSumPolY, 0, sizeof(Double_t) * 5);
+ memset(fSumPolZ, 0, sizeof(Double_t) * 2);
+}
+
+//_____________________________________________________________________________
+AliTRDtrackerV1::AliTRDtrackFitterRieman::~AliTRDtrackFitterRieman(){
+ //
+ // Destructor
+ //
+ if(fZfitter) delete fZfitter;
+ if(fCovarPolY) delete fCovarPolY;
+ if(fCovarPolZ) delete fCovarPolZ;
+}
+
+//_____________________________________________________________________________
+void AliTRDtrackerV1::AliTRDtrackFitterRieman::Reset(){
+ //
+ // Reset the Fitter
+ //
+ if(fTrackFitter){
+ fTrackFitter->StoreData(kTRUE);
+ fTrackFitter->ClearPoints();
+ }
+ if(fZfitter){
+ fZfitter->Reset();
+ }
+ fXref = 0.;
+ memset(fTracklets, 0, sizeof(AliTRDseedV1 *) * AliTRDgeometry::kNlayer);
+ memset(fParameters, 0, sizeof(Double_t) * 5);
+ memset(fSumPolY, 0, sizeof(Double_t) * 5);
+ memset(fSumPolZ, 0, sizeof(Double_t) * 2);
+ for(Int_t irow = 0; irow < fCovarPolY->GetNrows(); irow++)
+ for(Int_t icol = 0; icol < fCovarPolY->GetNcols(); icol++){
+ (*fCovarPolY)(irow, icol) = 0.;
+ if(irow < 2 && icol < 2)
+ (*fCovarPolZ)(irow, icol) = 0.;
+ }
+}
+
+//_____________________________________________________________________________
+void AliTRDtrackerV1::AliTRDtrackFitterRieman::SetTracklet(Int_t itr, AliTRDseedV1 *tracklet){
+ //
+ // Add tracklet into the fitter
+ //
+ if(itr >= AliTRDgeometry::kNlayer) return;
+ fTracklets[itr] = tracklet;
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDtrackerV1::AliTRDtrackFitterRieman::Eval(){
+ //
+ // Perform the fit
+ // 1. Apply linear transformation and store points in the fitter
+ // 2. Evaluate the fit
+ // 3. Check if the result of the fit in z-direction is reasonable
+ // if not
+ // 3a. Fix the parameters 3 and 4 with the results of a simple least
+ // square fit
+ // 3b. Redo the fit with the fixed parameters
+ // 4. Store fit results (parameters and errors)
+ //
+ if(!fTrackFitter){
+ return 1e10;
+ }
+ fXref = CalculateReferenceX();
+ for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++) UpdateFitters(fTracklets[il]);
+ if(!fTrackFitter->GetNpoints()) return 1e10;
+ // perform the fit
+ fTrackFitter->Eval();
+ fZfitter->Eval();
+ fParameters[3] = fTrackFitter->GetParameter(3);
+ fParameters[4] = fTrackFitter->GetParameter(4);
+ if(!CheckAcceptable(fParameters[3], fParameters[4])) {
+ fTrackFitter->FixParameter(3, fZfitter->GetFunctionValue(&fXref));
+ fTrackFitter->FixParameter(4, fZfitter->GetFunctionParameter(1));
+ fTrackFitter->Eval();
+ fTrackFitter->ReleaseParameter(3);
+ fTrackFitter->ReleaseParameter(4);
+ fParameters[3] = fTrackFitter->GetParameter(3);
+ fParameters[4] = fTrackFitter->GetParameter(4);
+ }
+ // Update the Fit Parameters and the errors
+ fParameters[0] = fTrackFitter->GetParameter(0);
+ fParameters[1] = fTrackFitter->GetParameter(1);
+ fParameters[2] = fTrackFitter->GetParameter(2);
+
+ // Prepare Covariance estimation
+ (*fCovarPolY)(0,0) = fSumPolY[0]; (*fCovarPolY)(1,1) = fSumPolY[2]; (*fCovarPolY)(2,2) = fSumPolY[4];
+ (*fCovarPolY)(1,0) = (*fCovarPolY)(0,1) = fSumPolY[1];
+ (*fCovarPolY)(2,0) = (*fCovarPolY)(0,2) = fSumPolY[2];
+ (*fCovarPolY)(2,1) = (*fCovarPolY)(1,2) = fSumPolY[3];
+ fCovarPolY->Invert();
+ (*fCovarPolZ)(0,0) = fSumPolZ[0]; (*fCovarPolZ)(1,1) = fSumPolZ[2];
+ (*fCovarPolZ)(1,0) = (*fCovarPolZ)(0,1) = fSumPolZ[1];
+ fCovarPolZ->Invert();
+ return fTrackFitter->GetChisquare() / fTrackFitter->GetNpoints();
+}
+
+//_____________________________________________________________________________
+void AliTRDtrackerV1::AliTRDtrackFitterRieman::UpdateFitters(AliTRDseedV1 * const tracklet){
+ //
+ // Does the transformations and updates the fitters
+ // The following transformation is applied
+ //
+ AliTRDcluster *cl = NULL;
+ Double_t x, y, z, dx, t, w, we, yerr, zerr;
+ Double_t uvt[4];
+ if(!tracklet || !tracklet->IsOK()) return;
+ Double_t tilt = tracklet->GetTilt();
+ for(Int_t itb = 0; itb < AliTRDseedV1::kNclusters; itb++){
+ if(!(cl = tracklet->GetClusters(itb))) continue;
+ if(!cl->IsInChamber()) continue;
+ if (!tracklet->IsUsable(itb)) continue;
+ x = cl->GetX();
+ y = cl->GetY();
+ z = cl->GetZ();
+ dx = x - fXref;
+ // Transformation
+ t = 1./(x*x + y*y);
+ uvt[0] = 2. * x * t;
+ uvt[1] = t;
+ uvt[2] = 2. * tilt * t;
+ uvt[3] = 2. * tilt * dx * t;
+ w = 2. * (y + tilt*z) * t;
+ // error definition changes for the different calls
+ we = 2. * t;
+ we *= TMath::Sqrt(cl->GetSigmaY2()+tilt*tilt*cl->GetSigmaZ2());
+ // Update sums for error calculation
+ yerr = 1./(TMath::Sqrt(cl->GetSigmaY2()) + fSysClusterError);
+ yerr *= yerr;
+ zerr = 1./cl->GetSigmaZ2();
+ for(Int_t ipol = 0; ipol < 5; ipol++){
+ fSumPolY[ipol] += yerr;
+ yerr *= x;
+ if(ipol < 3){
+ fSumPolZ[ipol] += zerr;
+ zerr *= x;
+ }
+ }
+ fTrackFitter->AddPoint(uvt, w, we);
+ fZfitter->AddPoint(&x, z, static_cast<Double_t>(TMath::Sqrt(cl->GetSigmaZ2())));
+ }
+}
+
+//_____________________________________________________________________________
+Bool_t AliTRDtrackerV1::AliTRDtrackFitterRieman::CheckAcceptable(Double_t offset, Double_t slope){
+ //
+ // Check whether z-results are acceptable
+ // Definition: Distance between tracklet fit and track fit has to be
+ // less then half a padlength
+ // Point of comparision is at the anode wire
+ //
+ Bool_t acceptablez = kTRUE;
+ Double_t zref = 0.0;
+ for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) {
+ if(!fTracklets[iLayer]->IsOK()) continue;
+ zref = offset + slope * (fTracklets[iLayer]->GetX0() - fXref);
+ if (TMath::Abs(fTracklets[iLayer]->GetZfit(0) - zref) > fTracklets[iLayer]->GetPadLength() * 0.5 + 1.0)
+ acceptablez = kFALSE;
+ }
+ return acceptablez;
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDtrackerV1::AliTRDtrackFitterRieman::GetYat(Double_t x) const {
+ //
+ // Calculate y position out of the track parameters
+ // y: R^2 = (x - x0)^2 + (y - y0)^2
+ // => y = y0 +/- Sqrt(R^2 - (x - x0)^2)
+ // R = Sqrt() = 1/Curvature
+ // => y = y0 +/- Sqrt(1/Curvature^2 - (x - x0)^2)
+ //
+ Double_t y = 0;
+ Double_t disc = (x * fParameters[0] + fParameters[1]);
+ disc = 1 - fParameters[0]*fParameters[2] + fParameters[1]*fParameters[1] - disc*disc;
+ if (disc >= 0) {
+ disc = TMath::Sqrt(disc);
+ y = (1.0 - disc) / fParameters[0];
+ }
+ return y;
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDtrackerV1::AliTRDtrackFitterRieman::GetZat(Double_t x) const {
+ //
+ // Return z position for a given x position
+ // Simple linear function
+ //
+ return fParameters[3] + fParameters[4] * (x - fXref);
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDtrackerV1::AliTRDtrackFitterRieman::GetDyDxAt(Double_t x) const {
+ //
+ // Calculate dydx at a given radial position out of the track parameters
+ // dy: R^2 = (x - x0)^2 + (y - y0)^2
+ // => y = +/- Sqrt(R^2 - (x - x0)^2) + y0
+ // => dy/dx = (x - x0)/Sqrt(R^2 - (x - x0)^2)
+ // Curvature: cr = 1/R = a/Sqrt(1 + b^2 - c*a)
+ // => dy/dx = (x - x0)/(1/(cr^2) - (x - x0)^2)
+ //
+ Double_t x0 = -fParameters[1] / fParameters[0];
+ Double_t curvature = GetCurvature();
+ Double_t dy = 0;
+ if (-fParameters[2] * fParameters[0] + fParameters[1] * fParameters[1] + 1 > 0) {
+ if (1.0/(curvature * curvature) - (x - x0) * (x - x0) > 0.0) {
+ Double_t yderiv = (x - x0) / TMath::Sqrt(1.0/(curvature * curvature) - (x - x0) * (x - x0));
+ if (fParameters[0] < 0) yderiv *= -1.0;
+ dy = yderiv;
+ }
+ }
+ return dy;
+}
+
+//_____________________________________________________________________________
+Double_t AliTRDtrackerV1::AliTRDtrackFitterRieman::GetCurvature() const {
+ //
+ // Calculate track curvature
+ //
+ //
+ Double_t curvature = 1.0 + fParameters[1]*fParameters[1] - fParameters[2]*fParameters[0];
+ if (curvature > 0.0)
+ curvature = fParameters[0] / TMath::Sqrt(curvature);
+ return curvature;
+}
+
+//_____________________________________________________________________________
+void AliTRDtrackerV1::AliTRDtrackFitterRieman::GetCovAt(Double_t x, Double_t *cov) const {
+ //
+ // Error Definition according to gauss error propagation
+ //
+ TMatrixD transform(3,3);
+ transform(0,0) = transform(1,1) = transform(2,2) = 1;
+ transform(0,1) = transform(1,2) = x;
+ transform(0,2) = x*x;
+ TMatrixD covariance(transform, TMatrixD::kMult, *fCovarPolY);
+ covariance *= transform.T();
+ cov[0] = covariance(0,0);
+ TMatrixD transformZ(2,2);
+ transformZ(0,0) = transformZ(1,1) = 1;
+ transformZ(0,1) = x;
+ TMatrixD covarZ(transformZ, TMatrixD::kMult, *fCovarPolZ);
+ covarZ *= transformZ.T();
+ cov[1] = covarZ(0,0);
+ cov[2] = 0;
+}
+
+//____________________________________________________________________
+Double_t AliTRDtrackerV1::AliTRDtrackFitterRieman::CalculateReferenceX(){
+ //
+ // Calculates the reference x-position for the tilted Rieman fit defined as middle
+ // of the stack (middle between layers 2 and 3). For the calculation all the tracklets
+ // are taken into account
+ //
+ // Parameters: - Array of tracklets(AliTRDseedV1)
+ //
+ // Output: - The reference x-position(Float_t)
+ //
+ Int_t nDistances = 0;
+ Float_t meanDistance = 0.;
+ Int_t startIndex = 5;
+ for(Int_t il =5; il > 0; il--){
+ if(fTracklets[il]->IsOK() && fTracklets[il -1]->IsOK()){
+ Float_t xdiff = fTracklets[il]->GetX0() - fTracklets[il -1]->GetX0();
+ meanDistance += xdiff;
+ nDistances++;
+ }
+ if(fTracklets[il]->IsOK()) startIndex = il;
+ }
+ if(fTracklets[0]->IsOK()) startIndex = 0;
+ if(!nDistances){
+ // We should normally never get here
+ Float_t xpos[2]; memset(xpos, 0, sizeof(Float_t) * 2);
+ Int_t iok = 0, idiff = 0;
+ // This attempt is worse and should be avoided:
+ // check for two chambers which are OK and repeat this without taking the mean value
+ // Strategy avoids a division by 0;
+ for(Int_t il = 5; il >= 0; il--){
+ if(fTracklets[il]->IsOK()){
+ xpos[iok] = fTracklets[il]->GetX0();
+ iok++;
+ startIndex = il;
+ }
+ if(iok) idiff++; // to get the right difference;
+ if(iok > 1) break;
+ }
+ if(iok > 1){
+ meanDistance = (xpos[0] - xpos[1])/idiff;
+ }
+ else{
+ // we have do not even have 2 layers which are OK? The we do not need to fit at all
+ return 331.;
+ }
+ }
+ else{
+ meanDistance /= nDistances;
+ }
+ return fTracklets[startIndex]->GetX0() + (2.5 - startIndex) * meanDistance - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
+}