TLinearFitter* AliTRDtrackerV1::fgTiltedRiemanConstrained = 0x0;
//____________________________________________________________________
-AliTRDtrackerV1::AliTRDtrackerV1()
+AliTRDtrackerV1::AliTRDtrackerV1(AliTRDReconstructor *rec)
:AliTracker()
+ ,fReconstructor(0x0)
,fGeom(new AliTRDgeometry())
,fClusters(0x0)
,fTracklets(0x0)
//
// Default constructor.
//
- if (!AliTRDcalibDB::Instance()) {
+ AliTRDcalibDB *trd = 0x0;
+ if (!(trd = AliTRDcalibDB::Instance())) {
AliFatal("Could not get calibration object");
}
- fgNTimeBins = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
- for (Int_t isector = 0; isector < AliTRDgeometry::kNsect; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector);
+ if(!fgNTimeBins) fgNTimeBins = trd->GetNumberOfTimeBins();
+
+ for (Int_t isector = 0; isector < AliTRDgeometry::kNsector; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector);
- if(AliTRDReconstructor::StreamLevel() > 1){
- TDirectory *savedir = gDirectory;
- fgDebugStreamer = new TTreeSRedirector("TRD.TrackerDebug.root");
- savedir->cd();
- }
+ for(Int_t isl =0; isl<kNSeedPlanes; isl++) fSeedTB[isl] = 0x0;
+
+ // Initialize debug stream
+ if(rec) SetReconstructor(rec);
}
//____________________________________________________________________
if(fgRieman) delete fgRieman;
if(fgTiltedRieman) delete fgTiltedRieman;
if(fgTiltedRiemanConstrained) delete fgTiltedRiemanConstrained;
+ for(Int_t isl =0; isl<kNSeedPlanes; isl++) if(fSeedTB[isl]) delete fSeedTB[isl];
if(fTracks) {fTracks->Delete(); delete fTracks;}
if(fTracklets) {fTracklets->Delete(); delete fTracklets;}
- if(fClusters) {fClusters->Delete(); delete fClusters;}
+ if(fClusters) {
+ fClusters->Delete(); delete fClusters;
+ }
if(fGeom) delete fGeom;
}
// See AliTRDtrackerV1::Clusters2TracksSM() for details.
//
- if(!AliTRDReconstructor::RecoParam()){
+ if(!fReconstructor->GetRecoParam() ){
AliError("Reconstruction configuration not initialized. Call first AliTRDReconstructor::SetRecoParam().");
return 0;
}
//AliInfo("Start Track Finder ...");
Int_t ntracks = 0;
- for(int ism=0; ism<AliTRDgeometry::kNsect; ism++){
+ for(int ism=0; ism<AliTRDgeometry::kNsector; ism++){
// for(int ism=1; ism<2; ism++){
//AliInfo(Form("Processing supermodule %i ...", ism));
ntracks += Clusters2TracksSM(ism, esd);
{
//AliInfo(Form("Asking for tracklet %d", index));
- if(index<0 || index == 0xffff) return kFALSE;
- AliTRDseedV1 *tracklet = 0x0;
- if(!(tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index))) return kFALSE;
+ AliTRDseedV1 *tracklet = GetTracklet(index);
+ if (!tracklet) return kFALSE;
// get detector for this tracklet
AliTRDcluster *cl = 0x0;
- Int_t ic = 0; do; while(!(cl = tracklet->GetClusters(ic++)));
+ Int_t ic = 0; do {} while(!(cl = tracklet->GetClusters(ic++)));
Int_t idet = cl->GetDetector();
Double_t local[3];
// setting volume id
AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
- switch (fGeom->GetPlane(idet)) {
+ switch (fGeom->GetLayer(idet)) {
case 0:
iLayer = AliGeomManager::kTRD1;
break;
iLayer = AliGeomManager::kTRD6;
break;
};
- Int_t modId = fGeom->GetSector(idet) * fGeom->Ncham() + fGeom->GetChamber(idet);
+ Int_t modId = fGeom->GetSector(idet) * fGeom->Nstack() + fGeom->GetStack(idet);
UShort_t volid = AliGeomManager::LayerToVolUID(iLayer, modId);
p.SetVolumeID(volid);
//____________________________________________________________________
AliRieman* AliTRDtrackerV1::GetRiemanFitter()
{
- if(!fgRieman) fgRieman = new AliRieman(AliTRDtrackingChamber::kNTimeBins * AliTRDgeometry::kNplan);
+ if(!fgRieman) fgRieman = new AliRieman(AliTRDtrackingChamber::kNTimeBins * AliTRDgeometry::kNlayer);
return fgRieman;
}
Int_t found = 0; // number of tracks found
Float_t foundMin = 20.0;
+ Float_t *quality = 0x0;
+ Int_t *index = 0x0;
Int_t nSeed = event->GetNumberOfTracks();
- if(!nSeed){
- // run stand alone tracking
- if (AliTRDReconstructor::SeedingOn()) Clusters2Tracks(event);
- return 0;
- }
-
- Float_t *quality = new Float_t[nSeed];
- Int_t *index = new Int_t[nSeed];
- for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) {
- AliESDtrack *seed = event->GetTrack(iSeed);
- Double_t covariance[15];
- seed->GetExternalCovariance(covariance);
- quality[iSeed] = covariance[0] + covariance[2];
+ if(nSeed){
+ quality = new Float_t[nSeed];
+ index = new Int_t[nSeed];
+ for (Int_t iSeed = 0; iSeed < nSeed; 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);
}
- // Sort tracks according to covariance of local Y and Z
- TMath::Sort(nSeed,quality,index,kFALSE);
// Backpropagate all seeds
Int_t expectedClr;
// Do the back prolongation
new(&track) AliTRDtrackV1(*seed);
- //track->Print();
+ track.SetReconstructor(fReconstructor);
+
//Int_t lbl = seed->GetLabel();
//track.SetSeedLabel(lbl);
- seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup); // Make backup
+
+ // Make backup and mark entrance in the TRD
+ seed->UpdateTrackParams(&track, AliESDtrack::kTRDin);
+ seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup);
Float_t p4 = track.GetC();
- if((expectedClr = FollowBackProlongation(track))){
+ expectedClr = FollowBackProlongation(track);
+
+ if (expectedClr<0) continue; // Back prolongation failed
+
+ if(expectedClr){
+ found++;
// 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*/){
+ AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
+ calibTrack->SetOwner();
+ seed->AddCalibObject(calibTrack);
+ }
+ //update ESD track
if ((track.GetNumberOfClusters() > 15) && (track.GetNumberOfClusters() > 0.5*expectedClr)) {
seed->UpdateTrackParams(&track, AliESDtrack::kTRDout);
-
track.UpdateESDtrack(seed);
-
- // Add TRD track to ESDfriendTrack
- if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
- AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
- calibTrack->SetOwner();
- seed->AddCalibObject(calibTrack);
- }
}
}
//track.CookdEdxTimBin(seed->GetID());
track.CookLabel(1. - fgkLabelFraction);
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;
Double_t c2 = track.GetSnp() + track.GetC() * (xtof - track.GetX());
if (TMath::Abs(c2) >= 0.99) continue;
- PropagateToX(track, xTOF0, fgkMaxStep);
+ if (!PropagateToX(track, xTOF0, fgkMaxStep)) continue;
// Energy losses taken to the account - check one more time
c2 = track.GetSnp() + track.GetC() * (xtof - track.GetX());
if (track.PropagateTo(xtof)) {
seed->UpdateTrackParams(&track, AliESDtrack::kTRDout);
track.UpdateESDtrack(seed);
-
- // Add TRD track to ESDfriendTrack
-// if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
-// AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
-// calibTrack->SetOwner();
-// seed->AddCalibObject(calibTrack);
-// }
- found++;
}
} else {
if ((track.GetNumberOfClusters() > 15) && (track.GetNumberOfClusters() > 0.5*expectedClr)) {
seed->UpdateTrackParams(&track, AliESDtrack::kTRDout);
track.UpdateESDtrack(seed);
-
- // Add TRD track to ESDfriendTrack
-// if (AliTRDReconstructor::StreamLevel() > 0 /*&& quality TODO*/){
-// AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track);
-// calibTrack->SetOwner();
-// seed->AddCalibObject(calibTrack);
-// }
- found++;
}
}
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));
- delete [] index;
- delete [] quality;
+ // run stand alone tracking
+ if (fReconstructor->IsSeeding()) Clusters2Tracks(event);
return 0;
}
}
ULong_t status = seed->GetStatus();
+ // reject tracks which failed propagation in the TRD
if((status & AliESDtrack::kTRDout) == 0) continue;
- if((status & AliESDtrack::kTRDin) != 0) continue;
+
+ // reject tracks which are produced by the TRD stand alone track finder.
+ if((status & AliESDtrack::kTRDin) == 0) continue;
nseed++;
track.ResetCovariance(50.0);
return 0;
}
-
//____________________________________________________________________
Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
{
// Get material budget
Double_t param[7];
- AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
+ if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) break;
Double_t xrho= param[0]*param[4];
Double_t xx0 = param[1]; // Get mean propagation parameters
}
}
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
Int_t index;
for(int iplane=0; iplane<6; iplane++){
AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index);
AliTRDtrackingChamber *chamber = 0x0;
AliTRDseedV1 tracklet, *ptrTracklet = 0x0;
+ // in case of stand alone tracking we store all the pointers to the tracklets in a temporary array
+ AliTRDseedV1 *tracklets[kNPlanes];
+ memset(tracklets, 0, sizeof(AliTRDseedV1 *) * kNPlanes);
+ for(Int_t ip = 0; ip < kNPlanes; ip++){
+ tracklets[ip] = t.GetTracklet(ip);
+ t.UnsetTracklet(ip);
+ }
- // Loop through the TRD planes
- for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) {
+ // Loop through the TRD layers
+ for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) {
// BUILD TRACKLET IF NOT ALREADY BUILT
Double_t x = 0., y, z, alpha;
- ptrTracklet = t.GetTracklet(iplane);
+ ptrTracklet = tracklets[ilayer];
if(!ptrTracklet){
- ptrTracklet = new(&tracklet) AliTRDseedV1(iplane);
+ ptrTracklet = new(&tracklet) AliTRDseedV1(ilayer);
+ ptrTracklet->SetReconstructor(fReconstructor);
alpha = t.GetAlpha();
- Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsect));
+ Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsector));
if(!fTrSec[sector].GetNChambers()) continue;
- if((x = fTrSec[sector].GetX(iplane)) < 1.) continue;
+ if((x = fTrSec[sector].GetX(ilayer)) < 1.) continue;
- if (!t.GetProlongation(x, y, z)) break;
- Int_t stack = fGeom->GetChamber(z, iplane);
+ 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; icham<nCandidates; icham++, z+=8){
- if((stack = fGeom->GetChamber(z, iplane)) < 0) continue;
+ if((stack = fGeom->GetStack(z, ilayer)) < 0) continue;
- if(!(chamber = fTrSec[sector].GetChamber(stack, iplane))) continue;
+ if(!(chamber = fTrSec[sector].GetChamber(stack, ilayer))) continue;
- if(chamber->GetNClusters() < fgNTimeBins*AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
+ if(chamber->GetNClusters() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()) continue;
x = chamber->GetX();
- AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, stack);
- tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()));
+ AliTRDpadPlane *pp = fGeom->GetPadPlane(ilayer, stack);
+ tracklet.SetTilt(TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle()));
tracklet.SetPadLength(pp->GetLengthIPad());
- tracklet.SetPlane(iplane);
+ tracklet.SetDetector(chamber->GetDetector());
tracklet.SetX0(x);
if(!tracklet.Init(&t)){
t.SetStopped(kTRUE);
return nClustersExpected;
}
- if(!tracklet.AttachClustersIter(chamber, 1000.)) continue;
+ if(!tracklet.AttachClustersIter(chamber, 1000./*, kTRUE*/)) continue;
tracklet.Init(&t);
- if(tracklet.GetN() < fgNTimeBins * AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
+ if(tracklet.GetN() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()) continue;
break;
}
+ //ptrTracklet->UseClusters();
}
if(!ptrTracklet->IsOK()){
if(x < 1.) continue; //temporary
- if(!PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) break;
- if(!AdjustSector(&t)) break;
- if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
+ if(!PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) return -1/*nClustersExpected*/;
+ if(!AdjustSector(&t)) return -1/*nClustersExpected*/;
+ if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/;
continue;
}
// Propagate closer to the current chamber if neccessary
x -= clength;
- if (x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) break;
- if (!AdjustSector(&t)) break;
- if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break;
+ 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 = SetTracklet(ptrTracklet);
t.GetXYZ(xyz0);
alpha = t.GetAlpha();
x = ptrTracklet->GetX0();
- if (!t.GetProlongation(x, y, z)) break;
+ 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];
- AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
+ 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
- t.PropagateTo(x, xx0, xrho);
- if (!AdjustSector(&t)) break;
+ if (!t.PropagateTo(x, xx0, xrho)) return -1/*nClustersExpected*/;
+ if (!AdjustSector(&t)) return -1/*nClustersExpected*/;
Double_t maxChi2 = t.GetPredictedChi2(ptrTracklet);
- if (maxChi2<1e+10 && t.Update(ptrTracklet, maxChi2)){
+ if (!t.Update(ptrTracklet, maxChi2)) return -1/*nClustersExpected*/;
+ if (maxChi2<1e+10) {
nClustersExpected += ptrTracklet->GetN();
//t.SetTracklet(&tracklet, index);
}
// Reset material budget if 2 consecutive gold
- if(iplane>0 && t.GetTracklet(iplane-1) && ptrTracklet->GetN() + t.GetTracklet(iplane-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.
(TMath::Abs(t.GetSnp()) < 0.85) &&
(t.GetNumberOfClusters() > 20)) t.MakeBackupTrack();
- } // end planes loop
+ } // end layers loop
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
TTreeSRedirector &cstreamer = *fgDebugStreamer;
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
//AliTRDtrackV1 *debugTrack = new AliTRDtrackV1(t);
Float_t x, y, z, w, t, error, tilt;
Double_t uvt[2];
Int_t nPoints = 0;
- for(Int_t ipl = 0; ipl < AliTRDgeometry::kNplan; ipl++){
- if(!tracklets[ipl].IsOK()) continue;
+ for(Int_t ilr = 0; ilr < AliTRDgeometry::kNlayer; ilr++){
+ if(!tracklets[ilr].IsOK()) continue;
for(Int_t itb = 0; itb < fgNTimeBins; itb++){
- if(!tracklets[ipl].IsUsable(itb)) continue;
- cl = tracklets[ipl].GetClusters(itb);
+ if(!tracklets[ilr].IsUsable(itb)) continue;
+ cl = tracklets[ilr].GetClusters(itb);
x = cl->GetX();
y = cl->GetY();
z = cl->GetZ();
- tilt = tracklets[ipl].GetTilt();
+ tilt = tracklets[ilr].GetTilt();
// Transformation
t = 1./(x * x + y * y);
uvt[0] = 2. * x * t;
for(Int_t ip = 0; ip < AliTRDtrackerV1::kNPlanes; ip++)
tracklets[ip].SetCC(curvature);
- if(AliTRDReconstructor::StreamLevel() >= 5){
+/* if(fReconstructor->GetStreamLevel() >= 5){
//Linear Model on z-direction
Double_t xref = CalculateReferenceX(tracklets); // Relative to the middle of the stack
Double_t slope = fitter->GetParameter(2);
<< "Chi2Z=" << chi2Z
<< "zref=" << zref
<< "\n";
- }
+ }*/
return chi2track;
}
tracklets[iLayer].SetChi2(chi2track);
}
- if(AliTRDReconstructor::StreamLevel() >=5){
+/* if(fReconstructor->GetStreamLevel() >=5){
TTreeSRedirector &cstreamer = *fgDebugStreamer;
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
<< "xref=" << xref
<< "Chi2Z=" << chi2z
<< "\n";
- }
+ }*/
return chi2track;
}
+//____________________________________________________________________
+Double_t AliTRDtrackerV1::FitLine(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t err, Int_t np, AliTrackPoint *points)
+{
+ AliTRDLeastSquare yfitter, zfitter;
+ AliTRDcluster *cl = 0x0;
+
+ AliTRDseedV1 work[kNPlanes], *tracklet = 0x0;
+ if(!tracklets){
+ for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+ if(!(tracklet = track->GetTracklet(ipl))) continue;
+ if(!tracklet->IsOK()) continue;
+ new(&work[ipl]) AliTRDseedV1(*tracklet);
+ }
+ tracklets = &work[0];
+ }
+
+ Double_t xref = CalculateReferenceX(tracklets);
+ Double_t x, y, z, dx, ye, yr, tilt;
+ for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+ if(!tracklets[ipl].IsOK()) continue;
+ for(Int_t itb = 0; itb < fgNTimeBins; itb++){
+ if(!(cl = tracklets[ipl].GetClusters(itb))) continue;
+ if (!tracklets[ipl].IsUsable(itb)) continue;
+ x = cl->GetX();
+ z = cl->GetZ();
+ dx = x - xref;
+ zfitter.AddPoint(&dx, z, static_cast<Double_t>(TMath::Sqrt(cl->GetSigmaZ2())));
+ }
+ }
+ zfitter.Eval();
+ Double_t z0 = zfitter.GetFunctionParameter(0);
+ Double_t dzdx = zfitter.GetFunctionParameter(1);
+ for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
+ if(!tracklets[ipl].IsOK()) continue;
+ for(Int_t itb = 0; itb < fgNTimeBins; itb++){
+ if(!(cl = tracklets[ipl].GetClusters(itb))) continue;
+ if (!tracklets[ipl].IsUsable(itb)) continue;
+ x = cl->GetX();
+ y = cl->GetY();
+ z = cl->GetZ();
+ tilt = tracklets[ipl].GetTilt();
+ dx = x - xref;
+ yr = y + tilt*(z - z0 - dzdx*dx);
+ // error definition changes for the different calls
+ ye = tilt*TMath::Sqrt(cl->GetSigmaZ2());
+ ye += err ? tracklets[ipl].GetSigmaY() : 0.2;
+ yfitter.AddPoint(&dx, yr, ye);
+ }
+ }
+ yfitter.Eval();
+ Double_t y0 = yfitter.GetFunctionParameter(0);
+ Double_t dydx = yfitter.GetFunctionParameter(1);
+ Double_t chi2 = 0.;//yfitter.GetChisquare()/Double_t(nPoints);
+
+ //update track points array
+ if(np && points){
+ Float_t xyz[3];
+ for(int ip=0; ip<np; ip++){
+ points[ip].GetXYZ(xyz);
+ xyz[1] = y0 + dydx * (xyz[0] - xref);
+ xyz[2] = z0 + dzdx * (xyz[0] - xref);
+ points[ip].SetXYZ(xyz);
+ }
+ }
+ return chi2;
+}
+
+
//_________________________________________________________________________
Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t sigError, Int_t np, AliTrackPoint *points)
{
nPoints++;
}
}
- fitter->Eval();
+ if(fitter->Eval()) return 1.E10;
+
Double_t z0 = fitter->GetParameter(3);
Double_t dzdx = fitter->GetParameter(4);
Double_t c = fitter->GetParameter(2);
Double_t y0 = 1. / a;
Double_t x0 = -b * y0;
- Double_t R = TMath::Sqrt(y0*y0 + x0*x0 - c*y0);
+ Double_t tmp = y0*y0 + x0*x0 - c*y0;
+ if(tmp<=0.) return 1.E10;
+ Double_t R = TMath::Sqrt(tmp);
Double_t C = 1.0 + b*b - c*a;
if (C > 0.0) C = a / TMath::Sqrt(C);
if(!track){
for(Int_t ip = 0; ip < kNPlanes; ip++) {
x = tracklets[ip].GetX0();
- Double_t tmp = TMath::Sqrt(R*R-(x-x0)*(x-x0));
+ tmp = R*R-(x-x0)*(x-x0);
+ if(tmp <= 0.) continue;
+ tmp = TMath::Sqrt(tmp);
// y: R^2 = (x - x0)^2 + (y - y0)^2
// => y = y0 +/- Sqrt(R^2 - (x - x0)^2)
Float_t xyz[3];
for(int ip=0; ip<np; ip++){
points[ip].GetXYZ(xyz);
- xyz[1] = y0 - (y0>0.?1.:-1)*TMath::Sqrt(R*R-(xyz[0]-x0)*(xyz[0]-x0));
+ xyz[1] = y0 - (y0>0.?1.:-1.)*TMath::Sqrt(R*R-(xyz[0]-x0)*(xyz[0]-x0));
xyz[2] = z0 + dzdx * (xyz[0] - xref);
points[ip].SetXYZ(xyz);
}
}
- if(AliTRDReconstructor::StreamLevel() >=5){
- TTreeSRedirector &cstreamer = *fgDebugStreamer;
- Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
- Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
- Double_t chi2z = CalculateChi2Z(tracklets, z0, dzdx, xref);
- cstreamer << "FitRiemanTilt"
- << "EventNumber=" << eventNumber
- << "CandidateNumber=" << candidateNumber
- << "xref=" << xref
- << "Chi2Z=" << chi2z
- << "\n";
- }
return chi2;
}
+//____________________________________________________________________
+Double_t AliTRDtrackerV1::FitKalman(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t up, Int_t np, AliTrackPoint *points)
+{
+// Kalman filter implementation for the TRD.
+// It returns the positions of the fit in the array "points"
+//
+// Author : A.Bercuci@gsi.de
+
+ //printf("Start track @ x[%f]\n", track->GetX());
+
+ //prepare marker points along the track
+ Int_t ip = np ? 0 : 1;
+ while(ip<np){
+ if((up?-1:1) * (track->GetX() - points[ip].GetX()) > 0.) break;
+ //printf("AliTRDtrackerV1::FitKalman() : Skip track marker x[%d] = %7.3f. Before track start ( %7.3f ).\n", ip, points[ip].GetX(), track->GetX());
+ ip++;
+ }
+ //if(points) printf("First marker point @ x[%d] = %f\n", ip, points[ip].GetX());
+
+
+ AliTRDseedV1 tracklet, *ptrTracklet = 0x0;
+
+ //Loop through the TRD planes
+ for (Int_t jplane = 0; jplane < kNPlanes; jplane++) {
+ // GET TRACKLET OR BUILT IT
+ Int_t iplane = up ? jplane : kNPlanes - 1 - jplane;
+ if(tracklets){
+ if(!(ptrTracklet = &tracklets[iplane])) continue;
+ }else{
+ if(!(ptrTracklet = track->GetTracklet(iplane))){
+ /*AliTRDtrackerV1 *tracker = 0x0;
+ if(!(tracker = dynamic_cast<AliTRDtrackerV1*>( AliTRDReconstructor::Tracker()))) continue;
+ ptrTracklet = new(&tracklet) AliTRDseedV1(iplane);
+ if(!tracker->MakeTracklet(ptrTracklet, track)) */
+ continue;
+ }
+ }
+ if(!ptrTracklet->IsOK()) continue;
+
+ Double_t x = ptrTracklet->GetX0();
+
+ while(ip < np){
+ //don't do anything if next marker is after next update point.
+ if((up?-1:1) * (points[ip].GetX() - x) - fgkMaxStep < 0) break;
+
+ //printf("Propagate to x[%d] = %f\n", ip, points[ip].GetX());
+
+ if(((up?-1:1) * (points[ip].GetX() - track->GetX()) < 0) && !PropagateToX(*track, points[ip].GetX(), fgkMaxStep)) return -1.;
+
+ Double_t xyz[3]; // should also get the covariance
+ track->GetXYZ(xyz); points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]);
+ ip++;
+ }
+ //printf("plane[%d] tracklet[%p] x[%f]\n", iplane, ptrTracklet, x);
+
+ //Propagate closer to the next update point
+ if(((up?-1:1) * (x - track->GetX()) + fgkMaxStep < 0) && !PropagateToX(*track, x + (up?-1:1)*fgkMaxStep, fgkMaxStep)) return -1.;
+
+ if(!AdjustSector(track)) return -1;
+ if(TMath::Abs(track->GetSnp()) > fgkMaxSnp) return -1;
+
+ //load tracklet to the tracker and the track
+/* Int_t index;
+ if((index = FindTracklet(ptrTracklet)) < 0){
+ ptrTracklet = SetTracklet(&tracklet);
+ index = fTracklets->GetEntriesFast()-1;
+ }
+ track->SetTracklet(ptrTracklet, index);*/
+
+
+ // register tracklet to track with tracklet creation !!
+ // PropagateBack : loaded tracklet to the tracker and update index
+ // RefitInward : update index
+ // MakeTrack : loaded tracklet to the tracker and update index
+ if(!tracklets) track->SetTracklet(ptrTracklet, -1);
+
+
+ //Calculate the mean material budget along the path inside the chamber
+ Double_t xyz0[3]; track->GetXYZ(xyz0);
+ Double_t alpha = track->GetAlpha();
+ Double_t xyz1[3], y, z;
+ if(!track->GetProlongation(x, y, z)) return -1;
+ xyz1[0] = x * TMath::Cos(alpha) - y * TMath::Sin(alpha);
+ xyz1[1] = +x * TMath::Sin(alpha) + y * TMath::Cos(alpha);
+ xyz1[2] = z;
+ Double_t param[7];
+ if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) break;
+ Double_t xrho = param[0]*param[4]; // density*length
+ Double_t xx0 = param[1]; // radiation length
+
+ //Propagate the track
+ track->PropagateTo(x, xx0, xrho);
+ if (!AdjustSector(track)) break;
+
+ //Update track
+ Double_t chi2 = track->GetPredictedChi2(ptrTracklet);
+ if(chi2<1e+10) track->Update(ptrTracklet, chi2);
+
+ if(!up) continue;
+
+ //Reset material budget if 2 consecutive gold
+ if(iplane>0 && track->GetTracklet(iplane-1) && ptrTracklet->GetN() + track->GetTracklet(iplane-1)->GetN() > 20) track->SetBudget(2, 0.);
+ } // end planes loop
+
+ // extrapolation
+ while(ip < np){
+ if(((up?-1:1) * (points[ip].GetX() - track->GetX()) < 0) && !PropagateToX(*track, points[ip].GetX(), fgkMaxStep)) return -1.;
+
+ Double_t xyz[3]; // should also get the covariance
+ track->GetXYZ(xyz); points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]);
+ ip++;
+ }
+
+ return track->GetChi2();
+}
//_________________________________________________________________________
Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope, Double_t xref)
// Output: - The Chi2 value of the track in z-Direction
//
Float_t chi2Z = 0, nLayers = 0;
- for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNplan; iLayer++) {
+ 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);
// Calculate the mean material budget between start and
// end point of this prolongation step
- AliTracker::MeanMaterialBudget(xyz0, xyz1, param);
+ if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) return 0;
// Propagate the track to the X-position after the next step
if (!t.PropagateTo(x,param[1],param[0]*param[4])) {
branch->SetAddress(&clusterArray);
if(!fClusters){
- array = new TClonesArray("AliTRDcluster", nsize);
+ Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters();
+ if(fReconstructor->IsHLT()) nclusters /= AliTRDgeometry::kNsector;
+ array = new TClonesArray("AliTRDcluster", Int_t(nclusters));
array->SetOwner(kTRUE);
}
Int_t AliTRDtrackerV1::LoadClusters(TTree *cTree)
{
//
- // Fills clusters into TRD tracking_sectors
- // Note that the numbering scheme for the TRD tracking_sectors
- // differs from that of TRD sectors
+ // Fills clusters into TRD tracking sectors
//
-
- if (ReadClusters(fClusters, cTree)) {
- AliError("Problem with reading the clusters !");
+ if(!fReconstructor->IsWritingClusters()){
+ fClusters = AliTRDReconstructor::GetClusters();
+ } else {
+ if (ReadClusters(fClusters, cTree)) {
+ AliError("Problem with reading the clusters !");
+ return 1;
+ }
+ }
+ SetClustersOwner();
+
+ if(!fClusters || !fClusters->GetEntriesFast()){
+ AliInfo("No TRD clusters");
return 1;
}
- Int_t ncl = fClusters->GetEntriesFast(), nin = 0;
- if(!ncl){
- AliInfo("Clusters 0");
+
+ //Int_t nin =
+ BuildTrackingContainers();
+
+ //Int_t ncl = fClusters->GetEntriesFast();
+ //AliInfo(Form("Clusters %d [%6.2f %% in the active volume]", ncl, 100.*float(nin)/ncl));
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t AliTRDtrackerV1::LoadClusters(TClonesArray *clusters)
+{
+ //
+ // Fills clusters into TRD tracking sectors
+ // Function for use in the HLT
+
+ if(!clusters || !clusters->GetEntriesFast()){
+ AliInfo("No TRD clusters");
return 1;
}
- Int_t icl = ncl;
+ fClusters = clusters;
+ SetClustersOwner();
+
+ //Int_t nin =
+ BuildTrackingContainers();
+
+ //Int_t ncl = fClusters->GetEntriesFast();
+ //AliInfo(Form("Clusters %d [%6.2f %% in the active volume]", ncl, 100.*float(nin)/ncl));
+
+ return 0;
+}
+
+
+//____________________________________________________________________
+Int_t AliTRDtrackerV1::BuildTrackingContainers()
+{
+// Building tracking containers for clusters
+
+ Int_t nin =0, icl = fClusters->GetEntriesFast();
while (icl--) {
AliTRDcluster *c = (AliTRDcluster *) fClusters->UncheckedAt(icl);
if(c->IsInChamber()) nin++;
Int_t detector = c->GetDetector();
Int_t sector = fGeom->GetSector(detector);
- Int_t stack = fGeom->GetChamber(detector);
- Int_t plane = fGeom->GetPlane(detector);
+ Int_t stack = fGeom->GetStack(detector);
+ Int_t layer = fGeom->GetLayer(detector);
- fTrSec[sector].GetChamber(stack, plane, kTRUE)->InsertCluster(c, icl);
+ fTrSec[sector].GetChamber(stack, layer, kTRUE)->InsertCluster(c, icl);
}
- AliInfo(Form("Clusters %d in %6.2f %%", ncl, 100.*float(nin)/ncl));
-
- for(int isector =0; isector<AliTRDgeometry::kNsect; isector++){
+
+ const AliTRDCalDet *cal = AliTRDcalibDB::Instance()->GetT0Det();
+ for(int isector =0; isector<AliTRDgeometry::kNsector; isector++){
if(!fTrSec[isector].GetNChambers()) continue;
- fTrSec[isector].Init();
+ fTrSec[isector].Init(fReconstructor, cal);
}
-
- return 0;
+
+ return nin;
}
+
//____________________________________________________________________
void AliTRDtrackerV1::UnloadClusters()
{
if(fTracks) fTracks->Delete();
if(fTracklets) fTracklets->Delete();
- if(fClusters) fClusters->Delete();
+ if(fClusters){
+ if(IsClustersOwner()) fClusters->Delete();
+
+ // save clusters array in the reconstructor for further use.
+ if(!fReconstructor->IsWritingClusters()){
+ AliTRDReconstructor::SetClusters(fClusters);
+ SetClustersOwner(kFALSE);
+ } else AliTRDReconstructor::SetClusters(0x0);
+ }
- for (int i = 0; i < AliTRDgeometry::kNsect; i++) fTrSec[i].Clear();
+ for (int i = 0; i < AliTRDgeometry::kNsector; i++) fTrSec[i].Clear();
// Increment the Event Number
AliTRDtrackerDebug::SetEventNumber(AliTRDtrackerDebug::GetEventNumber() + 1);
// and adds the new tracklet to the list.
//
if(!fTracklets){
- fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsect()*kMaxTracksStack);
+ fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsector()*kMaxTracksStack);
fTracklets->SetOwner(kTRUE);
}
Int_t nentries = fTracklets->GetEntriesFast();
return new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet);
}
+//____________________________________________________________________
+AliTRDtrackV1* AliTRDtrackerV1::SetTrack(AliTRDtrackV1 *track)
+{
+ // Add this track to the list of tracks stored in the tracker
+ //
+ // Parameters
+ // - track : pointer to the track to be added to the list
+ //
+ // Output
+ // - the pointer added
+ //
+ // Detailed description
+ // Build the tracks list if it is not yet created (late initialization)
+ // and adds the new track to the list.
+ //
+ if(!fTracks){
+ fTracks = new TClonesArray("AliTRDtrackV1", AliTRDgeometry::Nsector()*kMaxTracksStack);
+ fTracks->SetOwner(kTRUE);
+ }
+ Int_t nentries = fTracks->GetEntriesFast();
+ return new ((*fTracks)[nentries]) AliTRDtrackV1(*track);
+}
+
+
//____________________________________________________________________
Int_t AliTRDtrackerV1::Clusters2TracksSM(Int_t sector, AliESDEvent *esd)
Int_t nTracks = 0;
Int_t nChambers = 0;
AliTRDtrackingChamber **stack = 0x0, *chamber = 0x0;
- for(int istack = 0; istack<AliTRDgeometry::kNcham; istack++){
+ for(int istack = 0; istack<AliTRDgeometry::kNstack; istack++){
if(!(stack = fTrSec[sector].GetStack(istack))) continue;
nChambers = 0;
- for(int iplane=0; iplane<AliTRDgeometry::kNplan; iplane++){
- if(!(chamber = stack[iplane])) continue;
- if(chamber->GetNClusters() < fgNTimeBins * AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue;
+ for(int ilayer=0; ilayer<AliTRDgeometry::kNlayer; ilayer++){
+ if(!(chamber = stack[ilayer])) continue;
+ if(chamber->GetNClusters() < fgNTimeBins * fReconstructor->GetRecoParam() ->GetFindableClusters()) continue;
nChambers++;
- //AliInfo(Form("sector %d stack %d plane %d clusters %d", sector, istack, iplane, chamber->GetNClusters()));
+ //AliInfo(Form("sector %d stack %d layer %d clusters %d", sector, istack, ilayer, chamber->GetNClusters()));
}
if(nChambers < 4) continue;
//AliInfo(Form("Doing stack %d", istack));
// 8. Build ESD track and register it to the output list
//
+ const AliTRDCalDet *cal = AliTRDcalibDB::Instance()->GetT0Det();
AliTRDtrackingChamber *chamber = 0x0;
AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized
Int_t pars[4]; // MakeSeeds parameters
// Build initial seeding configurations
Double_t quality = BuildSeedingConfigs(stack, configs);
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
AliInfo(Form("Plane config %d %d %d Quality %f"
, configs[0], configs[1], configs[2], quality));
}
+
// Initialize contors
Int_t ntracks, // number of TRD track candidates
ntracks1, // number of registered TRD tracks/iter
ntracks2 = 0; // number of all registered TRD tracks in stack
fSieveSeeding = 0;
+
+ // Get stack index
+ Int_t ic = 0; AliTRDtrackingChamber **cIter = &stack[0];
+ while(ic<kNPlanes && !(*cIter)){ic++; cIter++;}
+ if(!(*cIter)) return ntracks2;
+ Int_t istack = fGeom->GetStack((*cIter)->GetDetector());
+
do{
// Loop over seeding configurations
ntracks = 0; ntracks1 = 0;
for (Int_t iconf = 0; iconf<3; iconf++) {
pars[0] = configs[iconf];
pars[1] = ntracks;
+ pars[2] = istack;
ntracks = MakeSeeds(stack, &sseed[6*ntracks], pars);
if(ntracks == kMaxTracksStack) break;
}
- if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding));
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding));
if(!ntracks) break;
// Check track candidates
candidates = 0;
for (Int_t itrack = 0; itrack < ntracks; itrack++) {
- Int_t trackIndex = sort[itrack];
- if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
+ Int_t trackIndex = sort[itrack];
+ if (signedTrack[trackIndex] || fakeTrack[trackIndex]) continue;
- // Calculate track parameters from tracklets seeds
- Int_t labelsall[1000];
- Int_t nlabelsall = 0;
- Int_t naccepted = 0;
- Int_t ncl = 0;
- Int_t nused = 0;
- Int_t nlayers = 0;
- Int_t findable = 0;
- 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++;
-
- sseed[jseed].UpdateUsed();
- ncl += sseed[jseed].GetN2();
- nused += sseed[jseed].GetNUsed();
- nlayers++;
-
- // 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++;
- }
- }
- }
+ // Calculate track parameters from tracklets seeds
+ Int_t ncl = 0;
+ Int_t nused = 0;
+ Int_t nlayers = 0;
+ Int_t findable = 0;
+ 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++;
+
+ 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);
}
signedTrack[trackIndex] = kTRUE;
-
- // Build track label - what happens if measured data ???
- Int_t labels[1000];
- Int_t outlab[1000];
- Int_t nlab = 0;
- for (Int_t iLayer = 0; iLayer < 6; iLayer++) {
- Int_t jseed = kNPlanes*trackIndex+iLayer;
- 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++;
- }
- }
- 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);
-
// Sign clusters
AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1;
- for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
+ 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){
- Int_t ic = 0;
+ ic = 0;
while(!(cl = sseed[jseed].GetClusters(ic))) ic++;
clusterIndex = sseed[jseed].GetIndexes(ic);
}
// Build track parameters
AliTRDseedV1 *lseed =&sseed[trackIndex*6];
- Int_t idx = 0;
+/* Int_t idx = 0;
while(idx<3 && !lseed->IsOK()) {
idx++;
lseed++;
- }
- Double_t cR = lseed->GetC();
- Double_t x = lseed->GetX0() - 3.5;
+ }*/
+ Double_t x = lseed->GetX0();// - 3.5;
trackParams[0] = x; //NEW AB
- trackParams[1] = lseed->GetYat(x);//lseed->GetYref(0);
- trackParams[2] = lseed->GetZat(x);//lseed->GetZref(0);
- trackParams[3] = lseed->GetX0() * cR - TMath::Sin(TMath::ATan(lseed->GetYref(1)));
+ 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] = cR;
+ trackParams[5] = lseed->GetC();
Int_t ich = 0; while(!(chamber = stack[ich])) ich++;
trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *alpha+shift; // Supermodule*/
- if(AliTRDReconstructor::StreamLevel() > 1){
+ 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];
- for(int is=0; is<6; is++){
- dseed[is] = new AliTRDseedV1(sseed[trackIndex*6+is]);
- dseed[is]->SetOwner();
- nclusters += sseed[is].GetN2();
+
+ // 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++;
+ }
+ }
}
+ 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();
<< "Ncl=" << ncl
<< "NLayers=" << nlayers
<< "Findable=" << findable
-
<< "NUsed=" << nused
<< "\n";
}
AliWarning("Fail to build a TRD Track.");
continue;
}
+
//AliInfo("End of MakeTrack()");
- AliESDtrack esdTrack;
- esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout);
- esdTrack.SetLabel(track->GetLabel());
- track->UpdateESDtrack(&esdTrack);
+ AliESDtrack *esdTrack = new ((*esdTrackList)[ntracks0++]) AliESDtrack();
+ esdTrack->UpdateTrackParams(track, AliESDtrack::kTRDout);
+ esdTrack->SetLabel(track->GetLabel());
+ track->UpdateESDtrack(esdTrack);
// write ESD-friends if neccessary
- if (AliTRDReconstructor::StreamLevel() > 0){
- //printf("Creating Calibrations Object\n");
+ if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0){
AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(*track);
calibTrack->SetOwner();
- esdTrack.AddCalibObject(calibTrack);
+ esdTrack->AddCalibObject(calibTrack);
}
- new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack);
ntracks1++;
AliTRDtrackerDebug::SetTrackNumber(AliTRDtrackerDebug::GetTrackNumber() + 1);
}
// increment counters
ntracks2 += ntracks1;
+
+ if(fReconstructor->IsHLT()) break;
fSieveSeeding++;
// Rebuild plane configurations and indices taking only unused clusters into account
quality = BuildSeedingConfigs(stack, configs);
- if(quality < 1.E-7) break; //AliTRDReconstructor::RecoParam()->GetPlaneQualityThreshold()) break;
+ if(quality < 1.E-7) break; //fReconstructor->GetRecoParam() ->GetPlaneQualityThreshold()) break;
for(Int_t ip = 0; ip < kNPlanes; ip++){
if(!(chamber = stack[ip])) continue;
- chamber->Build(fGeom);//Indices(fSieveSeeding);
+ chamber->Build(fGeom, cal);//Indices(fSieveSeeding);
}
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
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
//
AliTRDtrackingChamber *chamber = 0x0;
- AliTRDcluster *c[4] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
+ AliTRDcluster *c[kNSeedPlanes] = {0x0, 0x0, 0x0, 0x0}; // initilize seeding clusters
AliTRDseedV1 *cseed = &sseed[0]; // initialize tracklets for first track
Int_t ncl, mcl; // working variable for looping over clusters
Int_t index[AliTRDchamberTimeBin::kMaxClustersLayer], jndex[AliTRDchamberTimeBin::kMaxClustersLayer];
// 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
Double_t seedQuality[kMaxTracksStack];
// unpack control parameters
Int_t config = ipar[0];
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);
+
+
// Init chambers geometry
- Int_t ic = 0; while(!(chamber = stack[ic])) ic++;
- Int_t istack = fGeom->GetChamber(chamber->GetDetector());
Double_t hL[kNPlanes]; // Tilting angle
Float_t padlength[kNPlanes]; // pad lenghts
AliTRDpadPlane *pp = 0x0;
for(int iplane=0; iplane<kNPlanes; iplane++){
pp = fGeom->GetPadPlane(iplane, istack);
- hL[iplane] = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle());
+ hL[iplane] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle());
padlength[iplane] = pp->GetLengthIPad();
}
- if(AliTRDReconstructor::StreamLevel() > 1){
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks));
}
+ // Build seeding layers
+ ResetSeedTB();
Int_t nlayers = 0;
- AliTRDchamberTimeBin *layer[] = {0x0, 0x0, 0x0, 0x0};
for(int isl=0; isl<kNSeedPlanes; isl++){
if(!(chamber = stack[planes[isl]])) continue;
- if(!(layer[isl] = chamber->GetSeedingLayer(fGeom))) continue;
+ if(!chamber->GetSeedingLayer(fSeedTB[isl], fGeom, fReconstructor)) continue;
nlayers++;
- //AliInfo(Form("seeding plane %d clusters %d", planes[isl], Int_t(*layer[isl])));
}
if(nlayers < 4) return 0;
// Start finding seeds
Double_t cond0[4], cond1[4], cond2[4];
Int_t icl = 0;
- while((c[3] = (*layer[3])[icl++])){
+ while((c[3] = (*fSeedTB[3])[icl++])){
if(!c[3]) continue;
- layer[0]->BuildCond(c[3], cond0, 0);
- layer[0]->GetClusters(cond0, index, ncl);
+ fSeedTB[0]->BuildCond(c[3], cond0, 0);
+ fSeedTB[0]->GetClusters(cond0, index, ncl);
//printf("Found c[3] candidates 0 %d\n", ncl);
Int_t jcl = 0;
while(jcl<ncl) {
- c[0] = (*layer[0])[index[jcl++]];
+ 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;
- layer[1]->BuildCond(c[0], cond1, 1, theta, phi);
- layer[1]->GetClusters(cond1, jndex, mcl);
+ fSeedTB[1]->BuildCond(c[0], cond1, 1, theta, phi);
+ fSeedTB[1]->GetClusters(cond1, jndex, mcl);
//printf("Found c[0] candidates 1 %d\n", mcl);
Int_t kcl = 0;
while(kcl<mcl) {
- c[1] = (*layer[1])[jndex[kcl++]];
- if(!c[1]) continue;
- layer[2]->BuildCond(c[1], cond2, 2, theta, phi);
- c[2] = layer[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());
-
- for (Int_t il = 0; il < 6; il++) cseed[il].Reset();
-
- FitRieman(c, chi2);
-
- AliTRDseedV1 *tseed = 0x0;
- for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
- Int_t jLayer = planes[iLayer];
- tseed = &cseed[jLayer];
- tseed->SetPlane(jLayer);
- tseed->SetTilt(hL[jLayer]);
- tseed->SetPadLength(padlength[jLayer]);
- tseed->SetX0(stack[jLayer]->GetX());
- tseed->Init(GetRiemanFitter());
- }
-
- Bool_t isFake = kFALSE;
- if(AliTRDReconstructor::StreamLevel() >= 2){
- if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
- if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
- if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
-
- Double_t xpos[4];
- for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX();
- Float_t yref[4];
- for(int il=0; il<4; il++) yref[il] = cseed[planes[il]].GetYref(0);
- Int_t ll = c[3]->GetLabel(0);
- Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
- Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
- AliRieman *rim = GetRiemanFitter();
- TTreeSRedirector &cs0 = *fgDebugStreamer;
- cs0 << "MakeSeeds0"
- <<"EventNumber=" << eventNumber
- <<"CandidateNumber=" << candidateNumber
- <<"isFake=" << isFake
- <<"config=" << config
- <<"label=" << ll
- <<"chi2z=" << chi2[0]
- <<"chi2y=" << chi2[1]
- <<"Y2exp=" << cond2[0]
- <<"Z2exp=" << cond2[1]
- <<"X0=" << xpos[0] //layer[sLayer]->GetX()
- <<"X1=" << xpos[1] //layer[sLayer + 1]->GetX()
- <<"X2=" << xpos[2] //layer[sLayer + 2]->GetX()
- <<"X3=" << xpos[3] //layer[sLayer + 3]->GetX()
- <<"yref0=" << yref[0]
- <<"yref1=" << yref[1]
- <<"yref2=" << yref[2]
- <<"yref3=" << yref[3]
- <<"c0.=" << c[0]
- <<"c1.=" << c[1]
- <<"c2.=" << c[2]
- <<"c3.=" << c[3]
- <<"Seed0.=" << &cseed[planes[0]]
- <<"Seed1.=" << &cseed[planes[1]]
- <<"Seed2.=" << &cseed[planes[2]]
- <<"Seed3.=" << &cseed[planes[3]]
- <<"RiemanFitter.=" << rim
- <<"\n";
- }
-
- if(chi2[0] > AliTRDReconstructor::RecoParam()->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
- //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0]));
- AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
- continue;
- }
- if(chi2[1] > AliTRDReconstructor::RecoParam()->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){
- //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;
- for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
- Int_t jLayer = planes[iLayer];
- if(!cseed[jLayer].AttachClustersIter(stack[jLayer], 5., kFALSE, c[iLayer])) continue;
- nUsedCl += cseed[jLayer].GetNUsed();
- if(nUsedCl > 25) break;
- mlayers++;
- }
- if(mlayers < kNSeedPlanes){
- //AliInfo(Form("Failed updating all seeds %d [%d].", mlayers, kNSeedPlanes));
- AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
- continue;
- }
- // fit tracklets and cook likelihood
- FitTiltedRieman(&cseed[0], kTRUE);// Update Seeds and calculate Likelihood
- chi2[0] = GetChi2Y(&cseed[0]);
- chi2[1] = GetChi2Z(&cseed[0]);
- //Chi2 definitions in testing stage
- //chi2[0] = GetChi2YTest(&cseed[0]);
- //chi2[1] = GetChi2ZTest(&cseed[0]);
- Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked
-
- if (TMath::Log(1.E-9 + like) < AliTRDReconstructor::RecoParam()->GetTrackLikelihood()){
- //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
- AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
- continue;
- }
- //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
+ c[1] = (*fSeedTB[1])[jndex[kcl++]];
+ if(!c[1]) continue;
+ fSeedTB[2]->BuildCond(c[1], cond2, 2, theta, phi);
+ 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());
+
+ for (Int_t il = 0; il < kNPlanes; il++) cseed[il].Reset();
+
+ FitRieman(c, chi2);
+
+ AliTRDseedV1 *tseed = &cseed[0];
+ AliTRDtrackingChamber **cIter = &stack[0];
+ for(int iLayer=0; iLayer<kNPlanes; iLayer++, tseed++, cIter++){
+ tseed->SetDetector((*cIter) ? (*cIter)->GetDetector() : -1);
+ tseed->SetTilt(hL[iLayer]);
+ tseed->SetPadLength(padlength[iLayer]);
+ tseed->SetReconstructor(fReconstructor);
+ tseed->SetX0((*cIter) ? (*cIter)->GetX() : x_def[iLayer]);
+ tseed->Init(GetRiemanFitter());
+ }
+
+ Bool_t isFake = kFALSE;
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){
+ if (c[0]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
+ if (c[1]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
+ if (c[2]->GetLabel(0) != c[3]->GetLabel(0)) isFake = kTRUE;
+
+ Double_t xpos[4];
+ for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = fSeedTB[l]->GetX();
+ Float_t yref[4];
+ for(int il=0; il<4; il++) yref[il] = cseed[planes[il]].GetYref(0);
+ Int_t ll = c[3]->GetLabel(0);
+ Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+ Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+ AliRieman *rim = GetRiemanFitter();
+ TTreeSRedirector &cs0 = *fgDebugStreamer;
+ cs0 << "MakeSeeds0"
+ <<"EventNumber=" << eventNumber
+ <<"CandidateNumber=" << candidateNumber
+ <<"isFake=" << isFake
+ <<"config=" << config
+ <<"label=" << ll
+ <<"chi2z=" << chi2[0]
+ <<"chi2y=" << chi2[1]
+ <<"Y2exp=" << cond2[0]
+ <<"Z2exp=" << cond2[1]
+ <<"X0=" << xpos[0] //layer[sLayer]->GetX()
+ <<"X1=" << xpos[1] //layer[sLayer + 1]->GetX()
+ <<"X2=" << xpos[2] //layer[sLayer + 2]->GetX()
+ <<"X3=" << xpos[3] //layer[sLayer + 3]->GetX()
+ <<"yref0=" << yref[0]
+ <<"yref1=" << yref[1]
+ <<"yref2=" << yref[2]
+ <<"yref3=" << yref[3]
+ <<"c0.=" << c[0]
+ <<"c1.=" << c[1]
+ <<"c2.=" << c[2]
+ <<"c3.=" << c[3]
+ <<"Seed0.=" << &cseed[planes[0]]
+ <<"Seed1.=" << &cseed[planes[1]]
+ <<"Seed2.=" << &cseed[planes[2]]
+ <<"Seed3.=" << &cseed[planes[3]]
+ <<"RiemanFitter.=" << rim
+ <<"\n";
+ }
+ if(chi2[0] > fReconstructor->GetRecoParam() ->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){
+// //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]));
+ AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
+ continue;
+ }
+ //AliInfo("Passed chi2 filter.");
+
+ // try attaching clusters to tracklets
+ Int_t nUsedCl = 0;
+ Int_t mlayers = 0;
+ for(int iLayer=0; iLayer<kNSeedPlanes; iLayer++){
+ Int_t jLayer = planes[iLayer];
+ if(!cseed[jLayer].AttachClustersIter(stack[jLayer], 5., kFALSE, c[iLayer])) continue;
+ nUsedCl += cseed[jLayer].GetNUsed();
+ if(nUsedCl > 25) break;
+ mlayers++;
+ }
+
+ if(mlayers < kNSeedPlanes){
+ //AliInfo(Form("Failed updating all seeds %d [%d].", mlayers, kNSeedPlanes));
+ AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
+ continue;
+ }
- // book preliminary results
- seedQuality[ntracks] = like;
- fSeedLayer[ntracks] = config;/*sLayer;*/
+ // temporary exit door for the HLT
+ if(fReconstructor->IsHLT()){
+ // attach clusters to extrapolation chambers
+ for(int iLayer=0; iLayer<kNPlanes-kNSeedPlanes; iLayer++){
+ Int_t jLayer = planesExt[iLayer];
+ if(!(chamber = stack[jLayer])) continue;
+ cseed[jLayer].AttachClustersIter(chamber, 1000.);
+ }
+ fTrackQuality[ntracks] = 1.; // dummy value
+ ntracks++;
+ if(ntracks == kMaxTracksStack) return ntracks;
+ cseed += 6;
+ continue;
+ }
- // attach clusters to the extrapolation seeds
- Int_t lextrap[2];
- GetExtrapolationConfig(config, lextrap);
- Int_t nusedf = 0; // debug value
- for(int iLayer=0; iLayer<2; iLayer++){
- Int_t jLayer = lextrap[iLayer];
- if(!(chamber = stack[jLayer])) continue;
-
- // prepare extrapolated seed
- cseed[jLayer].Reset();
- cseed[jLayer].SetPlane(jLayer);
- cseed[jLayer].SetTilt(hL[jLayer]);
- cseed[jLayer].SetX0(chamber->GetX());
- cseed[jLayer].SetPadLength(padlength[jLayer]);
-
- // fit extrapolated seed
- if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
- if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
- AliTRDseedV1 pseed = cseed[jLayer];
- if(!pseed.AttachClustersIter(chamber, 1000.)) continue;
- cseed[jLayer] = pseed;
- nusedf += cseed[jLayer].GetNUsed(); // debug value
- FitTiltedRieman(cseed, kTRUE);
- }
-
- // AliInfo("Extrapolation done.");
- // Debug Stream containing all the 6 tracklets
- if(AliTRDReconstructor::StreamLevel() >= 2){
- TTreeSRedirector &cstreamer = *fgDebugStreamer;
- TLinearFitter *tiltedRieman = GetTiltedRiemanFitter();
- Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
- Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
- cstreamer << "MakeSeeds1"
- << "EventNumber=" << eventNumber
- << "CandidateNumber=" << candidateNumber
- << "S0.=" << &cseed[0]
- << "S1.=" << &cseed[1]
- << "S2.=" << &cseed[2]
- << "S3.=" << &cseed[3]
- << "S4.=" << &cseed[4]
- << "S5.=" << &cseed[5]
- << "FitterT.=" << tiltedRieman
- << "\n";
- }
-
- if(ImproveSeedQuality(stack, cseed) < 4){
- AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
- continue;
- }
- //AliInfo("Improve seed quality done.");
-
- // fit full track and cook likelihoods
- // Double_t curv = FitRieman(&cseed[0], chi2);
- // Double_t chi2ZF = chi2[0] / TMath::Max((mlayers - 3.), 1.);
- // Double_t chi2RF = chi2[1] / TMath::Max((mlayers - 3.), 1.);
-
- // do the final track fitting (Once with vertex constraint and once without vertex constraint)
- Double_t chi2Vals[3];
- chi2Vals[0] = FitTiltedRieman(&cseed[0], kFALSE);
- chi2Vals[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ());
- chi2Vals[2] = GetChi2Z(&cseed[0]) / TMath::Max((mlayers - 3.), 1.);
- // Chi2 definitions in testing stage
- //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]);
- }
+ // fit tracklets and cook likelihood
+ FitTiltedRieman(&cseed[0], kTRUE);// Update Seeds and calculate Likelihood
+ Double_t like = CookLikelihood(&cseed[0], planes); // to be checked
- if(AliTRDReconstructor::StreamLevel() >= 2){
- TTreeSRedirector &cstreamer = *fgDebugStreamer;
- Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
- Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
- TLinearFitter *fitterTC = GetTiltedRiemanFitterConstraint();
- TLinearFitter *fitterT = GetTiltedRiemanFitter();
- cstreamer << "MakeSeeds2"
- << "EventNumber=" << eventNumber
- << "CandidateNumber=" << candidateNumber
- << "Chi2TR=" << chi2Vals[0]
- << "Chi2TC=" << chi2Vals[1]
- << "Nlayers=" << mlayers
- << "NUsedS=" << nUsedCl
- << "NUsed=" << nusedf
- << "Like=" << like
- << "S0.=" << &cseed[0]
- << "S1.=" << &cseed[1]
- << "S2.=" << &cseed[2]
- << "S3.=" << &cseed[3]
- << "S4.=" << &cseed[4]
- << "S5.=" << &cseed[5]
- << "Label=" << label
- << "Freq=" << frequency
- << "FitterT.=" << fitterT
- << "FitterTC.=" << fitterTC
- << "\n";
- }
-
- ntracks++;
- AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
- if(ntracks == kMaxTracksStack){
- AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
- for(int isl=0; isl<4; isl++) delete layer[isl];
- return ntracks;
- }
- cseed += 6;
+ if (TMath::Log(1.E-9 + like) < fReconstructor->GetRecoParam() ->GetTrackLikelihood()){
+ //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
+ AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
+ continue;
+ }
+ //AliInfo(Form("Passed likelihood %f[%e].", TMath::Log(1.E-9 + like), like));
+
+ // book preliminary results
+ seedQuality[ntracks] = like;
+ fSeedLayer[ntracks] = config;/*sLayer;*/
+
+ // attach clusters to the extrapolation seeds
+ Int_t nusedf = 0; // debug value
+ for(int iLayer=0; iLayer<kNPlanes-kNSeedPlanes; iLayer++){
+ Int_t jLayer = planesExt[iLayer];
+ if(!(chamber = stack[jLayer])) continue;
+
+ // fit extrapolated seed
+ if ((jLayer == 0) && !(cseed[1].IsOK())) continue;
+ if ((jLayer == 5) && !(cseed[4].IsOK())) continue;
+ AliTRDseedV1 pseed = cseed[jLayer];
+ if(!pseed.AttachClustersIter(chamber, 1000.)) continue;
+ cseed[jLayer] = pseed;
+ nusedf += cseed[jLayer].GetNUsed(); // debug value
+ FitTiltedRieman(cseed, kTRUE);
+ }
+
+ // AliInfo("Extrapolation done.");
+ // Debug Stream containing all the 6 tracklets
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){
+ TTreeSRedirector &cstreamer = *fgDebugStreamer;
+ TLinearFitter *tiltedRieman = GetTiltedRiemanFitter();
+ Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+ Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+ cstreamer << "MakeSeeds1"
+ << "EventNumber=" << eventNumber
+ << "CandidateNumber=" << candidateNumber
+ << "S0.=" << &cseed[0]
+ << "S1.=" << &cseed[1]
+ << "S2.=" << &cseed[2]
+ << "S3.=" << &cseed[3]
+ << "S4.=" << &cseed[4]
+ << "S5.=" << &cseed[5]
+ << "FitterT.=" << tiltedRieman
+ << "\n";
+ }
+
+ if(fReconstructor->GetRecoParam()->HasImproveTracklets() && ImproveSeedQuality(stack, cseed) < 4){
+ AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
+ continue;
+ }
+ //AliInfo("Improve seed quality done.");
+
+ // fit full track and cook likelihoods
+ // Double_t curv = FitRieman(&cseed[0], chi2);
+ // Double_t chi2ZF = chi2[0] / TMath::Max((mlayers - 3.), 1.);
+ // Double_t chi2RF = chi2[1] / TMath::Max((mlayers - 3.), 1.);
+
+ // 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[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ()); // Do Vertex Constrained fit if desired
+ else
+ chi2Vals[1] = 1.;
+ chi2Vals[2] = GetChi2Z(&cseed[0]) / TMath::Max((mlayers - 3.), 1.);
+ // Chi2 definitions in testing stage
+ //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 = *fgDebugStreamer;
+ Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+ Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+ TLinearFitter *fitterTC = GetTiltedRiemanFitterConstraint();
+ TLinearFitter *fitterT = GetTiltedRiemanFitter();
+ Int_t ncls = 0;
+ for(Int_t iseed = 0; iseed < kNPlanes; iseed++){
+ ncls += cseed[iseed].IsOK() ? cseed[iseed].GetN2() : 0;
+ }
+ cstreamer << "MakeSeeds2"
+ << "EventNumber=" << eventNumber
+ << "CandidateNumber=" << candidateNumber
+ << "Chi2TR=" << chi2Vals[0]
+ << "Chi2TC=" << chi2Vals[1]
+ << "Nlayers=" << mlayers
+ << "NClusters=" << ncls
+ << "NUsedS=" << nUsedCl
+ << "NUsed=" << nusedf
+ << "Like=" << like
+ << "S0.=" << &cseed[0]
+ << "S1.=" << &cseed[1]
+ << "S2.=" << &cseed[2]
+ << "S3.=" << &cseed[3]
+ << "S4.=" << &cseed[4]
+ << "S5.=" << &cseed[5]
+ << "Label=" << label
+ << "Freq=" << frequency
+ << "FitterT.=" << fitterT
+ << "FitterTC.=" << fitterTC
+ << "\n";
+ }
+
+ ntracks++;
+ AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1);
+ if(ntracks == kMaxTracksStack){
+ AliWarning(Form("Number of seeds reached maximum allowed (%d) in stack.", kMaxTracksStack));
+ return ntracks;
+ }
+ cseed += 6;
}
}
}
- for(int isl=0; isl<4; isl++) delete layer[isl];
return ntracks;
}
// To be discussed with Marian !!
//
- AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
- if (!calibra) AliInfo("Could not get Calibra instance\n");
Double_t alpha = AliTRDgeometry::GetAlpha();
Double_t shift = AliTRDgeometry::GetAlpha()/2.0;
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;
- AliTRDtrackV1 *track = new AliTRDtrackV1(seeds, ¶ms[1], c, params[0], params[6]*alpha+shift);
- track->PropagateTo(params[0]-5.0);
- track->ResetCovariance(1);
- Int_t nc = FollowBackProlongation(*track);
- if (nc < 30) {
- delete track;
- track = 0x0;
- } else {
- //track->CookdEdx();
- //track->CookdEdxTimBin(-1);
- track->CookLabel(.9);
- // computes PID for track
- track->CookPID();
- // update calibration references using this track
- if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(track);
+ AliTRDtrackV1 track(seeds, ¶ms[1], c, params[0], params[6]*alpha+shift);
+ track.PropagateTo(params[0]-5.0);
+ if(fReconstructor->IsHLT()){
+ AliTRDseedV1 *ptrTracklet = 0x0;
+ for(Int_t ip=0; ip<kNPlanes; ip++){
+ track.UnsetTracklet(ip);
+ ptrTracklet = SetTracklet(&seeds[ip]);
+ track.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1);
+ }
+ return SetTrack(&track);
+ }
+
+ track.ResetCovariance(1);
+ Int_t nc = TMath::Abs(FollowBackProlongation(track));
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 5){
+ Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+ Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+ Double_t p[5]; // Track Params for the Debug Stream
+ track.GetExternalParameters(params[0], p);
+ TTreeSRedirector &cs = *fgDebugStreamer;
+ cs << "MakeTrack"
+ << "EventNumber=" << eventNumber
+ << "CandidateNumber=" << candidateNumber
+ << "nc=" << nc
+ << "X=" << params[0]
+ << "Y=" << p[0]
+ << "Z=" << p[1]
+ << "snp=" << p[2]
+ << "tnd=" << p[3]
+ << "crv=" << p[4]
+ << "Yin=" << params[1]
+ << "Zin=" << params[2]
+ << "snpin=" << params[3]
+ << "tndin=" << params[4]
+ << "crvin=" << params[5]
+ << "track.=" << &track
+ << "\n";
}
+ if (nc < 30) return 0x0;
- return track;
+ AliTRDtrackV1 *ptrTrack = SetTrack(&track);
+ ptrTrack->SetReconstructor(fReconstructor);
+ ptrTrack->CookLabel(.9);
+
+ // computes PID for track
+ ptrTrack->CookPID();
+ // update calibration references using this track
+ AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
+ if (!calibra){
+ AliInfo("Could not get Calibra instance\n");
+ if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(ptrTrack);
+ }
+ return ptrTrack;
}
Int_t sortindexes[6];
for (Int_t jLayer = 0; jLayer < 6; jLayer++) {
- squality[jLayer] = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.;
+ squality[jLayer] = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : 1000.;
sumquality += squality[jLayer];
}
if ((sumquality >= lastquality) || (chi2 > lastchi2)) break;
}
chi2 = FitTiltedRieman(bseed, kTRUE);
- if(AliTRDReconstructor::StreamLevel() >= 7){
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 7){
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
TLinearFitter *tiltedRieman = GetTiltedRiemanFitter();
sumdaf /= Float_t (nLayers - 2.0);
Double_t likeChi2Z = TMath::Exp(-chi2[2] * 0.14); // Chi2Z
- Double_t likeChi2TC = TMath::Exp(-chi2[1] * 0.677); // Constrained Tilted Riemann
+ Double_t likeChi2TC = (fReconstructor->GetRecoParam() ->IsVertexConstrained()) ?
+ 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;
- if(AliTRDReconstructor::StreamLevel() >= 2){
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
TTreeSRedirector &cstreamer = *fgDebugStreamer;
}
//____________________________________________________________________
-Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4]
- , Double_t *chi2)
+Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4])
{
//
// Calculate the probability of this track candidate.
//
// ratio of the total number of clusters/track which are expected to be found by the tracker.
- Float_t fgFindable = AliTRDReconstructor::RecoParam()->GetFindableClusters();
-
+ const AliTRDrecoParam *fRecoPars = fReconstructor->GetRecoParam();
- Int_t nclusters = 0;
+ Double_t chi2y = GetChi2Y(&cseed[0]);
+ Double_t chi2z = GetChi2Z(&cseed[0]);
+
+ Float_t nclusters = 0.;
Double_t sumda = 0.;
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));
}
- Double_t likea = TMath::Exp(-sumda*10.6);
+ nclusters *= .25;
+
+ Double_t likea = TMath::Exp(-sumda * fRecoPars->GetPhiCut());
Double_t likechi2y = 0.0000000001;
- if (chi2[0] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[0]) * 7.73);
- Double_t likechi2z = TMath::Exp(-chi2[1] * 0.088) / TMath::Exp(-chi2[1] * 0.019);
- Int_t enc = Int_t(fgFindable*4.*fgNTimeBins); // Expected Number Of Clusters, normally 72
- Double_t likeN = TMath::Exp(-(enc - nclusters) * 0.19);
-
+ if (fReconstructor->IsCosmic() || chi2y < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2y) * fRecoPars->GetChi2YCut());
+ Double_t likechi2z = TMath::Exp(-chi2z * fRecoPars->GetChi2ZCut());
+ 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(AliTRDReconstructor::StreamLevel() >= 2){
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){
Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber();
+ Int_t nTracklets = 0; Float_t mean_ncls = 0;
+ for(Int_t iseed=0; iseed < kNPlanes; iseed++){
+ if(!cseed[iseed].IsOK()) continue;
+ nTracklets++;
+ mean_ncls += cseed[iseed].GetN2();
+ }
+ if(nTracklets) mean_ncls /= nTracklets;
// The Debug Stream contains the seed
TTreeSRedirector &cstreamer = *fgDebugStreamer;
cstreamer << "CookLikelihood"
<< "tracklet4.=" << &cseed[4]
<< "tracklet5.=" << &cseed[5]
<< "sumda=" << sumda
- << "chi0=" << chi2[0]
- << "chi1=" << chi2[1]
+ << "chi2y=" << chi2y
+ << "chi2z=" << chi2z
<< "likea=" << likea
<< "likechi2y=" << likechi2y
<< "likechi2z=" << likechi2z
<< "nclusters=" << nclusters
<< "likeN=" << likeN
<< "like=" << like
+ << "meanncls=" << mean_ncls
<< "\n";
}
AliCluster* AliTRDtrackerV1::GetCluster(Int_t idx) const
{
Int_t ncls = fClusters->GetEntriesFast();
- return idx >= 0 || idx < ncls ? (AliCluster*)fClusters->UncheckedAt(idx) : 0x0;
+ return idx >= 0 && idx < ncls ? (AliCluster*)fClusters->UncheckedAt(idx) : 0x0;
}
//____________________________________________________________________
AliTRDseedV1* AliTRDtrackerV1::GetTracklet(Int_t idx) const
{
Int_t ntrklt = fTracklets->GetEntriesFast();
- return idx >= 0 || idx < ntrklt ? (AliTRDseedV1*)fTracklets->UncheckedAt(idx) : 0x0;
+ return idx >= 0 && idx < ntrklt ? (AliTRDseedV1*)fTracklets->UncheckedAt(idx) : 0x0;
}
//____________________________________________________________________
AliKalmanTrack* AliTRDtrackerV1::GetTrack(Int_t idx) const
{
Int_t ntrk = fTracks->GetEntriesFast();
- return idx >= 0 || idx < ntrk ? (AliKalmanTrack*)fTracks->UncheckedAt(idx) : 0x0;
+ return idx >= 0 && idx < ntrk ? (AliKalmanTrack*)fTracks->UncheckedAt(idx) : 0x0;
}
//____________________________________________________________________
}
+
+//____________________________________________________________________
+void AliTRDtrackerV1::SetReconstructor(const AliTRDReconstructor *rec)
+{
+ fReconstructor = rec;
+ if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){
+ if(!fgDebugStreamer){
+ TDirectory *savedir = gDirectory;
+ fgDebugStreamer = new TTreeSRedirector("TRD.TrackerDebug.root");
+ savedir->cd();
+ }
+ }
+}
+
//_____________________________________________________________________________
Float_t AliTRDtrackerV1::GetChi2Y(AliTRDseedV1 *tracklets) const
{
Float_t chi2 = 0;
for(Int_t ipl = 0; ipl < kNPlanes; ipl++){
if(!tracklets[ipl].IsOK()) continue;
- Double_t distLayer = tracklets[ipl].GetYfit(0) - tracklets[ipl].GetYref(0);
+ 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; isl<kNSeedPlanes; isl++){
+ if(!fSeedTB[isl]) fSeedTB[isl] = new AliTRDchamberTimeBin();
+ else fSeedTB[isl]->Clear();
+ }
+}
+
//_____________________________________________________________________________
Float_t AliTRDtrackerV1::GetChi2Z(AliTRDseedV1 *tracklets) const
{
- // Chi2 definition on z-direction
+ // Calculates normalized chi2 in z-direction
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++){
if(!tracklets[ipl].IsOK()) continue;
- Double_t distLayer = tracklets[ipl].GetMeanz() - tracklets[ipl].GetZref(0);
+ Double_t distLayer = (tracklets[ipl].GetMeanz() - tracklets[ipl].GetZref(0)); // /(tracklets[ipl].GetPadLength()/2);
chi2 += distLayer * distLayer;
}
return chi2;
//
Double_t denominator = fSums[0] * fSums[4] - fSums[1] *fSums[1];
+ if(denominator==0) return;
+
// for(Int_t isum = 0; isum < 5; isum++)
// printf("fSums[%d] = %f\n", isum, fSums[isum]);
// printf("denominator = %f\n", denominator);