X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=blobdiff_plain;f=TRD%2FAliTRDtrackerV1.cxx;h=b4ad87e8d0f6628ffae46e6c8dff18ff4f3a0466;hp=d33d9d10ee1b75cfd0686ef4eb98eff7f69d1109;hb=29f3223b1e0086ed7b3c95c45250140736fa8330;hpb=9c87a076c9794c98dfda2398e0c3d421f26bf3e4 diff --git a/TRD/AliTRDtrackerV1.cxx b/TRD/AliTRDtrackerV1.cxx index d33d9d10ee1..b4ad87e8d0f 100644 --- a/TRD/AliTRDtrackerV1.cxx +++ b/TRD/AliTRDtrackerV1.cxx @@ -35,8 +35,11 @@ #include #include #include +#include +#include #include "AliLog.h" +#include "AliMathBase.h" #include "AliESDEvent.h" #include "AliGeomManager.h" #include "AliRieman.h" @@ -68,10 +71,12 @@ const Double_t AliTRDtrackerV1::fgkMaxChi2 = 12.0; // const Double_t AliTRDtrackerV1::fgkMaxSnp = 0.95; // Maximum local sine of the azimuthal angle const Double_t AliTRDtrackerV1::fgkMaxStep = 2.0; // Maximal step size in propagation Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = { - 0.1112, 0.1112, 0.1112, 0.0786, 0.0786, + 0.5112, 0.5112, 0.5112, 0.0786, 0.0786, 0.0786, 0.0786, 0.0579, 0.0579, 0.0474, 0.0474, 0.0408, 0.0335, 0.0335, 0.0335 -}; +}; +const Double_t AliTRDtrackerV1::fgkX0[kNPlanes] = { + 300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; Int_t AliTRDtrackerV1::fgNTimeBins = 0; AliRieman* AliTRDtrackerV1::fgRieman = 0x0; TLinearFitter* AliTRDtrackerV1::fgTiltedRieman = 0x0; @@ -81,7 +86,7 @@ TLinearFitter* AliTRDtrackerV1::fgTiltedRiemanConstrained = 0x0; AliTRDtrackerV1::AliTRDtrackerV1(AliTRDReconstructor *rec) :AliTracker() ,fReconstructor(0x0) - ,fGeom(new AliTRDgeometry()) + ,fGeom(0x0) ,fClusters(0x0) ,fTracklets(0x0) ,fTracks(0x0) @@ -90,19 +95,44 @@ AliTRDtrackerV1::AliTRDtrackerV1(AliTRDReconstructor *rec) // // Default constructor. // + + SetReconstructor(rec); // initialize reconstructor + + // initialize geometry + if(!AliGeomManager::GetGeometry()){ + AliFatal("Could not get geometry."); + } + fGeom = new AliTRDgeometry(); + fGeom->CreateClusterMatrixArray(); + TGeoHMatrix *matrix = 0x0; + Double_t loc[] = {0., 0., 0.}; + Double_t glb[] = {0., 0., 0.}; + for(Int_t ily=kNPlanes; ily--;){ + Int_t ism = 0; + while(!(matrix = fGeom->GetClusterMatrix(AliTRDgeometry::GetDetector(ily, 2, ism)))) ism++; + if(!matrix){ + AliError(Form("Could not get transformation matrix for layer %d. Use default.", ily)); + fR[ily] = fgkX0[ily]; + continue; + } + matrix->LocalToMaster(loc, glb); + fR[ily] = glb[0]+ AliTRDgeometry::AnodePos()-.5*AliTRDgeometry::AmThick() - AliTRDgeometry::DrThick(); + } + + // initialize calibration values AliTRDcalibDB *trd = 0x0; if (!(trd = AliTRDcalibDB::Instance())) { - AliFatal("Could not get calibration object"); + AliFatal("Could not get calibration."); } - if(!fgNTimeBins) fgNTimeBins = trd->GetNumberOfTimeBins(); + // initialize cluster containers for (Int_t isector = 0; isector < AliTRDgeometry::kNsector; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector); - for(Int_t isl =0; islGetDetector(); - + Int_t det = tracklet->GetDetector(); + Int_t sec = fGeom->GetSector(det); + Double_t alpha = (sec+.5)*AliTRDgeometry::GetAlpha(), + sinA = TMath::Sin(alpha), + cosA = TMath::Cos(alpha); Double_t local[3]; - local[0] = tracklet->GetX0(); - local[1] = tracklet->GetYfit(0); - local[2] = tracklet->GetZfit(0); + local[0] = tracklet->GetX(); + local[1] = tracklet->GetY(); + local[2] = tracklet->GetZ(); Double_t global[3]; - fGeom->RotateBack(idet, local, global); - p.SetXYZ(global[0],global[1],global[2]); - + fGeom->RotateBack(det, local, global); + + Double_t cov2D[3]; Float_t cov[6]; + tracklet->GetCovAt(local[0], cov2D); + cov[0] = cov2D[0]*sinA*sinA; + cov[1] =-cov2D[0]*sinA*cosA; + cov[2] =-cov2D[1]*sinA; + cov[3] = cov2D[0]*cosA*cosA; + cov[4] = cov2D[1]*cosA; + cov[5] = cov2D[2]; + // store the global position of the tracklet and its covariance matrix in the track point + p.SetXYZ(global[0],global[1],global[2], cov); // setting volume id - AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1; - switch (fGeom->GetLayer(idet)) { - case 0: - iLayer = AliGeomManager::kTRD1; - break; - case 1: - iLayer = AliGeomManager::kTRD2; - break; - case 2: - iLayer = AliGeomManager::kTRD3; - break; - case 3: - iLayer = AliGeomManager::kTRD4; - break; - case 4: - iLayer = AliGeomManager::kTRD5; - break; - case 5: - iLayer = AliGeomManager::kTRD6; - break; - }; - Int_t modId = fGeom->GetSector(idet) * fGeom->Nstack() + fGeom->GetStack(idet); + AliGeomManager::ELayerID iLayer = AliGeomManager::ELayerID(AliGeomManager::kTRD1+fGeom->GetLayer(det)); + Int_t modId = fGeom->GetSector(det) * AliTRDgeometry::kNstack + fGeom->GetStack(det); UShort_t volid = AliGeomManager::LayerToVolUID(iLayer, modId); p.SetVolumeID(volid); @@ -227,81 +249,117 @@ TLinearFitter* AliTRDtrackerV1::GetTiltedRiemanFitterConstraint() //____________________________________________________________________ AliRieman* AliTRDtrackerV1::GetRiemanFitter() { - if(!fgRieman) fgRieman = new AliRieman(AliTRDtrackingChamber::kNTimeBins * AliTRDgeometry::kNlayer); + if(!fgRieman) fgRieman = new AliRieman(AliTRDseedV1::kNtb * AliTRDgeometry::kNlayer); return fgRieman; } //_____________________________________________________________________________ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) { - // - // Gets seeds from ESD event. The seeds are AliTPCtrack's found and - // backpropagated by the TPC tracker. Each seed is first propagated - // to the TRD, and then its prolongation is searched in the TRD. - // If sufficiently long continuation of the track is found in the TRD - // the track is updated, otherwise it's stored as originaly defined - // by the TPC tracker. - // - - // Calibration monitor - AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance(); +// Propagation of ESD tracks from TPC to TOF detectors and building of the TRD track. For building +// a TRD track an ESD track is used as seed. The informations obtained on the TRD track (measured points, +// covariance, PID, etc.) are than used to update the corresponding ESD track. +// Each track seed is first propagated to the geometrical limit of the TRD detector. +// Its prolongation is searched in the TRD and if corresponding clusters are found tracklets are +// constructed out of them (see AliTRDseedV1::AttachClusters()) and the track is updated. +// Otherwise the ESD track is left unchanged. +// +// The following steps are performed: +// 1. Selection of tracks based on the variance in the y-z plane. +// 2. Propagation to the geometrical limit of the TRD volume. If track propagation fails the AliESDtrack::kTRDStop is set. +// 3. Prolongation inside the fiducial volume (see AliTRDtrackerV1::FollowBackProlongation()) and marking +// the following status bits: +// - AliESDtrack::kTRDin - if the tracks enters the TRD fiducial volume +// - AliESDtrack::kTRDStop - if the tracks fails propagation +// - AliESDtrack::kTRDbackup - if the tracks fulfills chi2 conditions and qualify for refitting +// 4. Writting to friends, PID, MC label, quality etc. Setting status bit AliESDtrack::kTRDout. +// 5. Propagation to TOF. If track propagation fails the AliESDtrack::kTRDStop is set. +// + + AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance(); // Calibration monitor if (!calibra) AliInfo("Could not get Calibra instance\n"); - Int_t found = 0; // number of tracks found + // Define scalers + Int_t nFound = 0, // number of tracks found + nSeeds = 0, // total number of ESD seeds + nTRDseeds= 0, // number of seeds in the TRD acceptance + nTPCseeds= 0; // number of TPC seeds Float_t foundMin = 20.0; Float_t *quality = 0x0; Int_t *index = 0x0; - Int_t nSeed = event->GetNumberOfTracks(); - if(nSeed){ - quality = new Float_t[nSeed]; - index = new Int_t[nSeed]; - for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) { + nSeeds = event->GetNumberOfTracks(); + // Sort tracks according to quality + // (covariance in the yz plane) + if(nSeeds){ + quality = new Float_t[nSeeds]; + index = new Int_t[nSeeds]; + for (Int_t iSeed = nSeeds; iSeed--;) { AliESDtrack *seed = event->GetTrack(iSeed); Double_t covariance[15]; seed->GetExternalCovariance(covariance); quality[iSeed] = covariance[0] + covariance[2]; } - // Sort tracks according to covariance of local Y and Z - TMath::Sort(nSeed, quality, index,kFALSE); + TMath::Sort(nSeeds, quality, index,kFALSE); } - // Backpropagate all seeds + // Propagate all seeds Int_t expectedClr; AliTRDtrackV1 track; - for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) { + for (Int_t iSeed = 0; iSeed < nSeeds; iSeed++) { // Get the seeds in sorted sequence AliESDtrack *seed = event->GetTrack(index[iSeed]); + Float_t p4 = seed->GetC(seed->GetBz()); // Check the seed status ULong_t status = seed->GetStatus(); if ((status & AliESDtrack::kTPCout) == 0) continue; if ((status & AliESDtrack::kTRDout) != 0) continue; - - // Do the back prolongation + + // Propagate to the entrance in the TRD mother volume new(&track) AliTRDtrackV1(*seed); - track.SetReconstructor(fReconstructor); + if(AliTRDgeometry::GetXtrdBeg() > (fgkMaxStep + track.GetX()) && !PropagateToX(track, AliTRDgeometry::GetXtrdBeg(), fgkMaxStep)){ + seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop); + continue; + } + if(!AdjustSector(&track)){ + seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop); + continue; + } + if(TMath::Abs(track.GetSnp()) > fgkMaxSnp) { + seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop); + continue; + } - //Int_t lbl = seed->GetLabel(); - //track.SetSeedLabel(lbl); + nTPCseeds++; - // Make backup and mark entrance in the TRD - seed->UpdateTrackParams(&track, AliESDtrack::kTRDin); + // store track status at TRD entrance seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup); - Float_t p4 = track.GetC(track.GetBz()); - expectedClr = FollowBackProlongation(track); - if (expectedClr<0) continue; // Back prolongation failed + // prepare track and do propagation in the TRD + track.SetReconstructor(fReconstructor); + track.SetKink(Bool_t(seed->GetKinkIndex(0))); + expectedClr = FollowBackProlongation(track); + // check if track entered the TRD fiducial volume + if(track.GetTrackLow()){ + seed->UpdateTrackParams(&track, AliESDtrack::kTRDin); + nTRDseeds++; + } + // check if track was stopped in the TRD + if (expectedClr<0){ + seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop); + continue; + } if(expectedClr){ - found++; + nFound++; // computes PID for track track.CookPID(); // update calibration references using this track if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(&track); // save calibration object - if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0 /*&& quality TODO*/){ + if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0){ AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track); calibTrack->SetOwner(); seed->AddCalibObject(calibTrack); @@ -314,22 +372,16 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) } if ((TMath::Abs(track.GetC(track.GetBz()) - p4) / TMath::Abs(p4) < 0.2) ||(track.Pt() > 0.8)) { - // + // Make backup for back propagation - // Int_t foundClr = track.GetNumberOfClusters(); if (foundClr >= foundMin) { - //AliInfo(Form("Making backup track ncls [%d]...", foundClr)); - //track.CookdEdx(); - //track.CookdEdxTimBin(seed->GetID()); track.CookLabel(1. - fgkLabelFraction); - if(track.GetBackupTrack()) UseClusters(track.GetBackupTrack()); + //if(track.GetBackupTrack()) UseClusters(track.GetBackupTrack()); // Sign only gold tracks if (track.GetChi2() / track.GetNumberOfClusters() < 4) { - if ((seed->GetKinkIndex(0) == 0) && (track.Pt() < 1.5)){ - //UseClusters(&track); - } + //if ((seed->GetKinkIndex(0) == 0) && (track.Pt() < 1.5)) UseClusters(&track); } Bool_t isGold = kFALSE; @@ -354,64 +406,54 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) isGold = kTRUE; } } - - //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) { - //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup); - //} } } - // Propagation to the TOF (I.Belikov) - if (track.IsStopped() == kFALSE) { - Double_t xtof = 371.0; - Double_t xTOF0 = 370.0; - - Double_t c2 = track.GetSnp() + track.GetC(track.GetBz()) * (xtof - track.GetX()); - if (TMath::Abs(c2) >= 0.99) continue; - - if (!PropagateToX(track, xTOF0, fgkMaxStep)) continue; - - // Energy losses taken to the account - check one more time - c2 = track.GetSnp() + track.GetC(track.GetBz()) * (xtof - track.GetX()); - if (TMath::Abs(c2) >= 0.99) continue; - - //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) { - // fHBackfit->Fill(7); - //delete track; - // continue; - //} - - Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha()); - Double_t y; - track.GetYAt(xtof,GetBz(),y); - if (y > ymax) { - if (!track.Rotate( AliTRDgeometry::GetAlpha())) continue; - }else if (y < -ymax) { - if (!track.Rotate(-AliTRDgeometry::GetAlpha())) continue; + // Propagation to the TOF + if(!(seed->GetStatus()&AliESDtrack::kTRDStop)) { + Int_t sm = track.GetSector(); + // default value in case we have problems with the geometry. + Double_t xtof = 371.; + //Calculate radial position of the beginning of the TOF + //mother volume. In order to avoid mixing of the TRD + //and TOF modules some hard values are needed. This are: + //1. The path to the TOF module. + //2. The width of the TOF (29.05 cm) + //(with the help of Annalisa de Caro Mar-17-2009) + if(gGeoManager){ + gGeoManager->cd(Form("/ALIC_1/B077_1/BSEGMO%d_1/BTOF%d_1", sm, sm)); + TGeoHMatrix *m = 0x0; + Double_t loc[]={0., 0., -.5*29.05}, glob[3]; + + if((m=gGeoManager->GetCurrentMatrix())){ + m->LocalToMaster(loc, glob); + xtof = TMath::Sqrt(glob[0]*glob[0]+glob[1]*glob[1]); + } } - - if (track.PropagateTo(xtof)) { - seed->UpdateTrackParams(&track, AliESDtrack::kTRDout); - track.UpdateESDtrack(seed); + if(xtof > (fgkMaxStep + track.GetX()) && !PropagateToX(track, xtof, fgkMaxStep)){ + seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop); + continue; } - } else { - if ((track.GetNumberOfClusters() > 15) && (track.GetNumberOfClusters() > 0.5*expectedClr)) { - seed->UpdateTrackParams(&track, AliESDtrack::kTRDout); - - track.UpdateESDtrack(seed); + if(!AdjustSector(&track)){ + seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop); + continue; } + if(TMath::Abs(track.GetSnp()) > fgkMaxSnp){ + seed->UpdateTrackParams(&track, AliESDtrack::kTRDStop); + continue; + } + seed->UpdateTrackParams(&track, AliESDtrack::kTRDout); + // TODO obsolete - delete + seed->SetTRDQuality(track.StatusForTOF()); } - - seed->SetTRDQuality(track.StatusForTOF()); seed->SetTRDBudget(track.GetBudget(0)); } if(index) delete [] index; if(quality) delete [] quality; - - AliInfo(Form("Number of seeds: %d", nSeed)); - AliInfo(Form("Number of back propagated TRD tracks: %d", found)); - + AliInfo(Form("Number of seeds: TPCout[%d] TRDin[%d]", nTPCseeds, nTRDseeds)); + AliInfo(Form("Number of tracks: TRDout[%d]", nFound)); + // run stand alone tracking if (fReconstructor->IsSeeding()) Clusters2Tracks(event); @@ -444,12 +486,11 @@ Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event) continue; } + // reject tracks which failed propagation in the TRD or + // are produced by the TRD stand alone tracker ULong_t status = seed->GetStatus(); - // reject tracks which failed propagation in the TRD - if((status & AliESDtrack::kTRDout) == 0) continue; - - // reject tracks which are produced by the TRD stand alone track finder. - if((status & AliESDtrack::kTRDin) == 0) continue; + if(!(status & AliESDtrack::kTRDout)) continue; + if(!(status & AliESDtrack::kTRDin)) continue; nseed++; track.ResetCovariance(50.0); @@ -460,16 +501,26 @@ Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event) if(FollowProlongation(track)){ // Prolongate to TPC if (PropagateToX(track, xTPC, fgkMaxStep)) { // -with update - seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit); - found++; - kUPDATE = kTRUE; + seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit); + found++; + kUPDATE = kTRUE; + } + + // Update the friend track + if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0){ + TObject *o = 0x0; Int_t ic = 0; + AliTRDtrackV1 *calibTrack = 0x0; + while((o = seed->GetCalibObject(ic++))){ + if(!(calibTrack = dynamic_cast(o))) continue; + calibTrack->SetTrackHigh(track.GetTrackHigh()); + } } - } + } // Prolongate to TPC without update if(!kUPDATE) { AliTRDtrackV1 tt(*seed); - if (PropagateToX(tt, xTPC, fgkMaxStep)) seed->UpdateTrackParams(&tt, AliESDtrack::kTRDrefit); + if (PropagateToX(tt, xTPC, fgkMaxStep)) seed->UpdateTrackParams(&tt, AliESDtrack::kTRDbackup); } } AliInfo(Form("Number of loaded seeds: %d",nseed)); @@ -505,15 +556,15 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) // Debug level 2 // + Bool_t kStoreIn = kTRUE; Int_t nClustersExpected = 0; - Int_t lastplane = 5; //GetLastPlane(&t); - for (Int_t iplane = lastplane; iplane >= 0; iplane--) { + for (Int_t iplane = kNPlanes; iplane--;) { Int_t index = 0; AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index); if(!tracklet) continue; if(!tracklet->IsOK()) AliWarning("tracklet not OK"); - Double_t x = tracklet->GetXref();//GetX0(); + Double_t x = tracklet->GetX();//GetX0(); // reject tracklets which are not considered for inward refit if(x > t.GetX()+fgkMaxStep) continue; @@ -551,9 +602,15 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) t.PropagateTo(x, xx0, xrho); if (!AdjustSector(&t)) break; } - - Double_t maxChi2 = t.GetPredictedChi2(tracklet); - if (maxChi2 < 1e+10 && t.Update(tracklet, maxChi2)){ + if(kStoreIn){ + t.SetTrackHigh(); + kStoreIn = kFALSE; + } + + Double_t cov[3]; tracklet->GetCovAt(x, cov); + Double_t p[2] = { tracklet->GetY(), tracklet->GetZ()}; + Double_t chi2 = ((AliExternalTrackParam)t).GetPredictedChi2(p, cov); + if (chi2 < 1e+10 && t.Update(p, cov, chi2)){ nClustersExpected += tracklet->GetN(); } } @@ -568,10 +625,12 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); + AliTRDtrackV1 track(t); + track.SetOwner(); cstreamer << "FollowProlongation" << "EventNumber=" << eventNumber << "ncl=" << nClustersExpected - //<< "track.=" << &t + << "track.=" << &track << "\n"; } @@ -582,32 +641,47 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) //_____________________________________________________________________________ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) { - // Extrapolates the TRD track in the TOF direction. - // - // Parameters - // t : the TRD track which has to be extrapolated - // - // Output - // number of clusters attached to the track - // - // Detailed description - // - // Starting from current radial position of track this function - // extrapolates the track through the 6 TRD layers. The following steps - // are being performed for each plane: - // 1. prepare track: - // a. get plane limits in the local x direction - // b. check crossing sectors - // c. check track inclination - // 2. build tracklet (see AliTRDseed::AttachClusters() for details) - // 3. evaluate material budget using the geo manager - // 4. propagate and update track using the tracklet information. - // - // Debug level 2 - // - - Int_t nClustersExpected = 0; - Double_t clength = .5*AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick(); +// Extrapolates/Build the TRD track in the TOF direction. +// +// Parameters +// t : the TRD track which has to be extrapolated +// +// Output +// number of clusters attached to the track +// +// Starting from current radial position of track this function +// extrapolates the track through the 6 TRD layers. The following steps +// are being performed for each plane: +// 1. Propagate track to the entrance of the next chamber: +// - get chamber limits in the radial direction +// - check crossing sectors +// - check track inclination +// - check track prolongation against boundary conditions (see exclusion boundaries on AliTRDgeometry::IsOnBoundary()) +// 2. Build tracklet (see AliTRDseed::AttachClusters() for details) for this layer if needed. If only +// Kalman filter is needed and tracklets are already linked to the track this step is skipped. +// 3. Fit tracklet using the information from the Kalman filter. +// 4. Propagate and update track at reference radial position of the tracklet. +// 5. Register tracklet with the tracker and track; update pulls monitoring. +// +// Observation +// 1. During the propagation a bit map is filled detailing the status of the track in each TRD chamber. The following errors are being registered for each tracklet: +// - AliTRDtrackV1::kProlongation : track prolongation failed +// - AliTRDtrackV1::kPropagation : track prolongation failed +// - AliTRDtrackV1::kAdjustSector : failed during sector crossing +// - AliTRDtrackV1::kSnp : too large bending +// - AliTRDtrackV1::kTrackletInit : fail to initialize tracklet +// - AliTRDtrackV1::kUpdate : fail to attach clusters or fit the tracklet +// - AliTRDtrackV1::kUnknown : anything which is not covered before +// 2. By default the status of the track before first TRD update is saved. +// +// Debug level 2 +// +// Author +// Alexandru Bercuci +// + + Int_t n = 0; + Double_t driftLength = .5*AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick(); AliTRDtrackingChamber *chamber = 0x0; AliTRDseedV1 tracklet, *ptrTracklet = 0x0; @@ -618,149 +692,260 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) tracklets[ip] = t.GetTracklet(ip); t.UnsetTracklet(ip); } + Bool_t kStoreIn = kTRUE, kPropagateIn = kTRUE; // Loop through the TRD layers - for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) { - // BUILD TRACKLET IF NOT ALREADY BUILT - Double_t x = 0., x0, y, z, alpha; - ptrTracklet = tracklets[ilayer]; - if(!ptrTracklet){ - ptrTracklet = new(&tracklet) AliTRDseedV1(ilayer); - ptrTracklet->SetReconstructor(fReconstructor); - alpha = t.GetAlpha(); - Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsector)); + TGeoHMatrix *matrix = 0x0; + Double_t x, y, z; + for (Int_t ily=0, sm=-1, stk=-1, det=-1; ily < AliTRDgeometry::kNlayer; ily++) { + // rough estimate of the entry point + if (!t.GetProlongation(fR[ily], y, z)){ + n=-1; + t.SetStatus(AliTRDtrackV1::kProlongation); + break; + } - if(!fTrSec[sector].GetNChambers()) continue; - - if((x = fTrSec[sector].GetX(ilayer)) < 1.) continue; - - // Propagate closer to the current layer - x0 = x - 1.5*clength; - if (x0 > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x0-fgkMaxStep, fgkMaxStep)) return -1/*nClustersExpected*/; - if (!AdjustSector(&t)) return -1/*nClustersExpected*/; - if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/; - - if (!t.GetProlongation(x, y, z)) return -1/*nClustersExpected*/; - Int_t stack = fGeom->GetStack(z, ilayer); - Int_t nCandidates = stack >= 0 ? 1 : 2; - z -= stack >= 0 ? 0. : 4.; - - for(int icham=0; ichamGetStack(z, ilayer)) < 0) continue; - - if(!(chamber = fTrSec[sector].GetChamber(stack, ilayer))) continue; - - if(chamber->GetNClusters() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()) continue; - - x = chamber->GetX(); - - AliTRDpadPlane *pp = fGeom->GetPadPlane(ilayer, stack); - tracklet.SetTilt(TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle())); - tracklet.SetPadLength(pp->GetLengthIPad()); - tracklet.SetDetector(chamber->GetDetector()); - tracklet.SetX0(x); - tracklet.UpDate(&t); -// if(!tracklet.Init(&t)){ -// t.SetStopped(kTRUE); -// return nClustersExpected; -// } - if(!tracklet.AttachClusters(chamber, kTRUE)) continue; - //if(!tracklet.AttachClustersIter(chamber, 1000.)) continue; - //tracklet.Init(&t); - - if(tracklet.GetN() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()) continue; + // find sector / stack / detector + sm = t.GetSector(); + // TODO cross check with y value ! + stk = fGeom->GetStack(z, ily); + det = stk>=0 ? AliTRDgeometry::GetDetector(ily, stk, sm) : -1; + matrix = det>=0 ? fGeom->GetClusterMatrix(det) : 0x0; + + // check if supermodule/chamber is installed + if( !fGeom->GetSMstatus(sm) || + stk<0. || + fGeom->IsHole(ily, stk, sm) || + !matrix ){ + // propagate to the default radial position + if(fR[ily] > (fgkMaxStep + t.GetX()) && !PropagateToX(t, fR[ily], fgkMaxStep)){ + n=-1; + t.SetStatus(AliTRDtrackV1::kPropagation); + break; + } + if(!AdjustSector(&t)){ + n=-1; + t.SetStatus(AliTRDtrackV1::kAdjustSector); + break; + } + if(TMath::Abs(t.GetSnp()) > fgkMaxSnp){ + n=-1; + t.SetStatus(AliTRDtrackV1::kSnp); + break; + } + t.SetStatus(AliTRDtrackV1::kGeometry, ily); + continue; + } + + // retrieve rotation matrix for the current chamber + Double_t loc[] = {AliTRDgeometry::AnodePos()- driftLength, 0., 0.}; + Double_t glb[] = {0., 0., 0.}; + matrix->LocalToMaster(loc, glb); + + // Propagate to the radial distance of the current layer + x = glb[0] - fgkMaxStep; + if(x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x, fgkMaxStep)){ + n=-1; + t.SetStatus(AliTRDtrackV1::kPropagation); + break; + } + if(!AdjustSector(&t)){ + n=-1; + t.SetStatus(AliTRDtrackV1::kAdjustSector); + break; + } + if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) { + n=-1; + t.SetStatus(AliTRDtrackV1::kSnp); + break; + } + Bool_t RECALCULATE = kFALSE; + if(sm != t.GetSector()){ + sm = t.GetSector(); + RECALCULATE = kTRUE; + } + if(stk != fGeom->GetStack(z, ily)){ + stk = fGeom->GetStack(z, ily); + RECALCULATE = kTRUE; + } + if(RECALCULATE){ + det = AliTRDgeometry::GetDetector(ily, stk, sm); + if(!(matrix = fGeom->GetClusterMatrix(det))){ + t.SetStatus(AliTRDtrackV1::kGeometry, ily); + continue; + } + matrix->LocalToMaster(loc, glb); + x = glb[0] - fgkMaxStep; + } + + // check if track is well inside fiducial volume + if (!t.GetProlongation(x+fgkMaxStep, y, z)) { + n=-1; + t.SetStatus(AliTRDtrackV1::kProlongation); + break; + } + if(fGeom->IsOnBoundary(det, y, z, .5)){ + t.SetStatus(AliTRDtrackV1::kBoundary, ily); + continue; + } + // mark track as entering the FIDUCIAL volume of TRD + if(kStoreIn){ + t.SetTrackLow(); + kStoreIn = kFALSE; + } + + ptrTracklet = tracklets[ily]; + if(!ptrTracklet){ // BUILD TRACKLET + // check data in supermodule + if(!fTrSec[sm].GetNChambers()){ + t.SetStatus(AliTRDtrackV1::kNoClusters, ily); + continue; + } + if(fTrSec[sm].GetX(ily) < 1.){ + t.SetStatus(AliTRDtrackV1::kNoClusters, ily); + continue; + } + // check data in chamber + if(!(chamber = fTrSec[sm].GetChamber(stk, ily))){ + t.SetStatus(AliTRDtrackV1::kNoClusters, ily); + continue; + } + if(chamber->GetNClusters() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()){ + t.SetStatus(AliTRDtrackV1::kNoClusters, ily); + continue; + } + // build tracklet + ptrTracklet = new(&tracklet) AliTRDseedV1(det); + ptrTracklet->SetReconstructor(fReconstructor); + ptrTracklet->SetKink(t.IsKink()); + ptrTracklet->SetPadPlane(fGeom->GetPadPlane(ily, stk)); + ptrTracklet->SetX0(glb[0]+driftLength); + if(!tracklet.Init(&t)){ + n=-1; + t.SetStatus(AliTRDtrackV1::kTrackletInit); break; } - //ptrTracklet->UseClusters(); - }// else ptrTracklet->Init(&t); - if(!ptrTracklet->IsOK()){ - if(x < 1.) continue; //temporary - if(!PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) return -1/*nClustersExpected*/; - if(!AdjustSector(&t)) return -1/*nClustersExpected*/; - if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/; + if(!tracklet.AttachClusters(chamber, kTRUE)){ + t.SetStatus(AliTRDtrackV1::kNoAttach, ily); + continue; + } + if(tracklet.GetN() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()){ + t.SetStatus(AliTRDtrackV1::kNoClustersTracklet, ily); + continue; + } + ptrTracklet->UpdateUsed(); + } + // propagate track to the radial position of the tracklet + ptrTracklet->UseClusters(); // TODO ? do we need this here ? + // fit tracklet no tilt correction + if(!ptrTracklet->Fit(kFALSE)){ + t.SetStatus(AliTRDtrackV1::kNoFit, ily); continue; + } + x = ptrTracklet->GetX(); //GetX0(); + if(x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x, fgkMaxStep)) { + n=-1; + t.SetStatus(AliTRDtrackV1::kPropagation); + break; } - - // Propagate closer to the current chamber if neccessary - x -= clength; - if (x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) return -1/*nClustersExpected*/; - if (!AdjustSector(&t)) return -1/*nClustersExpected*/; - if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/; - - // load tracklet to the tracker and the track - ptrTracklet->Fit(kFALSE); // no tilt correction + if(!AdjustSector(&t)) { + n=-1; + t.SetStatus(AliTRDtrackV1::kAdjustSector); + break; + } + if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) { + n=-1; + t.SetStatus(AliTRDtrackV1::kSnp); + break; + } + if(kPropagateIn){ + t.SetTrackLow(); + kPropagateIn = kFALSE; + } + Double_t cov[3]; ptrTracklet->GetCovAt(x, cov); + Double_t p[2] = { ptrTracklet->GetY(), ptrTracklet->GetZ()}; + Double_t chi2 = ((AliExternalTrackParam)t).GetPredictedChi2(p, cov); + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 2){ + Double_t ytrack = ptrTracklet->GetYref(0); + Double_t ztrack = ptrTracklet->GetZref(0); + Double_t ytracklet = ptrTracklet->GetYfit(0); + Double_t ztracklet = ptrTracklet->GetZfit(0); + Double_t phitrack = ptrTracklet->GetYref(1); + Double_t phitracklet = ptrTracklet->GetYfit(1); + Double_t thetatrack = ptrTracklet->GetZref(1); + Double_t thetatracklet = ptrTracklet->GetZfit(1); + + TTreeSRedirector &mystreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); + mystreamer << "FollowBackProlongation1" + << "il=" << ily + << "x=" << x + << "ytrack=" << ytrack + << "ztrack=" << ztrack + << "ytracklet=" << ytracklet + << "ztracklet=" << ztracklet + << "phitrack=" << phitrack + << "thetatrack=" << thetatrack + << "phitracklet=" << phitracklet + << "thetatracklet=" << thetatracklet + << "chi2=" << chi2 + << "\n"; + } + // update Kalman with the TRD measurement + if(chi2>1e+10){ // TODO + t.SetStatus(AliTRDtrackV1::kChi2, ily); + continue; + } + if(!t.Update(p, cov, chi2)) { + n=-1; + t.SetStatus(AliTRDtrackV1::kUpdate); + break; + } + // fill residuals ?! + AliTracker::FillResiduals(&t, p, cov, ptrTracklet->GetVolumeId()); + + + // load tracklet to the tracker + ptrTracklet->Update(&t); ptrTracklet = SetTracklet(ptrTracklet); t.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1); - - - // Calculate the mean material budget along the path inside the chamber - //Calculate global entry and exit positions of the track in chamber (only track prolongation) - Double_t xyz0[3]; // entry point - t.GetXYZ(xyz0); - alpha = t.GetAlpha(); - x = ptrTracklet->GetXref(); //GetX0(); - if (!t.GetProlongation(x, y, z)) return -1/*nClustersExpected*/; - Double_t xyz1[3]; // exit point - xyz1[0] = x * TMath::Cos(alpha) - y * TMath::Sin(alpha); - xyz1[1] = +x * TMath::Sin(alpha) + y * TMath::Cos(alpha); - xyz1[2] = z; - Double_t param[7]; - if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) return -1; - // The mean propagation parameters - Double_t xrho = param[0]*param[4]; // density*length - Double_t xx0 = param[1]; // radiation length - - // Propagate and update track - if (!t.PropagateTo(x, xx0, xrho)) return -1/*nClustersExpected*/; - if (!AdjustSector(&t)) return -1/*nClustersExpected*/; - Double_t maxChi2 = t.GetPredictedChi2(ptrTracklet); - if (!t.Update(ptrTracklet, maxChi2)) return -1/*nClustersExpected*/; - ptrTracklet->UpDate(&t); + n += ptrTracklet->GetN(); - if (maxChi2<1e+10) { - nClustersExpected += ptrTracklet->GetN(); - //t.SetTracklet(&tracklet, index); - } // Reset material budget if 2 consecutive gold - if(ilayer>0 && t.GetTracklet(ilayer-1) && ptrTracklet->GetN() + t.GetTracklet(ilayer-1)->GetN() > 20) t.SetBudget(2, 0.); +// if(ilayer>0 && t.GetTracklet(ilayer-1) && ptrTracklet->GetN() + t.GetTracklet(ilayer-1)->GetN() > 20) t.SetBudget(2, 0.); // Make backup of the track until is gold // TO DO update quality check of the track. // consider comparison with fTimeBinsRange Float_t ratio0 = ptrTracklet->GetN() / Float_t(fgNTimeBins); //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1); - //printf("tracklet.GetChi2() %f [< 18.0]\n", tracklet.GetChi2()); - //printf("ratio0 %f [> 0.8]\n", ratio0); - //printf("ratio1 %f [> 0.6]\n", ratio1); - //printf("ratio0+ratio1 %f [> 1.5]\n", ratio0+ratio1); - //printf("t.GetNCross() %d [== 0]\n", t.GetNCross()); - //printf("TMath::Abs(t.GetSnp()) %f [< 0.85]\n", TMath::Abs(t.GetSnp())); - //printf("t.GetNumberOfClusters() %d [> 20]\n", t.GetNumberOfClusters()); - if (//(tracklet.GetChi2() < 18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update + if( (chi2 < 18.0) && (ratio0 > 0.8) && //(ratio1 > 0.6) && //(ratio0+ratio1 > 1.5) && (t.GetNCross() == 0) && (TMath::Abs(t.GetSnp()) < 0.85) && - (t.GetNumberOfClusters() > 20)) t.MakeBackupTrack(); - + (t.GetNumberOfClusters() > 20)){ + t.MakeBackupTrack(); + } } // end layers loop + //printf("clusters[%d] chi2[%f] x[%f] status[%d ", n, t.GetChi2(), t.GetX(), t.GetStatusTRD()); + //for(int i=0; i<6; i++) printf("%d ", t.GetStatusTRD(i)); printf("]\n"); if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); - //AliTRDtrackV1 *debugTrack = new AliTRDtrackV1(t); - //debugTrack->SetOwner(); + AliTRDtrackV1 track(t); + track.SetOwner(); cstreamer << "FollowBackProlongation" - << "EventNumber=" << eventNumber - << "ncl=" << nClustersExpected - //<< "track.=" << debugTrack + << "EventNumber=" << eventNumber + << "ncl=" << n + << "track.=" << &track << "\n"; } - return nClustersExpected; + return n; } //_________________________________________________________________________ @@ -787,7 +972,7 @@ Float_t AliTRDtrackerV1::FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_ } for(Int_t il = 0; il < maxLayers; il++){ if(!tracklets[ppl[il]].IsOK()) continue; - fitter->AddPoint(tracklets[ppl[il]].GetX0(), tracklets[ppl[il]].GetYfitR(0), tracklets[ppl[il]].GetZProb(),1,10); + fitter->AddPoint(tracklets[ppl[il]].GetX0(), tracklets[ppl[il]].GetYfit(0), tracklets[ppl[il]].GetZfit(0),1,10); } fitter->Update(); // Set the reference position of the fit and calculate the chi2 values @@ -820,8 +1005,9 @@ void AliTRDtrackerV1::FitRieman(AliTRDcluster **seedcl, Double_t chi2[2]) // AliRieman *fitter = AliTRDtrackerV1::GetRiemanFitter(); fitter->Reset(); - for(Int_t i = 0; i < 4; i++) - fitter->AddPoint(seedcl[i]->GetX(), seedcl[i]->GetY(), seedcl[i]->GetZ(), 1, 10); + for(Int_t i = 0; i < 4; i++){ + fitter->AddPoint(seedcl[i]->GetX(), seedcl[i]->GetY(), seedcl[i]->GetZ(), 1., 10.); + } fitter->Update(); @@ -874,9 +1060,10 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub Int_t nPoints = 0; for(Int_t ilr = 0; ilr < AliTRDgeometry::kNlayer; ilr++){ if(!tracklets[ilr].IsOK()) continue; - for(Int_t itb = 0; itb < fgNTimeBins; itb++){ + for(Int_t itb = 0; itb < AliTRDseedV1::kNclusters; itb++){ if(!tracklets[ilr].IsUsable(itb)) continue; cl = tracklets[ilr].GetClusters(itb); + if(!cl->IsInChamber()) continue; x = cl->GetX(); y = cl->GetY(); z = cl->GetZ(); @@ -886,7 +1073,7 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub uvt[0] = 2. * x * t; uvt[1] = 2. * x * t * tilt ; w = 2. * (y + tilt * (z - zVertex)) * t; - error = 2. * 0.2 * t; + error = 2. * TMath::Sqrt(cl->GetSigmaY2()+tilt*tilt*cl->GetSigmaZ2()) * t; fitter->AddPoint(uvt, w, error); nPoints++; } @@ -900,7 +1087,7 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub Float_t chi2track = fitter->GetChisquare()/Double_t(nPoints); for(Int_t ip = 0; ip < AliTRDtrackerV1::kNPlanes; ip++) - tracklets[ip].SetCC(curvature); + tracklets[ip].SetC(curvature); /* if(fReconstructor->GetStreamLevel() >= 5){ //Linear Model on z-direction @@ -962,19 +1149,22 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro AliTRDcluster *cl = 0x0; Double_t xref = CalculateReferenceX(tracklets); - Double_t x, y, z, t, tilt, dx, w, we; - Double_t uvt[4]; + Double_t x, y, z, t, tilt, dx, w, we, erry, errz; + Double_t uvt[4], sumPolY[5], sumPolZ[3]; + memset(sumPolY, 0, sizeof(Double_t) * 5); + memset(sumPolZ, 0, sizeof(Double_t) * 3); 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 < fgNTimeBins; itb++){ + tilt = tracklets[ipl].GetTilt(); + for(Int_t itb = 0; itb < AliTRDseedV1::kNclusters; itb++){ if(!(cl = tracklets[ipl].GetClusters(itb))) continue; + if(!cl->IsInChamber()) 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); @@ -985,9 +1175,21 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro w = 2. * (y + tilt*z) * t; // error definition changes for the different calls we = 2. * t; - we *= sigError ? tracklets[ipl].GetSigmaY() : 0.2; + we *= sigError ? TMath::Sqrt(cl->GetSigmaY2()+tilt*tilt*cl->GetSigmaZ2()) : 0.2; fitter->AddPoint(uvt, w, we); zfitter.AddPoint(&x, z, static_cast(TMath::Sqrt(cl->GetSigmaZ2()))); + // adding points for covariance matrix estimation + erry = 1./(TMath::Sqrt(cl->GetSigmaY2()) + 0.1); // 0.1 is a systematic error (due to misalignment and miscalibration) + erry *= erry; + errz = 1./cl->GetSigmaZ2(); + for(Int_t ipol = 0; ipol < 5; ipol++){ + sumPolY[ipol] += erry; + erry *= x; + if(ipol < 3){ + sumPolZ[ipol] += errz; + errz *= x; + } + } nPoints++; } } @@ -1004,7 +1206,7 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) { if(!tracklets[iLayer].IsOK()) continue; zref = offset + slope * (tracklets[iLayer].GetX0() - xref); - if (TMath::Abs(tracklets[iLayer].GetZProb() - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0) + if (TMath::Abs(tracklets[iLayer].GetZfit(0) - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0) acceptablez = kFALSE; } if (!acceptablez) { @@ -1029,16 +1231,45 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro Double_t chi2track = fitter->GetChisquare()/Double_t(nPoints); + // Prepare error calculation + TMatrixD covarPolY(3,3); + covarPolY(0,0) = sumPolY[0]; covarPolY(1,1) = sumPolY[2]; covarPolY(2,2) = sumPolY[4]; + covarPolY(0,1) = covarPolY(1,0) = sumPolY[1]; + covarPolY(0,2) = covarPolY(2,0) = sumPolY[2]; + covarPolY(2,1) = covarPolY(1,2) = sumPolY[3]; + covarPolY.Invert(); + TMatrixD covarPolZ(2,2); + covarPolZ(0,0) = sumPolZ[0]; covarPolZ(1,1) = sumPolZ[2]; + covarPolZ(1,0) = covarPolZ(0,1) = sumPolZ[1]; + covarPolZ.Invert(); + // Update the tracklets - Double_t dy, dz; + Double_t x1, dy, dz; + Double_t cov[15]; + memset(cov, 0, sizeof(Double_t) * 15); for(Int_t iLayer = 0; iLayer < AliTRDtrackerV1::kNPlanes; iLayer++) { x = tracklets[iLayer].GetX0(); + x1 = x - xref; y = 0; z = 0; dy = 0; dz = 0; - + memset(cov, 0, sizeof(Double_t) * 3); + TMatrixD transform(3,3); + transform(0,0) = 1; + transform(0,1) = x; + transform(0,2) = x*x; + transform(1,1) = 1; + transform(1,2) = x; + transform(2,2) = 1; + TMatrixD covariance(transform, TMatrixD::kMult, covarPolY); + covariance *= transform.T(); + TMatrixD transformZ(2,2); + transformZ(0,0) = transformZ(1,1) = 1; + transformZ(0,1) = x; + TMatrixD covarZ(transformZ, TMatrixD::kMult, covarPolZ); + covarZ *= transformZ.T(); // y: R^2 = (x - x0)^2 + (y - y0)^2 // => y = y0 +/- Sqrt(R^2 - (x - x0)^2) // R = Sqrt() = 1/Curvature @@ -1050,6 +1281,9 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro res = TMath::Sqrt(res); y = (1.0 - res) / a; } + cov[0] = covariance(0,0); + cov[2] = covarZ(0,0); + cov[1] = 0.; // dy: R^2 = (x - x0)^2 + (y - y0)^2 // => y = +/- Sqrt(R^2 - (x - x0)^2) + y0 @@ -1059,9 +1293,9 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro Double_t x0 = -b / a; if (-c * a + b * b + 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 (a < 0) yderiv *= -1.0; - dy = yderiv; + Double_t yderiv = (x - x0) / TMath::Sqrt(1.0/(curvature * curvature) - (x - x0) * (x - x0)); + if (a < 0) yderiv *= -1.0; + dy = yderiv; } } z = offset + slope * (x - xref); @@ -1071,6 +1305,7 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro tracklets[iLayer].SetZref(0, z); tracklets[iLayer].SetZref(1, dz); tracklets[iLayer].SetC(curvature); + tracklets[iLayer].SetCovRef(cov); tracklets[iLayer].SetChi2(chi2track); } @@ -1161,35 +1396,56 @@ Double_t AliTRDtrackerV1::FitLine(const AliTRDtrackV1 *track, AliTRDseedV1 *trac //_________________________________________________________________________ Double_t AliTRDtrackerV1::FitRiemanTilt(const AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t sigError, Int_t np, AliTrackPoint *points) { - // - // 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) - // +// +// Performs a Riemann fit taking tilting pad correction into account +// +// Paramters: - Array of tracklets (connected to the track candidate) +// - Flag selecting the error definition +// Output: - Chi2 values of the track (in Parameter list) +// +// The equations which has to be solved simultaneously are: +// BEGIN_LATEX +// R^{2} = (x-x_{0})^{2} + (y^{*}-y_{0})^{2} +// y^{*} = y - tg(h)(z - z_{t}) +// z_{t} = z_{0}+dzdx*(x-x_{r}) +// END_LATEX +// with (x, y, z) the coordinate of the cluster, (x_0, y_0, z_0) the coordinate of the center of the Riemann circle, +// R its radius, x_r a constant refrence radial position in the middle of the TRD stack and dzdx the slope of the +// track in the x-z plane. Using the following transformations +// BEGIN_LATEX +// t = 1 / (x^{2} + y^{2}) +// u = 2 * x * t +// v = 2 * tan(h) * t +// w = 2 * tan(h) * (x - x_{r}) * t +// END_LATEX +// One gets the following linear equation +// BEGIN_LATEX +// a + b * u + c * t + d * v + e * w = 2 * (y + tg(h) * z) * t +// END_LATEX +// where the coefficients have the following meaning +// BEGIN_LATEX +// a = -1/y_{0} +// b = x_{0}/y_{0} +// c = (R^{2} -x_{0}^{2} - y_{0}^{2})/y_{0} +// d = z_{0} +// e = dz/dx +// END_LATEX +// The error calculation for the free term is thus +// BEGIN_LATEX +// #sigma = 2 * #sqrt{#sigma^{2}_{y} + (tilt corr ...) + tg^{2}(h) * #sigma^{2}_{z}} * t +// END_LATEX +// +// From this simple model one can compute chi^2 estimates and a rough approximation of pt from the curvature according +// to the formula: +// BEGIN_LATEX +// C = 1/R = a/(1 + b^{2} + c*a) +// END_LATEX +// +// Authors +// M.Ivanov +// A.Bercuci +// M.Fasel + TLinearFitter *fitter = GetTiltedRiemanFitter(); fitter->StoreData(kTRUE); fitter->ClearPoints(); @@ -1213,7 +1469,7 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(const AliTRDtrackV1 *track, AliTRDseedV1 // Containers for Least-square fitter for(Int_t ipl = 0; ipl < kNPlanes; ipl++){ if(!tracklets[ipl].IsOK()) continue; - for(Int_t itb = 0; itb < fgNTimeBins; itb++){ + for(Int_t itb = 0; itb < AliTRDseedV1::kNclusters; itb++){ if(!(cl = tracklets[ipl].GetClusters(itb))) continue; if (!tracklets[ipl].IsUsable(itb)) continue; x = cl->GetX(); @@ -1230,7 +1486,7 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(const AliTRDtrackV1 *track, AliTRDseedV1 w = 2. * (y + tilt*z) * t; // error definition changes for the different calls we = 2. * t; - we *= sigError ? tracklets[ipl].GetSigmaY() : 0.2; + we *= sigError ? TMath::Sqrt(cl->GetSigmaY2()) : 0.2; fitter->AddPoint(uvt, w, we); zfitter.AddPoint(&x, z, static_cast(TMath::Sqrt(cl->GetSigmaZ2()))); nPoints++; @@ -1249,7 +1505,7 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(const AliTRDtrackV1 *track, AliTRDseedV1 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].GetZProb() - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0) + if (TMath::Abs(tracklets[iLayer].GetZfit(0) - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0) accept = kFALSE; } if (!accept) { @@ -1304,7 +1560,7 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(const AliTRDtrackV1 *track, AliTRDseedV1 Float_t xyz[3]; for(int ip=0; ip R ? 100. : y0 - (y0>0.?1.:-1.)*TMath::Sqrt(R*R-(xyz[0]-x0)*(xyz[0]-x0)); + 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); } @@ -1409,8 +1665,10 @@ Double_t AliTRDtrackerV1::FitKalman(AliTRDtrackV1 *track, AliTRDseedV1 *tracklet if (!AdjustSector(track)) break; //Update track - Double_t chi2 = track->GetPredictedChi2(ptrTracklet); - if(chi2<1e+10) track->Update(ptrTracklet, chi2); + Double_t cov[3]; ptrTracklet->GetCovAt(x, cov); + Double_t p[2] = { ptrTracklet->GetY(), ptrTracklet->GetZ()}; + Double_t chi2 = ((AliExternalTrackParam*)track)->GetPredictedChi2(p, cov); + if(chi2<1e+10) track->Update(p, cov, chi2); if(!up) continue; //Reset material budget if 2 consecutive gold @@ -1448,7 +1706,7 @@ Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset 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].GetMeanz() - z); + chi2Z += TMath::Abs(tracklets[iLayer].GetZfit(0) - z); nLayers++; } chi2Z /= TMath::Max((nLayers - 3.0),1.0); @@ -1491,9 +1749,7 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m 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 - } + if(t.GetProlongation(x,y,z)<0) 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()); @@ -1559,7 +1815,6 @@ Int_t AliTRDtrackerV1::ReadClusters(TClonesArray* &array, TTree *clusterTree) co Int_t nCluster = clusterArray->GetEntriesFast(); for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue; - c->SetInChamber(); new((*fClusters)[ncl++]) AliTRDcluster(*c); delete (clusterArray->RemoveAt(iCluster)); } @@ -1657,12 +1912,19 @@ Int_t AliTRDtrackerV1::BuildTrackingContainers() //____________________________________________________________________ void AliTRDtrackerV1::UnloadClusters() { - // - // Clears the arrays of clusters and tracks. Resets sectors and timebins - // +// +// Clears the arrays of clusters and tracks. Resets sectors and timebins +// If option "force" is also set the containers are also deleted. This is useful +// in case of HLT - if(fTracks) fTracks->Delete(); - if(fTracklets) fTracklets->Delete(); + if(fTracks){ + fTracks->Delete(); + if(HasRemoveContainers()){delete fTracks; fTracks = 0x0;} + } + if(fTracklets){ + fTracklets->Delete(); + if(HasRemoveContainers()){delete fTracklets; fTracklets = 0x0;} + } if(fClusters){ if(IsClustersOwner()) fClusters->Delete(); @@ -1679,23 +1941,23 @@ void AliTRDtrackerV1::UnloadClusters() AliTRDtrackerDebug::SetEventNumber(AliTRDtrackerDebug::GetEventNumber() + 1); } -//____________________________________________________________________ -void AliTRDtrackerV1::UseClusters(const AliKalmanTrack *t, Int_t) const -{ - const AliTRDtrackV1 *track = dynamic_cast(t); - if(!track) return; - - AliTRDseedV1 *tracklet = 0x0; - for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){ - if(!(tracklet = track->GetTracklet(ily))) continue; - AliTRDcluster *c = 0x0; - for(Int_t ic=AliTRDseed::knTimebins; ic--;){ - if(!(c=tracklet->GetClusters(ic))) continue; - c->Use(); - } - } -} - +// //____________________________________________________________________ +// void AliTRDtrackerV1::UseClusters(const AliKalmanTrack *t, Int_t) const +// { +// const AliTRDtrackV1 *track = dynamic_cast(t); +// if(!track) return; +// +// AliTRDseedV1 *tracklet = 0x0; +// for(Int_t ily=AliTRDgeometry::kNlayer; ily--;){ +// if(!(tracklet = track->GetTracklet(ily))) continue; +// AliTRDcluster *c = 0x0; +// for(Int_t ic=AliTRDseed::kNclusters; ic--;){ +// if(!(c=tracklet->GetClusters(ic))) continue; +// c->Use(); +// } +// } +// } +// //_____________________________________________________________________________ Bool_t AliTRDtrackerV1::AdjustSector(AliTRDtrackV1 *track) @@ -1876,6 +2138,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon const AliTRDCalDet *cal = AliTRDcalibDB::Instance()->GetT0Det(); AliTRDtrackingChamber *chamber = 0x0; + AliTRDtrackingChamber **ci = 0x0; AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized Int_t pars[4]; // MakeSeeds parameters @@ -1883,9 +2146,16 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon //Double_t shift = .5 * alpha; Int_t configs[kNConfigs]; + // Purge used clusters from the containers + ci = &stack[0]; + for(Int_t ic = kNPlanes; ic--; ci++){ + if(!(*ci)) continue; + (*ci)->Update(); + } + // Build initial seeding configurations Double_t quality = BuildSeedingConfigs(stack, configs); - if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 10){ AliInfo(Form("Plane config %d %d %d Quality %f" , configs[0], configs[1], configs[2], quality)); } @@ -1898,10 +2168,10 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon fSieveSeeding = 0; // Get stack index - Int_t ic = 0; AliTRDtrackingChamber **cIter = &stack[0]; - while(icGetStack((*cIter)->GetDetector()); + Int_t ic = 0; ci = &stack[0]; + while(icGetStack((*ci)->GetDetector()); do{ // Loop over seeding configurations @@ -1911,6 +2181,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon pars[1] = ntracks; pars[2] = istack; ntracks = MakeSeeds(stack, &sseed[6*ntracks], pars); + //AliInfo(Form("Number of Tracks after iteration step %d: %d\n", iconf, ntracks)); if(ntracks == kMaxTracksStack) break; } if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding)); @@ -1953,196 +2224,133 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) { Int_t jseed = kNPlanes*trackIndex+jLayer; if(!sseed[jseed].IsOK()) continue; - if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.15) findable++; - + if (TMath::Abs(sseed[jseed].GetYref(0) / sseed[jseed].GetX0()) < 0.158) findable++; + // TODO here we get a sig fault which should never happen ! sseed[jseed].UpdateUsed(); ncl += sseed[jseed].GetN2(); nused += sseed[jseed].GetNUsed(); nlayers++; } - // Filter duplicated tracks - if (nused > 30){ - //printf("Skip %d nused %d\n", trackIndex, nused); - fakeTrack[trackIndex] = kTRUE; - continue; - } - if (Float_t(nused)/ncl >= .25){ - //printf("Skip %d nused/ncl >= .25\n", trackIndex); - fakeTrack[trackIndex] = kTRUE; - continue; - } - - // Classify tracks - Bool_t skip = kFALSE; - switch(jSieve){ - case 0: - if(nlayers < 6) {skip = kTRUE; break;} - if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;} - break; - - case 1: - if(nlayers < findable){skip = kTRUE; break;} - if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;} - break; - - case 2: - if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;} - if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;} - break; - - case 3: - if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;} - break; - - case 4: - if (nlayers == 3){skip = kTRUE; break;} - //if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;} - break; - } - if(skip){ - candidates++; - //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused); - continue; - } - signedTrack[trackIndex] = kTRUE; - - - // Sign clusters - AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1; - for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) { - Int_t jseed = kNPlanes*trackIndex+jLayer; - if(!sseed[jseed].IsOK()) continue; - if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian - sseed[jseed].UseClusters(); - if(!cl){ - ic = 0; - while(!(cl = sseed[jseed].GetClusters(ic))) ic++; - clusterIndex = sseed[jseed].GetIndexes(ic); - } - } - if(!cl) continue; - - - // Build track parameters - AliTRDseedV1 *lseed =&sseed[trackIndex*6]; -/* Int_t idx = 0; - while(idx<3 && !lseed->IsOK()) { - idx++; - lseed++; - }*/ - Double_t x = lseed->GetX0();// - 3.5; - trackParams[0] = x; //NEW AB - trackParams[1] = lseed->GetYref(0); // lseed->GetYat(x); - trackParams[2] = lseed->GetZref(0); // lseed->GetZat(x); - trackParams[3] = TMath::Sin(TMath::ATan(lseed->GetYref(1))); - trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1)); - trackParams[5] = lseed->GetC(); - Int_t ich = 0; while(!(chamber = stack[ich])) ich++; - trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *alpha+shift; // Supermodule*/ - - if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ - AliInfo(Form("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1])); - - Int_t nclusters = 0; - AliTRDseedV1 *dseed[6]; - - // Build track label - what happens if measured data ??? - Int_t labels[1000]; - Int_t outlab[1000]; - Int_t nlab = 0; - - Int_t labelsall[1000]; - Int_t nlabelsall = 0; - Int_t naccepted = 0; - - for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) { - Int_t jseed = kNPlanes*trackIndex+iLayer; - dseed[iLayer] = new AliTRDseedV1(sseed[jseed]); - dseed[iLayer]->SetOwner(); - nclusters += sseed[jseed].GetN2(); - if(!sseed[jseed].IsOK()) continue; - for(int ilab=0; ilab<2; ilab++){ - if(sseed[jseed].GetLabels(ilab) < 0) continue; - labels[nlab] = sseed[jseed].GetLabels(ilab); - nlab++; - } - - // Cooking label - for (Int_t itime = 0; itime < fgNTimeBins; itime++) { - if(!sseed[jseed].IsUsable(itime)) continue; - naccepted++; - Int_t tindex = 0, ilab = 0; - while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){ - labelsall[nlabelsall++] = tindex; - ilab++; + // Filter duplicated tracks + if (nused > 30){ + //printf("Skip %d nused %d\n", trackIndex, nused); + fakeTrack[trackIndex] = kTRUE; + continue; + } + if (Float_t(nused)/ncl >= .25){ + //printf("Skip %d nused/ncl >= .25\n", trackIndex); + fakeTrack[trackIndex] = kTRUE; + continue; } - } - } - Freq(nlab,labels,outlab,kFALSE); - Int_t label = outlab[0]; - Int_t frequency = outlab[1]; - Freq(nlabelsall,labelsall,outlab,kFALSE); - Int_t label1 = outlab[0]; - Int_t label2 = outlab[2]; - Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted); - //Int_t eventNrInFile = esd->GetEventNumberInFile(); - //AliInfo(Form("Number of clusters %d.", nclusters)); - Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); - Int_t trackNumber = AliTRDtrackerDebug::GetTrackNumber(); - Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); - TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); - cstreamer << "Clusters2TracksStack" - << "EventNumber=" << eventNumber - << "TrackNumber=" << trackNumber - << "CandidateNumber=" << candidateNumber - << "Iter=" << fSieveSeeding - << "Like=" << fTrackQuality[trackIndex] - << "S0.=" << dseed[0] - << "S1.=" << dseed[1] - << "S2.=" << dseed[2] - << "S3.=" << dseed[3] - << "S4.=" << dseed[4] - << "S5.=" << dseed[5] - << "p0=" << trackParams[0] - << "p1=" << trackParams[1] - << "p2=" << trackParams[2] - << "p3=" << trackParams[3] - << "p4=" << trackParams[4] - << "p5=" << trackParams[5] - << "p6=" << trackParams[6] - << "Label=" << label - << "Label1=" << label1 - << "Label2=" << label2 - << "FakeRatio=" << fakeratio - << "Freq=" << frequency - << "Ncl=" << ncl - << "NLayers=" << nlayers - << "Findable=" << findable - << "NUsed=" << nused - << "\n"; - } + // Classify tracks + Bool_t skip = kFALSE; + switch(jSieve){ + case 0: + if(nlayers < 6) {skip = kTRUE; break;} + if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;} + break; + + case 1: + if(nlayers < findable){skip = kTRUE; break;} + if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;} + break; + + case 2: + if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;} + if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;} + break; + + case 3: + if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;} + break; + + case 4: + if (nlayers == 3){skip = kTRUE; break;} + //if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;} + break; + } + if(skip){ + candidates++; + //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused); + continue; + } + signedTrack[trackIndex] = kTRUE; + + // Build track parameters + AliTRDseedV1 *lseed =&sseed[trackIndex*6]; + /* Int_t idx = 0; + while(idx<3 && !lseed->IsOK()) { + idx++; + lseed++; + }*/ + Double_t x = lseed->GetX0();// - 3.5; + trackParams[0] = x; //NEW AB + trackParams[1] = lseed->GetYref(0); // lseed->GetYat(x); + trackParams[2] = lseed->GetZref(0); // lseed->GetZat(x); + trackParams[3] = TMath::Sin(TMath::ATan(lseed->GetYref(1))); + trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1)); + trackParams[5] = lseed->GetC(); + Int_t ich = 0; while(!(chamber = stack[ich])) ich++; + trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *alpha+shift; // Supermodule*/ + + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ + //AliInfo(Form("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1])); + + AliTRDseedV1 *dseed[6]; + for(Int_t iseed = AliTRDgeometry::kNlayer; iseed--;) dseed[iseed] = new AliTRDseedV1(lseed[iseed]); + + //Int_t eventNrInFile = esd->GetEventNumberInFile(); + Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); + Int_t trackNumber = AliTRDtrackerDebug::GetTrackNumber(); + Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); + TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); + cstreamer << "Clusters2TracksStack" + << "EventNumber=" << eventNumber + << "TrackNumber=" << trackNumber + << "CandidateNumber=" << candidateNumber + << "Iter=" << fSieveSeeding + << "Like=" << fTrackQuality[trackIndex] + << "S0.=" << dseed[0] + << "S1.=" << dseed[1] + << "S2.=" << dseed[2] + << "S3.=" << dseed[3] + << "S4.=" << dseed[4] + << "S5.=" << dseed[5] + << "p0=" << trackParams[0] + << "p1=" << trackParams[1] + << "p2=" << trackParams[2] + << "p3=" << trackParams[3] + << "p4=" << trackParams[4] + << "p5=" << trackParams[5] + << "p6=" << trackParams[6] + << "Ncl=" << ncl + << "NLayers=" << nlayers + << "Findable=" << findable + << "NUsed=" << nused + << "\n"; + } + + AliTRDtrackV1 *track = MakeTrack(&sseed[trackIndex*kNPlanes], trackParams); + if(!track){ + AliWarning("Fail to build a TRD Track."); + continue; + } - AliTRDtrackV1 *track = MakeTrack(&sseed[trackIndex*kNPlanes], trackParams); - if(!track){ - AliWarning("Fail to build a TRD Track."); - continue; - } - - //AliInfo("End of MakeTrack()"); - AliESDtrack *esdTrack = new ((*esdTrackList)[ntracks0++]) AliESDtrack(); - esdTrack->UpdateTrackParams(track, AliESDtrack::kTRDout); - esdTrack->SetLabel(track->GetLabel()); - track->UpdateESDtrack(esdTrack); - // write ESD-friends if neccessary - if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0){ - AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(*track); - calibTrack->SetOwner(); - esdTrack->AddCalibObject(calibTrack); - } - ntracks1++; - AliTRDtrackerDebug::SetTrackNumber(AliTRDtrackerDebug::GetTrackNumber() + 1); + //AliInfo("End of MakeTrack()"); + AliESDtrack *esdTrack = new ((*esdTrackList)[ntracks0++]) AliESDtrack(); + esdTrack->UpdateTrackParams(track, AliESDtrack::kTRDout); + esdTrack->SetLabel(track->GetLabel()); + track->UpdateESDtrack(esdTrack); + // write ESD-friends if neccessary + if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0){ + AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(*track); + calibTrack->SetOwner(); + esdTrack->AddCalibObject(calibTrack); + } + ntracks1++; + AliTRDtrackerDebug::SetTrackNumber(AliTRDtrackerDebug::GetTrackNumber() + 1); } jSieve++; @@ -2164,7 +2372,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon chamber->Build(fGeom, cal);//Indices(fSieveSeeding); } - if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 10){ AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality)); } } while(fSieveSeeding<10); // end stack clusters sieve @@ -2203,15 +2411,15 @@ Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int // The overall chamber quality is given by the product of this 2 contributions. // - Double_t chamberQ[kNPlanes]; + Double_t chamberQ[kNPlanes];memset(chamberQ, 0, kNPlanes*sizeof(Double_t)); AliTRDtrackingChamber *chamber = 0x0; for(int iplane=0; iplaneGetQuality() : 0.; } - Double_t tconfig[kNConfigs]; - Int_t planes[4]; + Double_t tconfig[kNConfigs];memset(tconfig, 0, kNConfigs*sizeof(Double_t)); + Int_t planes[] = {0, 0, 0, 0}; for(int iconf=0; iconf seeding chambers configuration - // ipar[1] -> stack index - // ipar[2] -> number of track candidates found so far - // - // Output : - // Number of tracks candidates found. - // - // Detailed description - // - // The following steps are performed: - // 1. Select seeding layers from seeding chambers - // 2. Select seeding clusters from the seeding AliTRDpropagationLayerStack. - // The clusters are taken from layer 3, layer 0, layer 1 and layer 2, in - // this order. The parameters controling the range of accepted clusters in - // layer 0, 1, and 2 are defined in AliTRDchamberTimeBin::BuildCond(). - // 3. Helix fit of the cluster set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**)) - // 4. Initialize seeding tracklets in the seeding chambers. - // 5. Filter 0. - // Chi2 in the Y direction less than threshold ... (1./(3. - sLayer)) - // Chi2 in the Z direction less than threshold ... (1./(3. - sLayer)) - // 6. Attach clusters to seeding tracklets and find linear approximation of - // the tracklet (see AliTRDseedV1::AttachClustersIter()). The number of used - // clusters used by current seeds should not exceed ... (25). - // 7. Filter 1. - // All 4 seeding tracklets should be correctly constructed (see - // AliTRDseedV1::AttachClustersIter()) - // 8. Helix fit of the seeding tracklets - // 9. Filter 2. - // Likelihood calculation of the fit. (See AliTRDtrackerV1::CookLikelihood() for details) - // 10. Extrapolation of the helix fit to the other 2 chambers: - // a) Initialization of extrapolation tracklet with fit parameters - // b) Helix fit of tracklets - // c) Attach clusters and linear interpolation to extrapolated tracklets - // d) Helix fit of tracklets - // 11. Improve seeding tracklets quality by reassigning clusters. - // See AliTRDtrackerV1::ImproveSeedQuality() for details. - // 12. Helix fit of all 6 seeding tracklets and chi2 calculation - // 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details. - // 14. Cooking labels for tracklets. Should be done only for MC - // 15. Register seeds. - // +// +// Seed tracklets and build candidate TRD tracks. The procedure is used during barrel tracking to account for tracks which are +// either missed by TPC prolongation or conversions inside the TRD volume. +// For stand alone tracking the procedure is used to estimate all tracks measured by TRD. +// +// Parameters : +// layers : Array of stack propagation layers containing clusters +// sseed : Array of empty tracklet seeds. On exit they are filled. +// ipar : Control parameters: +// ipar[0] -> seeding chambers configuration +// ipar[1] -> stack index +// ipar[2] -> number of track candidates found so far +// +// Output : +// Number of tracks candidates found. +// +// The following steps are performed: +// 1. Build seeding layers by collapsing all time bins from each of the four seeding chambers along the +// radial coordinate. See AliTRDtrackingChamber::GetSeedingLayer() for details. The chambers selection for seeding +// is described in AliTRDtrackerV1::Clusters2TracksStack(). +// 2. Using the seeding clusters from the seeding layer (step 1) build combinatorics using the following algorithm: +// - for each seeding cluster in the lower seeding layer find +// - all seeding clusters in the upper seeding layer inside a road defined by a given phi angle. The angle +// is calculated on the minimum pt of tracks from vertex accesible to the stand alone tracker. +// - for each pair of two extreme seeding clusters select middle upper cluster using roads defined externally by the +// reco params +// - select last seeding cluster as the nearest to the linear approximation of the track described by the first three +// seeding clusters. +// The implementation of road calculation and cluster selection can be found in the functions AliTRDchamberTimeBin::BuildCond() +// and AliTRDchamberTimeBin::GetClusters(). +// 3. Helix fit of the seeding clusters set. (see AliTRDtrackerFitter::FitRieman(AliTRDcluster**)). No tilt correction is +// performed at this level +// 4. Initialize seeding tracklets in the seeding chambers. +// 5. *Filter 0* Chi2 cut on the Y and Z directions. The threshold is set externally by the reco params. +// 6. Attach (true) clusters to seeding tracklets (see AliTRDseedV1::AttachClusters()) and fit tracklet (see +// AliTRDseedV1::Fit()). The number of used clusters used by current seeds should not exceed ... (25). +// 7. *Filter 1* Check if all 4 seeding tracklets are correctly constructed. +// 8. Helix fit of the clusters from the seeding tracklets with tilt correction. Refit tracklets using the new +// approximation of the track. +// 9. *Filter 2* Calculate likelihood of the track. (See AliTRDtrackerV1::CookLikelihood()). The following quantities are +// checked against the Riemann fit: +// - position resolution in y +// - angular resolution in the bending plane +// - likelihood of the number of clusters attached to the tracklet +// 10. Extrapolation of the helix fit to the other 2 chambers *non seeding* chambers: +// - Initialization of extrapolation tracklets with the fit parameters +// - Attach clusters to extrapolated tracklets +// - Helix fit of tracklets +// 11. Improve seeding tracklets quality by reassigning clusters based on the last parameters of the track +// See AliTRDtrackerV1::ImproveSeedQuality() for details. +// 12. Helix fit of all 6 seeding tracklets and chi2 calculation +// 13. Hyperplane fit and track quality calculation. See AliTRDtrackerFitter::FitHyperplane() for details. +// 14. Cooking labels for tracklets. Should be done only for MC +// 15. Register seeds. +// +// Authors: +// Marian Ivanov +// Alexandru Bercuci +// Markus Fasel AliTRDtrackingChamber *chamber = 0x0; AliTRDcluster *c[kNSeedPlanes] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters @@ -2288,11 +2507,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss // chi2[1] = tracklet chi2 on the R direction Double_t chi2[4]; - // Default positions for the anode wire in all 6 Layers in case of a stack with missing clusters - // Positions taken using cosmic data taken with SM3 after rebuild - Double_t x_def[kNPlanes] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; - - // this should be data member of AliTRDtrack + // this should be data member of AliTRDtrack TODO Double_t seedQuality[kMaxTracksStack]; // unpack control parameters @@ -2300,20 +2515,39 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss Int_t ntracks = ipar[1]; Int_t istack = ipar[2]; Int_t planes[kNSeedPlanes]; GetSeedingConfig(config, planes); - Int_t planesExt[kNPlanes-kNSeedPlanes]; GetExtrapolationConfig(config, planesExt); + Int_t planesExt[kNPlanes-kNSeedPlanes]; GetExtrapolationConfig(config, planesExt); // Init chambers geometry Double_t hL[kNPlanes]; // Tilting angle Float_t padlength[kNPlanes]; // pad lenghts + Float_t padwidth[kNPlanes]; // pad widths AliTRDpadPlane *pp = 0x0; for(int iplane=0; iplaneGetPadPlane(iplane, istack); hL[iplane] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle()); padlength[iplane] = pp->GetLengthIPad(); + padwidth[iplane] = pp->GetWidthIPad(); + } + + // Init anode wire position for chambers + Double_t x0[kNPlanes], // anode wire position + driftLength = .5*AliTRDgeometry::AmThick() - AliTRDgeometry::DrThick(); // drift length + TGeoHMatrix *matrix = 0x0; + Double_t loc[] = {AliTRDgeometry::AnodePos(), 0., 0.}; + Double_t glb[] = {0., 0., 0.}; + AliTRDtrackingChamber **cIter = &stack[0]; + for(int iLayer=0; iLayerGetClusterMatrix((*cIter)->GetDetector()))){ + continue; + x0[iLayer] = fgkX0[iLayer]; + } + matrix->LocalToMaster(loc, glb); + x0[iLayer] = glb[0]; } - - if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ + + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 2){ AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks)); } @@ -2325,7 +2559,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss if(!chamber->GetSeedingLayer(fSeedTB[isl], fGeom, fReconstructor)) continue; nlayers++; } - if(nlayers < 4) return ntracks; + if(nlayers < kNSeedPlanes) return ntracks; // Start finding seeds @@ -2341,9 +2575,9 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss c[0] = (*fSeedTB[0])[index[jcl++]]; if(!c[0]) continue; Double_t dx = c[3]->GetX() - c[0]->GetX(); - Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx; - Double_t phi = (c[3]->GetY() - c[0]->GetY())/dx; - fSeedTB[1]->BuildCond(c[0], cond1, 1, theta, phi); + Double_t dzdx = (c[3]->GetZ() - c[0]->GetZ())/dx; + Double_t dydx = (c[3]->GetY() - c[0]->GetY())/dx; + fSeedTB[1]->BuildCond(c[0], cond1, 1, dzdx, dydx); fSeedTB[1]->GetClusters(cond1, jndex, mcl); //printf("Found c[0] candidates 1 %d\n", mcl); @@ -2351,27 +2585,30 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss while(kclBuildCond(c[1], cond2, 2, theta, phi); + fSeedTB[2]->BuildCond(c[1], cond2, 2, dzdx, dydx); c[2] = fSeedTB[2]->GetNearestCluster(cond2); //printf("Found c[1] candidate 2 %p\n", c[2]); if(!c[2]) continue; - // AliInfo("Seeding clusters found. Building seeds ..."); - // for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %6.3f, y = %6.3f, z = %6.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ()); + //AliInfo("Seeding clusters found. Building seeds ..."); + //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %6.3f, y = %6.3f, z = %6.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ()); for (Int_t il = 0; il < kNPlanes; il++) cseed[il].Reset(); FitRieman(c, chi2); AliTRDseedV1 *tseed = &cseed[0]; - AliTRDtrackingChamber **cIter = &stack[0]; + cIter = &stack[0]; for(int iLayer=0; iLayerSetDetector((*cIter) ? (*cIter)->GetDetector() : -1); + Int_t det = (*cIter) ? (*cIter)->GetDetector() : -1; + tseed->SetDetector(det); tseed->SetTilt(hL[iLayer]); tseed->SetPadLength(padlength[iLayer]); + tseed->SetPadWidth(padwidth[iLayer]); tseed->SetReconstructor(fReconstructor); - tseed->SetX0((*cIter) ? (*cIter)->GetX() : x_def[iLayer]); + tseed->SetX0(det<0 ? fR[iLayer]+driftLength : x0[iLayer]); tseed->Init(GetRiemanFitter()); + tseed->SetStandAlone(kTRUE); } Bool_t isFake = kFALSE; @@ -2419,25 +2656,32 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss <<"\n"; } if(chi2[0] > fReconstructor->GetRecoParam() ->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){ -// //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0])); + //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0])); AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1); continue; } if(chi2[1] > fReconstructor->GetRecoParam() ->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){ -// //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1])); + //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1])); AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1); continue; } //AliInfo("Passed chi2 filter."); // try attaching clusters to tracklets - Int_t nUsedCl = 0; - Int_t mlayers = 0; + Int_t mlayers = 0; + AliTRDcluster *cl = NULL; for(int iLayer=0; iLayer 25) break; + Int_t nNotInChamber = 0; + if(!cseed[jLayer].AttachClusters(stack[jLayer], kTRUE)) continue; + cseed[jLayer].Fit(); + cseed[jLayer].UpdateUsed(); + cseed[jLayer].ResetClusterIter(); + while((cl = cseed[jLayer].NextCluster())){ + if(!cl->IsInChamber()) nNotInChamber++; + } + //printf("clusters[%d], used[%d], not in chamber[%d]\n", cseed[jLayer].GetN(), cseed[jLayer].GetNUsed(), nNotInChamber); + if(cseed[jLayer].GetN() - (cseed[jLayer].GetNUsed() + nNotInChamber) < 5) continue; // checking for Cluster which are not in chamber is a much stronger restriction on real data mlayers++; } @@ -2453,7 +2697,8 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss for(int iLayer=0; iLayerGetRecoParam() ->GetTrackLikelihood()){ @@ -2479,7 +2729,6 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss fSeedLayer[ntracks] = config;/*sLayer;*/ // attach clusters to the extrapolation seeds - Int_t nusedf = 0; // debug value for(int iLayer=0; iLayerGetRecoParam()->HasImproveTracklets() && ImproveSeedQuality(stack, cseed) < 4){ + if(fReconstructor->HasImproveTracklets() && ImproveSeedQuality(stack, cseed) < 4){ AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1); continue; } @@ -2527,8 +2777,8 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss // do the final track fitting (Once with vertex constraint and once without vertex constraint) Double_t chi2Vals[3]; - chi2Vals[0] = FitTiltedRieman(&cseed[0], kFALSE); - if(fReconstructor->GetRecoParam()->IsVertexConstrained()) + chi2Vals[0] = FitTiltedRieman(&cseed[0], kTRUE); + if(fReconstructor->HasVertexConstrained()) chi2Vals[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ()); // Do Vertex Constrained fit if desired else chi2Vals[1] = 1.; @@ -2537,32 +2787,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss //chi2Vals[2] = GetChi2ZTest(&cseed[0]); fTrackQuality[ntracks] = CalculateTrackLikelihood(&cseed[0], &chi2Vals[0]); //AliInfo("Hyperplane fit done\n"); - - // finalize tracklets - Int_t labels[12]; - Int_t outlab[24]; - Int_t nlab = 0; - for (Int_t iLayer = 0; iLayer < 6; iLayer++) { - if (!cseed[iLayer].IsOK()) continue; - - if (cseed[iLayer].GetLabels(0) >= 0) { - labels[nlab] = cseed[iLayer].GetLabels(0); - nlab++; - } - - if (cseed[iLayer].GetLabels(1) >= 0) { - labels[nlab] = cseed[iLayer].GetLabels(1); - nlab++; - } - } - Freq(nlab,labels,outlab,kFALSE); - Int_t label = outlab[0]; - Int_t frequency = outlab[1]; - for (Int_t iLayer = 0; iLayer < 6; iLayer++) { - cseed[iLayer].SetFreq(frequency); - cseed[iLayer].SetChi2Z(chi2[1]); - } - + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){ TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); @@ -2580,8 +2805,6 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss << "Chi2TC=" << chi2Vals[1] << "Nlayers=" << mlayers << "NClusters=" << ncls - << "NUsedS=" << nUsedCl - << "NUsed=" << nusedf << "Like=" << like << "S0.=" << &cseed[0] << "S1.=" << &cseed[1] @@ -2589,8 +2812,6 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss << "S3.=" << &cseed[3] << "S4.=" << &cseed[4] << "S5.=" << &cseed[5] - << "Label=" << label - << "Freq=" << frequency << "FitterT.=" << fitterT << "FitterTC.=" << fitterTC << "\n"; @@ -2613,42 +2834,58 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss //_____________________________________________________________________________ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params) { - // - // Build a TRD track out of tracklet candidates - // - // Parameters : - // seeds : array of tracklets - // params : track parameters (see MakeSeeds() function body for a detailed description) - // - // Output : - // The TRD track. - // - // Detailed description - // - // To be discussed with Marian !! - // - +// +// Build a TRD track out of tracklet candidates +// +// Parameters : +// seeds : array of tracklets +// params : array of track parameters as they are estimated by stand alone tracker. 7 elements. +// [0] - radial position of the track at reference point +// [1] - y position of the fit at [0] +// [2] - z position of the fit at [0] +// [3] - snp of the first tracklet +// [4] - tgl of the first tracklet +// [5] - curvature of the Riemann fit - 1/pt +// [6] - sector rotation angle +// +// Output : +// The TRD track. +// +// Initialize the TRD track based on the parameters of the fit and a parametric covariance matrix +// (diagonal with constant variance terms TODO - correct parameterization) +// +// In case of HLT just register the tracklets in the tracker and return values of the Riemann fit. For the +// offline case perform a full Kalman filter on the already found tracklets (see AliTRDtrackerV1::FollowBackProlongation() +// for details). Do also MC label calculation and PID if propagation successfully. + Double_t alpha = AliTRDgeometry::GetAlpha(); Double_t shift = AliTRDgeometry::GetAlpha()/2.0; Double_t c[15]; - c[ 0] = 0.2; - c[ 1] = 0.0; c[ 2] = 2.0; - c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02; - c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1; - c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01; + c[ 0] = 0.2; // s^2_y + c[ 1] = 0.0; c[ 2] = 2.0; // s^2_z + c[ 3] = 0.0; c[ 4] = 0.0; c[ 5] = 0.02; // s^2_snp + c[ 6] = 0.0; c[ 7] = 0.0; c[ 8] = 0.0; c[ 9] = 0.1; // s^2_tgl + c[10] = 0.0; c[11] = 0.0; c[12] = 0.0; c[13] = 0.0; c[14] = params[5]*params[5]*0.01; // s^2_1/pt AliTRDtrackV1 track(seeds, ¶ms[1], c, params[0], params[6]*alpha+shift); track.PropagateTo(params[0]-5.0); + AliTRDseedV1 *ptrTracklet = 0x0; + + // skip Kalman filter for HLT if(fReconstructor->IsHLT()){ - AliTRDseedV1 *ptrTracklet = 0x0; - for(Int_t ip=0; ipIsOK()) continue; + if(TMath::Abs(ptrTracklet->GetYref(1) - ptrTracklet->GetYfit(1)) >= .2) continue; // check this condition with Marian + ptrTracklet = SetTracklet(ptrTracklet); + ptrTracklet->UseClusters(); track.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1); } AliTRDtrackV1 *ptrTrack = SetTrack(&track); + ptrTrack->CookPID(); ptrTrack->SetReconstructor(fReconstructor); return ptrTrack; } @@ -2707,7 +2944,7 @@ Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDs // layers : Array of propagation layers for a stack/supermodule // cseed : Array of 6 seeding tracklets which has to be improved // - // Output : + // Output : // cssed : Improved seeds // // Detailed description @@ -2749,7 +2986,8 @@ Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDs for (Int_t jLayer = 5; jLayer > 1; jLayer--) { Int_t bLayer = sortindexes[jLayer]; if(!(chamber = stack[bLayer])) continue; - bseed[bLayer].AttachClustersIter(chamber, squality[bLayer], kTRUE); + bseed[bLayer].AttachClusters(chamber, kTRUE); + bseed[bLayer].Fit(kTRUE); if(bseed[bLayer].IsOK()) nLayers++; } @@ -2773,7 +3011,6 @@ Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDs << "\n"; } } // Loop: iter - // we are sure that at least 2 tracklets are OK ! return nLayers+2; } @@ -2798,20 +3035,20 @@ Double_t AliTRDtrackerV1::CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Doub // debug level 2 // - Double_t sumdaf = 0, nLayers = 0; + Double_t chi2phi = 0, nLayers = 0; for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) { if(!tracklets[iLayer].IsOK()) continue; - sumdaf += TMath::Abs((tracklets[iLayer].GetYfit(1) - tracklets[iLayer].GetYref(1))/ tracklets[iLayer].GetSigmaY2()); + chi2phi += tracklets[iLayer].GetChi2Phi(); nLayers++; } - sumdaf /= Float_t (nLayers - 2.0); + chi2phi /= Float_t (nLayers - 2.0); Double_t likeChi2Z = TMath::Exp(-chi2[2] * 0.14); // Chi2Z - Double_t likeChi2TC = (fReconstructor->GetRecoParam() ->IsVertexConstrained()) ? + Double_t likeChi2TC = (fReconstructor->HasVertexConstrained()) ? TMath::Exp(-chi2[1] * 0.677) : 1; // Constrained Tilted Riemann - Double_t likeChi2TR = TMath::Exp(-chi2[0] * 0.78); // Non-constrained Tilted Riemann - Double_t likeAF = TMath::Exp(-sumdaf * 3.23); - Double_t trackLikelihood = likeChi2Z * likeChi2TR * likeAF; + Double_t likeChi2TR = TMath::Exp(-chi2[0] * 0.0078); // Non-constrained Tilted Riemann + Double_t likeChi2Phi= TMath::Exp(-chi2phi * 3.23);//3.23 + Double_t trackLikelihood = likeChi2Z * likeChi2TR * likeChi2Phi; if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){ Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); @@ -2823,11 +3060,11 @@ Double_t AliTRDtrackerV1::CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Doub << "LikeChi2Z=" << likeChi2Z << "LikeChi2TR=" << likeChi2TR << "LikeChi2TC=" << likeChi2TC - << "LikeAF=" << likeAF + << "LikeChi2Phi=" << likeChi2Phi << "TrackLikelihood=" << trackLikelihood << "\n"; } - + return trackLikelihood; } @@ -2869,7 +3106,7 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4]) for(UChar_t ilayer = 0; ilayer < 4; ilayer++){ Int_t jlayer = planes[ilayer]; nclusters += cseed[jlayer].GetN2(); - sumda += TMath::Abs(cseed[jlayer].GetYfitR(1) - cseed[jlayer].GetYref(1)); + sumda += TMath::Abs(cseed[jlayer].GetYfit(1) - cseed[jlayer].GetYref(1)); } nclusters *= .25; @@ -2880,7 +3117,6 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4]) Double_t likeN = TMath::Exp(-(fRecoPars->GetNMeanClusters() - nclusters) / fRecoPars->GetNSigmaClusters()); Double_t like = likea * likechi2y * likechi2z * likeN; - // AliInfo(Form("sumda(%f) chi2[0](%f) chi2[1](%f) likea(%f) likechi2y(%f) likechi2z(%f) nclusters(%d) likeN(%f)", sumda, chi2[0], chi2[1], likea, likechi2y, likechi2z, nclusters, likeN)); if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){ Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); @@ -3160,159 +3396,189 @@ AliKalmanTrack* AliTRDtrackerV1::GetTrack(Int_t idx) const return idx >= 0 && idx < ntrk ? (AliKalmanTrack*)fTracks->UncheckedAt(idx) : 0x0; } -//____________________________________________________________________ -Float_t AliTRDtrackerV1::CalculateReferenceX(AliTRDseedV1 *tracklets){ - // - // 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(tracklets[il].IsOK() && tracklets[il -1].IsOK()){ - Float_t xdiff = tracklets[il].GetX0() - tracklets[il -1].GetX0(); - meanDistance += xdiff; - nDistances++; - } - if(tracklets[il].IsOK()) startIndex = il; - } - if(tracklets[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(tracklets[il].IsOK()){ - xpos[iok] = tracklets[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 tracklets[startIndex].GetX0() + (2.5 - startIndex) * meanDistance - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick()); -} - -//_____________________________________________________________________________ -Int_t AliTRDtrackerV1::Freq(Int_t n, const Int_t *inlist - , Int_t *outlist, Bool_t down) -{ - // - // Sort eleements according occurancy - // The size of output array has is 2*n - // - if (n <= 0) { - return 0; - } - Int_t *sindexS = new Int_t[n]; // Temporary array for sorting - Int_t *sindexF = new Int_t[2*n]; - for (Int_t i = 0; i < n; i++) { - sindexF[i] = 0; - } - - TMath::Sort(n,inlist,sindexS,down); - - Int_t last = inlist[sindexS[0]]; - Int_t val = last; - sindexF[0] = 1; - sindexF[0+n] = last; - Int_t countPos = 0; +// //_____________________________________________________________________________ +// Int_t AliTRDtrackerV1::Freq(Int_t n, const Int_t *inlist +// , Int_t *outlist, Bool_t down) +// { +// // +// // Sort eleements according occurancy +// // The size of output array has is 2*n +// // +// +// if (n <= 0) { +// return 0; +// } +// +// Int_t *sindexS = new Int_t[n]; // Temporary array for sorting +// Int_t *sindexF = new Int_t[2*n]; +// for (Int_t i = 0; i < n; i++) { +// sindexF[i] = 0; +// } +// +// TMath::Sort(n,inlist,sindexS,down); +// +// Int_t last = inlist[sindexS[0]]; +// Int_t val = last; +// sindexF[0] = 1; +// sindexF[0+n] = last; +// Int_t countPos = 0; +// +// // Find frequency +// for (Int_t i = 1; i < n; i++) { +// val = inlist[sindexS[i]]; +// if (last == val) { +// sindexF[countPos]++; +// } +// else { +// countPos++; +// sindexF[countPos+n] = val; +// sindexF[countPos]++; +// last = val; +// } +// } +// if (last == val) { +// countPos++; +// } +// +// // Sort according frequency +// TMath::Sort(countPos,sindexF,sindexS,kTRUE); +// +// for (Int_t i = 0; i < countPos; i++) { +// outlist[2*i ] = sindexF[sindexS[i]+n]; +// outlist[2*i+1] = sindexF[sindexS[i]]; +// } +// +// delete [] sindexS; +// delete [] sindexF; +// +// return countPos; +// +// } - // Find frequency - for (Int_t i = 1; i < n; i++) { - val = inlist[sindexS[i]]; - if (last == val) { - sindexF[countPos]++; - } - else { - countPos++; - sindexF[countPos+n] = val; - sindexF[countPos]++; - last = val; - } - } - if (last == val) { - countPos++; - } - // Sort according frequency - TMath::Sort(countPos,sindexF,sindexS,kTRUE); +//____________________________________________________________________ +void AliTRDtrackerV1::ResetSeedTB() +{ +// reset buffer for seeding time bin layers. If the time bin +// layers are not allocated this function allocates them - for (Int_t i = 0; i < countPos; i++) { - outlist[2*i ] = sindexF[sindexS[i]+n]; - outlist[2*i+1] = sindexF[sindexS[i]]; + for(Int_t isl=0; islClear(); } - - delete [] sindexS; - delete [] sindexF; - - return countPos; - } -//____________________________________________________________________ - //_____________________________________________________________________________ Float_t AliTRDtrackerV1::GetChi2Y(AliTRDseedV1 *tracklets) const { - // Chi2 definition on y-direction + // Calculates normalized chi2 in y-direction + // chi2 = Sum chi2 / n_tracklets - Float_t chi2 = 0; - for(Int_t ipl = 0; ipl < kNPlanes; ipl++){ + Double_t chi2 = 0.; Int_t n = 0; + for(Int_t ipl = kNPlanes; ipl--;){ if(!tracklets[ipl].IsOK()) continue; - Double_t distLayer = (tracklets[ipl].GetYfit(0) - tracklets[ipl].GetYref(0));// /tracklets[ipl].GetSigmaY(); - chi2 += distLayer * distLayer; - } - return chi2; -} - -//____________________________________________________________________ -void AliTRDtrackerV1::ResetSeedTB() -{ -// reset buffer for seeding time bin layers. If the time bin -// layers are not allocated this function allocates them - - for(Int_t isl=0; islClear(); + chi2 += tracklets[ipl].GetChi2Y(); + n++; } + return n ? chi2/n : 0.; } //_____________________________________________________________________________ Float_t AliTRDtrackerV1::GetChi2Z(AliTRDseedV1 *tracklets) const { // Calculates normalized chi2 in z-direction + // chi2 = Sum chi2 / n_tracklets - Float_t chi2 = 0; - // chi2 = Sum ((z - zmu)/sigma)^2 - // Sigma for the z direction is defined as half of the padlength - for(Int_t ipl = 0; ipl < kNPlanes; ipl++){ + Double_t chi2 = 0; Int_t n = 0; + for(Int_t ipl = kNPlanes; ipl--;){ if(!tracklets[ipl].IsOK()) continue; - Double_t distLayer = (tracklets[ipl].GetMeanz() - tracklets[ipl].GetZref(0)); // /(tracklets[ipl].GetPadLength()/2); - chi2 += distLayer * distLayer; + chi2 += tracklets[ipl].GetChi2Z(); + n++; + } + return n ? chi2/n : 0.; +} + +//____________________________________________________________________ +Float_t AliTRDtrackerV1::CalculateReferenceX(AliTRDseedV1 *tracklets){ + // + // 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) + // Only kept for compatibility with the old code + // + Int_t nDistances = 0; + Float_t meanDistance = 0.; + Int_t startIndex = 5; + for(Int_t il =5; il > 0; il--){ + if(tracklets[il].IsOK() && tracklets[il -1].IsOK()){ + Float_t xdiff = tracklets[il].GetX0() - tracklets[il -1].GetX0(); + meanDistance += xdiff; + nDistances++; + } + if(tracklets[il].IsOK()) startIndex = il; + } + if(tracklets[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(tracklets[il].IsOK()){ + xpos[iok] = tracklets[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 tracklets[startIndex].GetX0() + (2.5 - startIndex) * meanDistance - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick()); +} + +//_____________________________________________________________________________ +Double_t AliTRDtrackerV1::FitTiltedRiemanV1(AliTRDseedV1 *tracklets){ + // + // Track Fitter Function using the new class implementation of + // the Rieman fit + // + AliTRDtrackFitterRieman fitter; + fitter.SetRiemanFitter(GetTiltedRiemanFitter()); + fitter.Reset(); + for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++) fitter.SetTracklet(il, &tracklets[il]); + Double_t chi2 = fitter.Eval(); + // Update the tracklets + Double_t cov[15]; Double_t x0; + memset(cov, 0, sizeof(Double_t) * 15); + for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){ + x0 = tracklets[il].GetX0(); + tracklets[il].SetYref(0, fitter.GetYat(x0)); + tracklets[il].SetZref(0, fitter.GetZat(x0)); + tracklets[il].SetYref(1, fitter.GetDyDxAt(x0)); + tracklets[il].SetZref(1, fitter.GetDzDx()); + tracklets[il].SetC(fitter.GetCurvature()); + fitter.GetCovAt(x0, cov); + tracklets[il].SetCovRef(cov); + tracklets[il].SetChi2(chi2); } return chi2; } @@ -3339,7 +3605,8 @@ void AliTRDtrackerV1::AliTRDLeastSquare::AddPoint(Double_t *x, Double_t y, Doubl // // Adding Point to the fitter // - Double_t weight = 1/(sigmaY * sigmaY); + Double_t weight = 1/(sigmaY > 1e-9 ? sigmaY : 1e-9); + weight *= weight; Double_t &xpt = *x; // printf("Adding point x = %f, y = %f, sigma = %f\n", xpt, y, sigmaY); fSums[0] += weight; @@ -3355,7 +3622,9 @@ void AliTRDtrackerV1::AliTRDLeastSquare::RemovePoint(Double_t *x, Double_t y, Do // // Remove Point from the sample // - Double_t weight = 1/(sigmaY * sigmaY); + + Double_t weight = 1/(sigmaY > 1e-9 ? sigmaY : 1e-9); + weight *= weight; Double_t &xpt = *x; fSums[0] -= weight; fSums[1] -= weight * xpt; @@ -3384,9 +3653,9 @@ void AliTRDtrackerV1::AliTRDLeastSquare::Eval(){ // printf("fParams[0] = %f, fParams[1] = %f\n", fParams[0], fParams[1]); // Covariance matrix - fCovarianceMatrix[0] = fSums[4] - fSums[1] * fSums[1] / fSums[0]; - fCovarianceMatrix[1] = fSums[5] - fSums[2] * fSums[2] / fSums[0]; - fCovarianceMatrix[2] = fSums[3] - fSums[1] * fSums[2] / fSums[0]; + fCovarianceMatrix[0] = fSums[4] / fSums[0] - fSums[1] * fSums[1] / (fSums[0] * fSums[0]); + fCovarianceMatrix[1] = fSums[5] / fSums[0] - fSums[2] * fSums[2] / (fSums[0] * fSums[0]); + fCovarianceMatrix[2] = fSums[3] / fSums[0] - fSums[1] * fSums[2] / (fSums[0] * fSums[0]); } //_____________________________________________________________________________ @@ -3405,3 +3674,335 @@ void AliTRDtrackerV1::AliTRDLeastSquare::GetCovarianceMatrix(Double_t *storage) memcpy(storage, fCovarianceMatrix, sizeof(Double_t) * 3); } +//_____________________________________________________________________________ +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 *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(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()); +}