+ // Performs a Riemann fit taking tilting pad correction into account
+ // The equation of a Riemann circle, where the y position is substituted by the
+ // measured y-position taking pad tilting into account, has to be transformed
+ // into a 4-dimensional hyperplane equation
+ // Riemann circle: (x-x0)^2 + (y-y0)^2 -R^2 = 0
+ // Measured y-Position: ymeas = y - tan(phiT)(zc - zt)
+ // zc: center of the pad row
+ // zt: z-position of the track
+ // The z-position of the track is assumed to be linear dependent on the x-position
+ // Transformed equation: a + b * u + c * t + d * v + e * w - 2 * (ymeas + tan(phiT) * zc) * t = 0
+ // Transformation: u = 2 * x * t
+ // v = 2 * tan(phiT) * t
+ // w = 2 * tan(phiT) * (x - xref) * t
+ // t = 1 / (x^2 + ymeas^2)
+ // Parameters: a = -1/y0
+ // b = x0/y0
+ // c = (R^2 -x0^2 - y0^2)/y0
+ // d = offset
+ // e = dz/dx
+ // If the offset respectively the slope in z-position is impossible, the parameters are fixed using
+ // results from the simple riemann fit. Afterwards the fit is redone.
+ // The curvature is calculated according to the formula:
+ // curv = a/(1 + b^2 + c*a) = 1/R
+ //
+ // Paramters: - Array of tracklets (connected to the track candidate)
+ // - Flag selecting the error definition
+ // Output: - Chi2 values of the track (in Parameter list)
+ //
+ TLinearFitter *fitter = GetTiltedRiemanFitter();
+ fitter->StoreData(kTRUE);
+ fitter->ClearPoints();
+ AliTRDLeastSquare zfitter;
+ AliTRDcluster *cl = 0x0;
+
+ AliTRDseedV1 work[kNPlanes], *tracklet = 0x0;
+ if(!tracklets){
+ for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+ if(!(tracklet = track->GetTracklet(ipl))) continue;
+ if(!tracklet->IsOK()) continue;
+ new(&work[ipl]) AliTRDseedV1(*tracklet);
+ }
+ tracklets = &work[0];
+ }
+
+ Double_t xref = CalculateReferenceX(tracklets);
+ Double_t x, y, z, t, tilt, dx, w, we;
+ Double_t uvt[4];
+ Int_t nPoints = 0;
+ // Containers for Least-square fitter
+ for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+ if(!tracklets[ipl].IsOK()) continue;
+ for(Int_t itb = 0; itb < AliTRDseedV1::kNTimeBins; itb++){
+ if(!(cl = tracklets[ipl].GetClusters(itb))) continue;
+ if (!tracklets[ipl].IsUsable(itb)) continue;
+ x = cl->GetX();
+ y = cl->GetY();
+ z = cl->GetZ();
+ tilt = tracklets[ipl].GetTilt();
+ dx = x - xref;
+ // 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 *= sigError ? TMath::Sqrt(cl->GetSigmaY2()) : 0.2;
+ fitter->AddPoint(uvt, w, we);
+ zfitter.AddPoint(&x, z, static_cast<Double_t>(TMath::Sqrt(cl->GetSigmaZ2())));
+ nPoints++;
+ }
+ }
+ if(fitter->Eval()) return 1.E10;
+
+ Double_t z0 = fitter->GetParameter(3);
+ Double_t dzdx = fitter->GetParameter(4);
+
+
+ // Linear fitter - not possible to make boundaries
+ // Do not accept non possible z and dzdx combinations
+ Bool_t accept = kTRUE;
+ Double_t zref = 0.0;
+ for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) {
+ if(!tracklets[iLayer].IsOK()) continue;
+ zref = z0 + dzdx * (tracklets[iLayer].GetX0() - xref);
+ if (TMath::Abs(tracklets[iLayer].GetZfit(0) - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0)
+ accept = kFALSE;
+ }
+ if (!accept) {
+ zfitter.Eval();
+ Double_t dzmf = zfitter.GetFunctionParameter(1);
+ Double_t zmf = zfitter.GetFunctionValue(&xref);
+ fitter->FixParameter(3, zmf);
+ fitter->FixParameter(4, dzmf);
+ fitter->Eval();
+ fitter->ReleaseParameter(3);
+ fitter->ReleaseParameter(4);
+ z0 = fitter->GetParameter(3); // = zmf ?
+ dzdx = fitter->GetParameter(4); // = dzmf ?
+ }
+
+ // Calculate Curvature
+ Double_t a = fitter->GetParameter(0);
+ Double_t b = fitter->GetParameter(1);
+ Double_t c = fitter->GetParameter(2);
+ Double_t y0 = 1. / a;
+ Double_t x0 = -b * y0;
+ Double_t tmp = y0*y0 + x0*x0 - c*y0;
+ if(tmp<=0.) return 1.E10;
+ Double_t R = TMath::Sqrt(tmp);
+ Double_t C = 1.0 + b*b - c*a;
+ if (C > 0.0) C = a / TMath::Sqrt(C);
+
+ // Calculate chi2 of the fit
+ Double_t chi2 = fitter->GetChisquare()/Double_t(nPoints);
+
+ // Update the tracklets
+ if(!track){
+ for(Int_t ip = 0; ip < kNPlanes; ip++) {
+ x = tracklets[ip].GetX0();
+ tmp = R*R-(x-x0)*(x-x0);
+ if(tmp <= 0.) continue;
+ tmp = TMath::Sqrt(tmp);
+
+ // y: R^2 = (x - x0)^2 + (y - y0)^2
+ // => y = y0 +/- Sqrt(R^2 - (x - x0)^2)
+ tracklets[ip].SetYref(0, y0 - (y0>0.?1.:-1)*tmp);
+ // => dy/dx = (x - x0)/Sqrt(R^2 - (x - x0)^2)
+ tracklets[ip].SetYref(1, (x - x0) / tmp);
+ tracklets[ip].SetZref(0, z0 + dzdx * (x - xref));
+ tracklets[ip].SetZref(1, dzdx);
+ tracklets[ip].SetC(C);
+ tracklets[ip].SetChi2(chi2);
+ }
+ }
+ //update track points array
+ if(np && points){
+ Float_t xyz[3];
+ for(int ip=0; ip<np; ip++){
+ points[ip].GetXYZ(xyz);
+ xyz[1] = TMath::Abs(xyz[0] - x0) > R ? 100. : y0 - (y0>0.?1.:-1.)*TMath::Sqrt((R-(xyz[0]-x0))*(R+(xyz[0]-x0)));
+ xyz[2] = z0 + dzdx * (xyz[0] - xref);
+ points[ip].SetXYZ(xyz);
+ }
+ }
+
+ return chi2;
+}
+
+
+//____________________________________________________________________
+Double_t AliTRDtrackerV1::FitKalman(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t up, Int_t np, AliTrackPoint *points)
+{
+// Kalman filter implementation for the TRD.
+// It returns the positions of the fit in the array "points"
+//
+// Author : A.Bercuci@gsi.de
+
+ // printf("Start track @ x[%f]\n", track->GetX());
+
+ //prepare marker points along the track
+ Int_t ip = np ? 0 : 1;
+ while(ip<np){
+ if((up?-1:1) * (track->GetX() - points[ip].GetX()) > 0.) break;
+ //printf("AliTRDtrackerV1::FitKalman() : Skip track marker x[%d] = %7.3f. Before track start ( %7.3f ).\n", ip, points[ip].GetX(), track->GetX());
+ ip++;
+ }
+ //if(points) printf("First marker point @ x[%d] = %f\n", ip, points[ip].GetX());
+
+
+ AliTRDseedV1 tracklet, *ptrTracklet = 0x0;
+
+ //Loop through the TRD planes
+ for (Int_t jplane = 0; jplane < kNPlanes; jplane++) {
+ // GET TRACKLET OR BUILT IT
+ Int_t iplane = up ? jplane : kNPlanes - 1 - jplane;
+ if(tracklets){
+ if(!(ptrTracklet = &tracklets[iplane])) continue;
+ }else{
+ if(!(ptrTracklet = track->GetTracklet(iplane))){
+ /*AliTRDtrackerV1 *tracker = 0x0;
+ if(!(tracker = dynamic_cast<AliTRDtrackerV1*>( AliTRDReconstructor::Tracker()))) continue;
+ ptrTracklet = new(&tracklet) AliTRDseedV1(iplane);
+ if(!tracker->MakeTracklet(ptrTracklet, track)) */
+ continue;
+ }
+ }
+ if(!ptrTracklet->IsOK()) continue;
+
+ Double_t x = ptrTracklet->GetX0();
+
+ while(ip < np){
+ //don't do anything if next marker is after next update point.
+ if((up?-1:1) * (points[ip].GetX() - x) - fgkMaxStep < 0) break;
+ if(((up?-1:1) * (points[ip].GetX() - track->GetX()) < 0) && !PropagateToX(*track, points[ip].GetX(), fgkMaxStep)) return -1.;
+
+ Double_t xyz[3]; // should also get the covariance
+ track->GetXYZ(xyz);
+ track->Global2LocalPosition(xyz, track->GetAlpha());
+ points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]);
+ ip++;
+ }
+ // printf("plane[%d] tracklet[%p] x[%f]\n", iplane, ptrTracklet, x);
+
+ // Propagate closer to the next update point
+ if(((up?-1:1) * (x - track->GetX()) + fgkMaxStep < 0) && !PropagateToX(*track, x + (up?-1:1)*fgkMaxStep, fgkMaxStep)) return -1.;
+
+ if(!AdjustSector(track)) return -1;
+ if(TMath::Abs(track->GetSnp()) > fgkMaxSnp) return -1;
+
+ //load tracklet to the tracker and the track
+/* Int_t index;
+ if((index = FindTracklet(ptrTracklet)) < 0){
+ ptrTracklet = SetTracklet(&tracklet);
+ index = fTracklets->GetEntriesFast()-1;
+ }
+ track->SetTracklet(ptrTracklet, index);*/
+
+
+ // register tracklet to track with tracklet creation !!
+ // PropagateBack : loaded tracklet to the tracker and update index
+ // RefitInward : update index
+ // MakeTrack : loaded tracklet to the tracker and update index
+ if(!tracklets) track->SetTracklet(ptrTracklet, -1);
+
+
+ //Calculate the mean material budget along the path inside the chamber
+ Double_t xyz0[3]; track->GetXYZ(xyz0);
+ Double_t alpha = track->GetAlpha();
+ Double_t xyz1[3], y, z;
+ if(!track->GetProlongation(x, y, z)) return -1;
+ xyz1[0] = x * TMath::Cos(alpha) - y * TMath::Sin(alpha);
+ xyz1[1] = +x * TMath::Sin(alpha) + y * TMath::Cos(alpha);
+ xyz1[2] = z;
+ if((xyz0[0] - xyz1[9] < 1e-3) && (xyz0[0] - xyz1[9] < 1e-3)) continue; // check wheter we are at the same global x position
+ Double_t param[7];
+ if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param) <=0.) break;
+ Double_t xrho = param[0]*param[4]; // density*length
+ Double_t xx0 = param[1]; // radiation length
+
+ //Propagate the track
+ track->PropagateTo(x, xx0, xrho);
+ if (!AdjustSector(track)) break;
+
+ //Update track
+ Double_t chi2 = track->GetPredictedChi2(ptrTracklet);
+ if(chi2<1e+10) track->Update(ptrTracklet, chi2);
+ if(!up) continue;
+
+ //Reset material budget if 2 consecutive gold
+ if(iplane>0 && track->GetTracklet(iplane-1) && ptrTracklet->GetN() + track->GetTracklet(iplane-1)->GetN() > 20) track->SetBudget(2, 0.);
+ } // end planes loop
+
+ // extrapolation
+ while(ip < np){
+ if(((up?-1:1) * (points[ip].GetX() - track->GetX()) < 0) && !PropagateToX(*track, points[ip].GetX(), fgkMaxStep)) return -1.;
+
+ Double_t xyz[3]; // should also get the covariance
+ track->GetXYZ(xyz);
+ track->Global2LocalPosition(xyz, track->GetAlpha());
+ points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]);
+ ip++;
+ }
+
+ return track->GetChi2();
+}
+
+//_________________________________________________________________________
+Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope, Double_t xref)
+{
+ //
+ // Calculates the chi2-value of the track in z-Direction including tilting pad correction.
+ // A linear dependence on the x-value serves as a model.
+ // The parameters are related to the tilted Riemann fit.
+ // Parameters: - Array of tracklets (AliTRDseedV1) related to the track candidate
+ // - the offset for the reference x
+ // - the slope
+ // - the reference x position
+ // Output: - The Chi2 value of the track in z-Direction
+ //
+ Float_t chi2Z = 0, nLayers = 0;
+ for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNlayer; iLayer++) {
+ if(!tracklets[iLayer].IsOK()) continue;
+ Double_t z = offset + slope * (tracklets[iLayer].GetX0() - xref);
+ chi2Z += TMath::Abs(tracklets[iLayer].GetZfit(0) - z);
+ nLayers++;
+ }
+ chi2Z /= TMath::Max((nLayers - 3.0),1.0);
+ return chi2Z;
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t maxStep)
+{
+ //
+ // Starting from current X-position of track <t> this function
+ // extrapolates the track up to radial position <xToGo>.
+ // Returns 1 if track reaches the plane, and 0 otherwise
+ //
+
+ const Double_t kEpsilon = 0.00001;
+
+ // Current track X-position
+ Double_t xpos = t.GetX();
+
+ // Direction: inward or outward
+ Double_t dir = (xpos < xToGo) ? 1.0 : -1.0;
+
+ while (((xToGo - xpos) * dir) > kEpsilon) {
+
+ Double_t xyz0[3];
+ Double_t xyz1[3];
+ Double_t param[7];
+ Double_t x;
+ Double_t y;
+ Double_t z;
+
+ // The next step size
+ Double_t step = dir * TMath::Min(TMath::Abs(xToGo-xpos),maxStep);
+
+ // Get the global position of the starting point
+ t.GetXYZ(xyz0);
+
+ // X-position after next step
+ x = xpos + step;
+
+ // Get local Y and Z at the X-position of the next step
+ if (!t.GetProlongation(x,y,z)) {
+ return 0; // No prolongation possible
+ }
+
+ // The global position of the end point of this prolongation step
+ xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha());
+ xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha());
+ xyz1[2] = z;
+
+ // Calculate the mean material budget between start and
+ // end point of this prolongation step
+ if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) return 0;
+
+ // Propagate the track to the X-position after the next step
+ if (!t.PropagateTo(x, param[1], param[0]*param[4])) return 0;
+
+ // Rotate the track if necessary
+ AdjustSector(&t);
+
+ // New track X-position
+ xpos = t.GetX();
+
+ }
+
+ return 1;
+
+}
+
+
+//_____________________________________________________________________________
+Int_t AliTRDtrackerV1::ReadClusters(TClonesArray* &array, TTree *clusterTree) const
+{
+ //
+ // Reads AliTRDclusters from the file.
+ // The names of the cluster tree and branches
+ // should match the ones used in AliTRDclusterizer::WriteClusters()