X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDtrackerV1.cxx;h=d32fee1b4d22bbbbfda616b5509fac86a5b4d04d;hb=77bbf113a85d8cf07857ee1ad2c0cb455f8bb7ae;hp=53e3188d005058b795d15832c76650f6f2bebc32;hpb=ae3fbe1fba81468406887c67b40e70b86e90bfbd;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDtrackerV1.cxx b/TRD/AliTRDtrackerV1.cxx index 53e3188d005..d32fee1b4d2 100644 --- a/TRD/AliTRDtrackerV1.cxx +++ b/TRD/AliTRDtrackerV1.cxx @@ -1,4 +1,3 @@ - /************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * @@ -38,6 +37,7 @@ #include #include "AliLog.h" +#include "AliMathBase.h" #include "AliESDEvent.h" #include "AliGeomManager.h" #include "AliRieman.h" @@ -74,14 +74,14 @@ Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = { 0.0474, 0.0408, 0.0335, 0.0335, 0.0335 }; Int_t AliTRDtrackerV1::fgNTimeBins = 0; -TTreeSRedirector *AliTRDtrackerV1::fgDebugStreamer = 0x0; AliRieman* AliTRDtrackerV1::fgRieman = 0x0; TLinearFitter* AliTRDtrackerV1::fgTiltedRieman = 0x0; TLinearFitter* AliTRDtrackerV1::fgTiltedRiemanConstrained = 0x0; //____________________________________________________________________ -AliTRDtrackerV1::AliTRDtrackerV1() +AliTRDtrackerV1::AliTRDtrackerV1(AliTRDReconstructor *rec) :AliTracker() + ,fReconstructor(0x0) ,fGeom(new AliTRDgeometry()) ,fClusters(0x0) ,fTracklets(0x0) @@ -99,24 +99,11 @@ AliTRDtrackerV1::AliTRDtrackerV1() if(!fgNTimeBins) fgNTimeBins = trd->GetNumberOfTimeBins(); for (Int_t isector = 0; isector < AliTRDgeometry::kNsector; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector); - - // retrive reco params - AliTRDrecoParam *rec = 0x0; - if (!(rec = AliTRDReconstructor::RecoParam())){ - if(!(rec = trd->GetRecoParam(0))){ - AliInfo("Using default RecoParams = LowFluxParam."); - rec = AliTRDrecoParam::GetLowFluxParam(); - } - AliTRDReconstructor::SetRecoParam(rec); - } for(Int_t isl =0; islGetStreamLevel() > 1){ - TDirectory *savedir = gDirectory; - fgDebugStreamer = new TTreeSRedirector("TRD.TrackerDebug.root"); - savedir->cd(); - } + // Initialize debug stream + if(rec) SetReconstructor(rec); } //____________________________________________________________________ @@ -126,14 +113,15 @@ AliTRDtrackerV1::~AliTRDtrackerV1() // Destructor // - if(fgDebugStreamer) delete fgDebugStreamer; - if(fgRieman) delete fgRieman; - if(fgTiltedRieman) delete fgTiltedRieman; - if(fgTiltedRiemanConstrained) delete fgTiltedRiemanConstrained; + if(fgRieman) delete fgRieman; fgRieman = 0x0; + if(fgTiltedRieman) delete fgTiltedRieman; fgTiltedRieman = 0x0; + if(fgTiltedRiemanConstrained) delete fgTiltedRiemanConstrained; fgTiltedRiemanConstrained = 0x0; for(Int_t isl =0; islDelete(); delete fTracks;} if(fTracklets) {fTracklets->Delete(); delete fTracklets;} - if(fClusters) {fClusters->Delete(); delete fClusters;} + if(fClusters) { + fClusters->Delete(); delete fClusters; + } if(fGeom) delete fGeom; } @@ -155,7 +143,7 @@ Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd) // See AliTRDtrackerV1::Clusters2TracksSM() for details. // - if(!AliTRDReconstructor::RecoParam()){ + if(!fReconstructor->GetRecoParam() ){ AliError("Reconstruction configuration not initialized. Call first AliTRDReconstructor::SetRecoParam()."); return 0; } @@ -177,13 +165,13 @@ Bool_t AliTRDtrackerV1::GetTrackPoint(Int_t index, AliTrackPoint &p) const { //AliInfo(Form("Asking for tracklet %d", index)); + // reset position of the point before using it + p.SetXYZ(0., 0., 0.); 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 idet = cl->GetDetector(); + Int_t idet = tracklet->GetDetector(); Double_t local[3]; local[0] = tracklet->GetX0(); @@ -240,7 +228,7 @@ 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; } @@ -276,7 +264,7 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) quality[iSeed] = covariance[0] + covariance[2]; } // Sort tracks according to covariance of local Y and Z - TMath::Sort(nSeed,quality,index,kFALSE); + TMath::Sort(nSeed, quality, index,kFALSE); } // Backpropagate all seeds @@ -294,11 +282,15 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) // Do the back prolongation new(&track) AliTRDtrackV1(*seed); - //track->Print(); + track.SetReconstructor(fReconstructor); + track.SetKink(Bool_t(seed->GetKinkIndex(0))); //Int_t lbl = seed->GetLabel(); //track.SetSeedLabel(lbl); - seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup); // Make backup - Float_t p4 = track.GetC(); + + // Make backup and mark entrance in the TRD + seed->UpdateTrackParams(&track, AliESDtrack::kTRDin); + seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup); + Float_t p4 = track.GetC(track.GetBz()); expectedClr = FollowBackProlongation(track); if (expectedClr<0) continue; // Back prolongation failed @@ -310,36 +302,29 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) // 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::RecoParam()->GetStreamLevel() > 0 /*&& quality TODO*/){ - AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(track); - calibTrack->SetOwner(); - seed->AddCalibObject(calibTrack); - } } } - if ((TMath::Abs(track.GetC() - p4) / TMath::Abs(p4) < 0.2) ||(track.Pt() > 0.8)) { - // + 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; @@ -376,13 +361,13 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) Double_t xtof = 371.0; Double_t xTOF0 = 370.0; - Double_t c2 = track.GetSnp() + track.GetC() * (xtof - track.GetX()); + 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() * (xtof - track.GetX()); + c2 = track.GetSnp() + track.GetC(track.GetBz()) * (xtof - track.GetX()); if (TMath::Abs(c2) >= 0.99) continue; //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) { @@ -419,11 +404,11 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) if(quality) delete [] quality; - AliInfo(Form("Number of seeds: %d", nSeed)); + AliInfo(Form("Number of TPC seeds: %d", nSeed)); AliInfo(Form("Number of back propagated TRD tracks: %d", found)); // run stand alone tracking - if (AliTRDReconstructor::RecoParam()->IsSeeding()) Clusters2Tracks(event); + if (fReconstructor->IsSeeding()) Clusters2Tracks(event); return 0; } @@ -454,9 +439,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(); - if((status & AliESDtrack::kTRDout) == 0) continue; - if((status & AliESDtrack::kTRDin) != 0) continue; + if(!(status & AliESDtrack::kTRDout)) continue; + if(!(status & AliESDtrack::kTRDin)) continue; nseed++; track.ResetCovariance(50.0); @@ -467,16 +454,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)); @@ -485,7 +482,6 @@ Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event) return 0; } - //____________________________________________________________________ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) { @@ -513,15 +509,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->GetX0(); + Double_t x = tracklet->GetX();//GetX0(); // reject tracklets which are not considered for inward refit if(x > t.GetX()+fgkMaxStep) continue; @@ -543,32 +539,43 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) xyz1[1] = x * TMath::Sin(alpha) + y * TMath::Cos(alpha); xyz1[2] = z; - // Get material budget - Double_t param[7]; - AliTracker::MeanMaterialBudget(xyz0, xyz1, param); - Double_t xrho= param[0]*param[4]; - Double_t xx0 = param[1]; // Get mean propagation parameters + Double_t length = TMath::Sqrt( + (xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0]) + + (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1]) + + (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2]) + ); + if(length>0.){ + // Get material budget + Double_t param[7]; + if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) break; + Double_t xrho= param[0]*param[4]; + Double_t xx0 = param[1]; // Get mean propagation parameters + + // Propagate and update + t.PropagateTo(x, xx0, xrho); + if (!AdjustSector(&t)) break; + } + if(kStoreIn){ + t.SetTrackHigh(); + kStoreIn = kFALSE; + } - // Propagate and update - t.PropagateTo(x, xx0, xrho); - if (!AdjustSector(&t)) break; - Double_t maxChi2 = t.GetPredictedChi2(tracklet); if (maxChi2 < 1e+10 && t.Update(tracklet, maxChi2)){ nClustersExpected += tracklet->GetN(); } } - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){ + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ Int_t index; - for(int iplane=0; iplane<6; iplane++){ + for(int iplane=0; iplaneGetDebugStream(AliTRDReconstructor::kTracker); cstreamer << "FollowProlongation" << "EventNumber=" << eventNumber << "ncl=" << nClustersExpected @@ -608,18 +615,29 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) // Int_t nClustersExpected = 0; - Double_t clength = AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick(); + Double_t clength = .5*AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick(); 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); + } + Bool_t kStoreIn = kTRUE; + // Loop through the TRD layers - for (Int_t ilayer = 0; ilayer < AliTRDgeometry::Nlayer(); ilayer++) { + for (Int_t ilayer = 0; ilayer < kNPlanes; ilayer++) { // BUILD TRACKLET IF NOT ALREADY BUILT - Double_t x = 0., y, z, alpha; - ptrTracklet = t.GetTracklet(ilayer); + Double_t x = 0., x0, y, z, alpha; + ptrTracklet = tracklets[ilayer]; if(!ptrTracklet){ ptrTracklet = new(&tracklet) AliTRDseedV1(ilayer); + ptrTracklet->SetReconstructor(fReconstructor); + ptrTracklet->SetKink(t.IsKink()); alpha = t.GetAlpha(); Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsector)); @@ -627,7 +645,13 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) if((x = fTrSec[sector].GetX(ilayer)) < 1.) continue; - if (!t.GetProlongation(x, y, z)) return -1; + // 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.; @@ -637,42 +661,49 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) 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(ilayer, stack); - tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle())); + tracklet.SetTilt(TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle())); tracklet.SetPadLength(pp->GetLengthIPad()); - tracklet.SetPlane(ilayer); + tracklet.SetPadWidth(pp->GetWidthIPad()); + tracklet.SetDetector(chamber->GetDetector()); tracklet.SetX0(x); - if(!tracklet.Init(&t)){ - t.SetStopped(kTRUE); - return nClustersExpected; - } - if(!tracklet.AttachClustersIter(chamber, 1000.)) continue; - tracklet.Init(&t); + 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 * AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue; + if(tracklet.GetN() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()) continue; break; } + ptrTracklet->UpdateUsed(); } if(!ptrTracklet->IsOK()){ + ptrTracklet->Reset(); if(x < 1.) continue; //temporary - if(!PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) return -1; - if(!AdjustSector(&t)) return -1; - if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1; + 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)) return -1; - if (!AdjustSector(&t)) return -1; - if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1; + 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->UseClusters(); + ptrTracklet->Fit(kFALSE); // no tilt correction ptrTracklet = SetTracklet(ptrTracklet); t.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1); @@ -682,23 +713,30 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) Double_t xyz0[3]; // entry point t.GetXYZ(xyz0); alpha = t.GetAlpha(); - x = ptrTracklet->GetX0(); - if (!t.GetProlongation(x, y, z)) return -1; + x = ptrTracklet->GetX(); //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]; - 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 - if (!t.PropagateTo(x, xx0, xrho)) return -1; - if (!AdjustSector(&t)) return -1; + if (!t.PropagateTo(x, xx0, xrho)) return -1/*nClustersExpected*/; + if (!AdjustSector(&t)) return -1/*nClustersExpected*/; + + if(kStoreIn){ + t.SetTrackLow(); + kStoreIn = kFALSE; + } Double_t maxChi2 = t.GetPredictedChi2(ptrTracklet); - if (!t.Update(ptrTracklet, maxChi2)) return -1; + if (!t.Update(ptrTracklet, maxChi2)) return -1/*nClustersExpected*/; + ptrTracklet->UpDate(&t); + if (maxChi2<1e+10) { nClustersExpected += ptrTracklet->GetN(); //t.SetTracklet(&tracklet, index); @@ -725,12 +763,13 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) //(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 - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 1){ - TTreeSRedirector &cstreamer = *fgDebugStreamer; + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ + TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); //AliTRDtrackV1 *debugTrack = new AliTRDtrackV1(t); //debugTrack->SetOwner(); @@ -768,7 +807,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 @@ -855,7 +894,7 @@ 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); x = cl->GetX(); @@ -867,7 +906,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()) * t; fitter->AddPoint(uvt, w, error); nPoints++; } @@ -881,9 +920,9 @@ 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(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 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); @@ -891,7 +930,7 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub Float_t chi2Z = CalculateChi2Z(tracklets, zref, slope, xref); Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); - TTreeSRedirector &treeStreamer = *fgDebugStreamer; + TTreeSRedirector &treeStreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); treeStreamer << "FitTiltedRiemanConstraint" << "EventNumber=" << eventNumber << "CandidateNumber=" << candidateNumber @@ -900,7 +939,7 @@ Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Doub << "Chi2Z=" << chi2Z << "zref=" << zref << "\n"; - } + }*/ return chi2track; } @@ -949,13 +988,13 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro // 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 (!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); @@ -966,7 +1005,7 @@ 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()) : 0.2; fitter->AddPoint(uvt, w, we); zfitter.AddPoint(&x, z, static_cast(TMath::Sqrt(cl->GetSigmaZ2()))); nPoints++; @@ -985,7 +1024,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) { @@ -1055,8 +1094,8 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro tracklets[iLayer].SetChi2(chi2track); } - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >=5){ - TTreeSRedirector &cstreamer = *fgDebugStreamer; +/* if(fReconstructor->GetStreamLevel() >=5){ + TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); Double_t chi2z = CalculateChi2Z(tracklets, offset, slope, xref); @@ -1066,13 +1105,13 @@ Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigErro << "xref=" << xref << "Chi2Z=" << chi2z << "\n"; - } + }*/ return chi2track; } //____________________________________________________________________ -Double_t AliTRDtrackerV1::FitLine(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t err, Int_t np, AliTrackPoint *points) +Double_t AliTRDtrackerV1::FitLine(const AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t err, Int_t np, AliTrackPoint *points) { AliTRDLeastSquare yfitter, zfitter; AliTRDcluster *cl = 0x0; @@ -1140,7 +1179,7 @@ Double_t AliTRDtrackerV1::FitLine(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, //_________________________________________________________________________ -Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t sigError, Int_t np, AliTrackPoint *points) +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 @@ -1194,7 +1233,7 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *trac // 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(); @@ -1211,13 +1250,14 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *trac 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++; } } - fitter->Eval(); + if(fitter->Eval()) return 1.E10; + Double_t z0 = fitter->GetParameter(3); Double_t dzdx = fitter->GetParameter(4); @@ -1229,7 +1269,7 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *trac 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) { @@ -1251,7 +1291,9 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *trac 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); @@ -1262,7 +1304,9 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *trac 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) @@ -1275,34 +1319,137 @@ Double_t AliTRDtrackerV1::FitRiemanTilt(AliTRDtrackV1 *track, AliTRDseedV1 *trac tracklets[ip].SetChi2(chi2); } } - //update track points array if(np && points){ Float_t xyz[3]; for(int ip=0; ip0.?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); } } - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >=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(ipGetX() - 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( AliTRDReconstructor::Tracker()))) continue; + ptrTracklet = new(&tracklet) AliTRDseedV1(iplane); + if(!tracker->MakeTracklet(ptrTracklet, track)) */ + continue; + } + } + if(!ptrTracklet->IsOK()) continue; + + Double_t x = ptrTracklet->GetX0(); + + while(ip < np){ + //don't do anything if next marker is after next update point. + if((up?-1:1) * (points[ip].GetX() - x) - fgkMaxStep < 0) break; + if(((up?-1:1) * (points[ip].GetX() - track->GetX()) < 0) && !PropagateToX(*track, points[ip].GetX(), fgkMaxStep)) return -1.; + + Double_t xyz[3]; // should also get the covariance + track->GetXYZ(xyz); + track->Global2LocalPosition(xyz, track->GetAlpha()); + points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]); + ip++; + } + // printf("plane[%d] tracklet[%p] x[%f]\n", iplane, ptrTracklet, x); + + // Propagate closer to the next update point + if(((up?-1:1) * (x - track->GetX()) + fgkMaxStep < 0) && !PropagateToX(*track, x + (up?-1:1)*fgkMaxStep, fgkMaxStep)) return -1.; + + if(!AdjustSector(track)) return -1; + if(TMath::Abs(track->GetSnp()) > fgkMaxSnp) return -1; + + //load tracklet to the tracker and the track +/* Int_t index; + if((index = FindTracklet(ptrTracklet)) < 0){ + ptrTracklet = SetTracklet(&tracklet); + index = fTracklets->GetEntriesFast()-1; + } + track->SetTracklet(ptrTracklet, index);*/ + + + // register tracklet to track with tracklet creation !! + // PropagateBack : loaded tracklet to the tracker and update index + // RefitInward : update index + // MakeTrack : loaded tracklet to the tracker and update index + if(!tracklets) track->SetTracklet(ptrTracklet, -1); + + + //Calculate the mean material budget along the path inside the chamber + Double_t xyz0[3]; track->GetXYZ(xyz0); + Double_t alpha = track->GetAlpha(); + Double_t xyz1[3], y, z; + if(!track->GetProlongation(x, y, z)) return -1; + xyz1[0] = x * TMath::Cos(alpha) - y * TMath::Sin(alpha); + xyz1[1] = +x * TMath::Sin(alpha) + y * TMath::Cos(alpha); + xyz1[2] = z; + if((xyz0[0] - xyz1[9] < 1e-3) && (xyz0[0] - xyz1[9] < 1e-3)) continue; // check wheter we are at the same global x position + Double_t param[7]; + if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param) <=0.) break; + Double_t xrho = param[0]*param[4]; // density*length + Double_t xx0 = param[1]; // radiation length + + //Propagate the track + track->PropagateTo(x, xx0, xrho); + if (!AdjustSector(track)) break; + + //Update track + Double_t chi2 = track->GetPredictedChi2(ptrTracklet); + if(chi2<1e+10) track->Update(ptrTracklet, chi2); + if(!up) continue; + + //Reset material budget if 2 consecutive gold + if(iplane>0 && track->GetTracklet(iplane-1) && ptrTracklet->GetN() + track->GetTracklet(iplane-1)->GetN() > 20) track->SetBudget(2, 0.); + } // end planes loop + + // extrapolation + while(ip < np){ + if(((up?-1:1) * (points[ip].GetX() - track->GetX()) < 0) && !PropagateToX(*track, points[ip].GetX(), fgkMaxStep)) return -1.; + + Double_t xyz[3]; // should also get the covariance + track->GetXYZ(xyz); + track->Global2LocalPosition(xyz, track->GetAlpha()); + points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]); + ip++; + } + + return track->GetChi2(); +} //_________________________________________________________________________ Float_t AliTRDtrackerV1::CalculateChi2Z(AliTRDseedV1 *tracklets, Double_t offset, Double_t slope, Double_t xref) @@ -1321,7 +1468,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); @@ -1375,12 +1522,10 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m // 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])) { - return 0; - } + if (!t.PropagateTo(x, param[1], param[0]*param[4])) return 0; // Rotate the track if necessary AdjustSector(&t); @@ -1415,7 +1560,9 @@ Int_t AliTRDtrackerV1::ReadClusters(TClonesArray* &array, TTree *clusterTree) co 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); } @@ -1447,23 +1594,64 @@ Int_t AliTRDtrackerV1::ReadClusters(TClonesArray* &array, TTree *clusterTree) co 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++; @@ -1474,17 +1662,18 @@ Int_t AliTRDtrackerV1::LoadClusters(TTree *cTree) fTrSec[sector].GetChamber(stack, layer, kTRUE)->InsertCluster(c, icl); } - AliInfo(Form("Clusters %d in %6.2f %%", ncl, 100.*float(nin)/ncl)); - + + const AliTRDCalDet *cal = AliTRDcalibDB::Instance()->GetT0Det(); for(int isector =0; isectorDelete(); 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::kNsector; i++) fTrSec[i].Clear(); @@ -1502,6 +1699,24 @@ 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::kNclusters; ic--;){ +// if(!(c=tracklet->GetClusters(ic))) continue; +// c->Use(); +// } +// } +// } +// + //_____________________________________________________________________________ Bool_t AliTRDtrackerV1::AdjustSector(AliTRDtrackV1 *track) { @@ -1512,7 +1727,7 @@ Bool_t AliTRDtrackerV1::AdjustSector(AliTRDtrackV1 *track) Double_t alpha = AliTRDgeometry::GetAlpha(); Double_t y = track->GetY(); Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha); - + if (y > ymax) { if (!track->Rotate( alpha)) { return kFALSE; @@ -1631,7 +1846,7 @@ Int_t AliTRDtrackerV1::Clusters2TracksSM(Int_t sector, AliESDEvent *esd) nChambers = 0; for(int ilayer=0; ilayerGetNClusters() < fgNTimeBins * AliTRDReconstructor::RecoParam()->GetFindableClusters()) continue; + if(chamber->GetNClusters() < fgNTimeBins * fReconstructor->GetRecoParam() ->GetFindableClusters()) continue; nChambers++; //AliInfo(Form("sector %d stack %d layer %d clusters %d", sector, istack, ilayer, chamber->GetNClusters())); } @@ -1679,7 +1894,9 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon // 8. Build ESD track and register it to the output list // + 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 @@ -1687,28 +1904,44 @@ 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(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 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; ci = &stack[0]; + while(icGetStack((*ci)->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::RecoParam()->GetStreamLevel() > 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; @@ -1736,209 +1969,146 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon // 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++; - } - } - } - // 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; - - - // 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); - + // 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.158) findable++; - // Sign clusters - AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1; - for (Int_t jLayer = 0; jLayer < 6; 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; - while(!(cl = sseed[jseed].GetClusters(ic))) ic++; - clusterIndex = sseed[jseed].GetIndexes(ic); - } - } - if(!cl) continue; + sseed[jseed].UpdateUsed(); + ncl += sseed[jseed].GetN2(); + nused += sseed[jseed].GetNUsed(); + nlayers++; + } - - // Build track parameters - AliTRDseedV1 *lseed =&sseed[trackIndex*6]; - Int_t idx = 0; - while(idx<3 && !lseed->IsOK()) { - idx++; - lseed++; - } - Double_t cR = lseed->GetC(); - 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[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1)); - trackParams[5] = cR; - Int_t ich = 0; while(!(chamber = stack[ich])) ich++; - trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *alpha+shift; // Supermodule*/ - - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 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(); - } - //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 = *fgDebugStreamer; - 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"; - } + // 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; + + // 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(); + //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] + << "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; - esdTrack.UpdateTrackParams(track, AliESDtrack::kTRDout); - esdTrack.SetLabel(track->GetLabel()); - track->UpdateESDtrack(&esdTrack); - // write ESD-friends if neccessary - if (AliTRDReconstructor::RecoParam()->GetStreamLevel() > 0){ - AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(*track); - calibTrack->SetOwner(); - esdTrack.AddCalibObject(calibTrack); - } - new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack); - 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++; @@ -1947,18 +2117,20 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClon // 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::RecoParam()->GetStreamLevel() > 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 @@ -2073,7 +2245,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss // 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]; @@ -2082,6 +2254,9 @@ 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 Double_t seedQuality[kMaxTracksStack]; @@ -2089,32 +2264,36 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss // 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->GetStack(chamber->GetDetector()); 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()); + hL[iplane] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle()); padlength[iplane] = pp->GetLengthIPad(); + padwidth[iplane] = pp->GetWidthIPad(); } - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() > 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; for(int isl=0; islGetSeedingLayer(fSeedTB[isl], fGeom)) continue; + if(!chamber->GetSeedingLayer(fSeedTB[isl], fGeom, fReconstructor)) continue; nlayers++; } - if(nlayers < 4) return 0; + if(nlayers < 4) return ntracks; // Start finding seeds @@ -2148,23 +2327,25 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss // 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(); + for (Int_t il = 0; il < kNPlanes; il++) cseed[il].Reset(); FitRieman(c, chi2); - AliTRDseedV1 *tseed = 0x0; - for(int iLayer=0; iLayerSetPlane(jLayer); - tseed->SetTilt(hL[jLayer]); - tseed->SetPadLength(padlength[jLayer]); - tseed->SetX0(stack[jLayer]->GetX()); + AliTRDseedV1 *tseed = &cseed[0]; + AliTRDtrackingChamber **cIter = &stack[0]; + for(int iLayer=0; iLayerSetDetector((*cIter) ? (*cIter)->GetDetector() : -1); + tseed->SetTilt(hL[iLayer]); + tseed->SetPadLength(padlength[iLayer]); + tseed->SetPadWidth(padwidth[iLayer]); + tseed->SetReconstructor(fReconstructor); + tseed->SetX0((*cIter) ? (*cIter)->GetX() : x_def[iLayer]); tseed->Init(GetRiemanFitter()); + tseed->SetStandAlone(kTRUE); } Bool_t isFake = kFALSE; - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){ + 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; @@ -2177,7 +2358,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); AliRieman *rim = GetRiemanFitter(); - TTreeSRedirector &cs0 = *fgDebugStreamer; + TTreeSRedirector &cs0 = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); cs0 << "MakeSeeds0" <<"EventNumber=" << eventNumber <<"CandidateNumber=" << candidateNumber @@ -2207,44 +2388,60 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss <<"RiemanFitter.=" << rim <<"\n"; } - - if(chi2[0] > AliTRDReconstructor::RecoParam()->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){ - //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0])); + 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] > AliTRDReconstructor::RecoParam()->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){ - //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1])); + 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 25) break; + if(!cseed[jLayer].AttachClusters(stack[jLayer], kTRUE)) continue; + cseed[jLayer].UpdateUsed(); + if(!cseed[jLayer].IsOK()) continue; mlayers++; } + if(mlayers < kNSeedPlanes){ //AliInfo(Form("Failed updating all seeds %d [%d].", mlayers, kNSeedPlanes)); AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1); continue; } + + // temporary exit door for the HLT + if(fReconstructor->IsHLT()){ + // attach clusters to extrapolation chambers + for(int iLayer=0; iLayerGetTrackLikelihood()){ + 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; @@ -2256,34 +2453,25 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss fSeedLayer[ntracks] = config;/*sLayer;*/ // 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]; + for(int iLayer=0; iLayerGetX()); - 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; + if(!pseed.AttachClusters(chamber, kTRUE)) continue; + pseed.Fit(kTRUE); cseed[jLayer] = pseed; - nusedf += cseed[jLayer].GetNUsed(); // debug value FitTiltedRieman(cseed, kTRUE); + cseed[jLayer].Fit(kTRUE); } // AliInfo("Extrapolation done."); // Debug Stream containing all the 6 tracklets - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){ - TTreeSRedirector &cstreamer = *fgDebugStreamer; + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){ + TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); TLinearFitter *tiltedRieman = GetTiltedRiemanFitter(); Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); @@ -2300,7 +2488,7 @@ Int_t AliTRDtrackerV1::MakeSeeds(AliTRDtrackingChamber **stack, AliTRDseedV1 *ss << "\n"; } - if(ImproveSeedQuality(stack, cseed) < 4){ + if(fReconstructor->GetRecoParam()->HasImproveTracklets() && ImproveSeedQuality(stack, cseed) < 4){ AliTRDtrackerDebug::SetCandidateNumber(AliTRDtrackerDebug::GetCandidateNumber() + 1); continue; } @@ -2314,7 +2502,7 @@ 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(AliTRDReconstructor::RecoParam()->IsVertexConstrained()) + if(fReconstructor->GetRecoParam()->IsVertexConstrained()) chi2Vals[1] = FitTiltedRiemanConstraint(&cseed[0], GetZ()); // Do Vertex Constrained fit if desired else chi2Vals[1] = 1.; @@ -2323,46 +2511,24 @@ 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(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){ - TTreeSRedirector &cstreamer = *fgDebugStreamer; + + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){ + TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); 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 - << "NUsedS=" << nUsedCl - << "NUsed=" << nusedf + << "NClusters=" << ncls << "Like=" << like << "S0.=" << &cseed[0] << "S1.=" << &cseed[1] @@ -2370,8 +2536,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"; @@ -2409,8 +2573,6 @@ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params) // 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; @@ -2424,17 +2586,66 @@ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params) AliTRDtrackV1 track(seeds, ¶ms[1], c, params[0], params[6]*alpha+shift); track.PropagateTo(params[0]-5.0); + AliTRDseedV1 *ptrTracklet = 0x0; + // Sign clusters + for (Int_t jLayer = 0; jLayer < AliTRDgeometry::kNlayer; jLayer++) { + ptrTracklet = &seeds[jLayer]; + if(!ptrTracklet->IsOK()) continue; + if(TMath::Abs(ptrTracklet->GetYref(1) - ptrTracklet->GetYfit(1)) >= .2) continue; // check this condition with Marian + } + // + if(fReconstructor->IsHLT()){ + for(Int_t ip=0; ipUseClusters(); + track.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1); + } + AliTRDtrackV1 *ptrTrack = SetTrack(&track); + ptrTrack->SetReconstructor(fReconstructor); + return ptrTrack; + } + track.ResetCovariance(1); - Int_t nc = FollowBackProlongation(track); + 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 = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); + 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; AliTRDtrackV1 *ptrTrack = SetTrack(&track); + ptrTrack->SetReconstructor(fReconstructor); ptrTrack->CookLabel(.9); + // computes PID for track ptrTrack->CookPID(); // update calibration references using this track - if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(ptrTrack); - + AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance(); + if (!calibra){ + AliInfo("Could not get Calibra instance\n"); + if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(ptrTrack); + } return ptrTrack; } @@ -2449,7 +2660,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 @@ -2477,7 +2688,7 @@ Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDs 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; @@ -2491,16 +2702,17 @@ 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++; } chi2 = FitTiltedRieman(bseed, kTRUE); - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 7){ + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 7){ Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); TLinearFitter *tiltedRieman = GetTiltedRiemanFitter(); - TTreeSRedirector &cstreamer = *fgDebugStreamer; + TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); cstreamer << "ImproveSeedQuality" << "EventNumber=" << eventNumber << "CandidateNumber=" << candidateNumber @@ -2515,7 +2727,6 @@ Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDs << "\n"; } } // Loop: iter - // we are sure that at least 2 tracklets are OK ! return nLayers+2; } @@ -2540,32 +2751,32 @@ 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 = (AliTRDReconstructor::RecoParam()->IsVertexConstrained()) ? + 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; + Double_t likeChi2Phi= TMath::Exp(-chi2phi * 3.23); + Double_t trackLikelihood = likeChi2Z * likeChi2TR * likeChi2Phi; - if(AliTRDReconstructor::RecoParam()->GetStreamLevel() >= 2){ + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){ Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); - TTreeSRedirector &cstreamer = *fgDebugStreamer; + TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); cstreamer << "CalculateTrackLikelihood0" << "EventNumber=" << eventNumber << "CandidateNumber=" << candidateNumber << "LikeChi2Z=" << likeChi2Z << "LikeChi2TR=" << likeChi2TR << "LikeChi2TC=" << likeChi2TC - << "LikeAF=" << likeAF + << "LikeChi2Phi=" << likeChi2Phi << "TrackLikelihood=" << trackLikelihood << "\n"; } @@ -2574,8 +2785,7 @@ Double_t AliTRDtrackerV1::CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Doub } //____________________________________________________________________ -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. @@ -2602,31 +2812,40 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4] // // 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)); + sumda += TMath::Abs(cseed[jlayer].GetYfit(1) - cseed[jlayer].GetYref(1)); } - Double_t likea = TMath::Exp(-sumda*10.6); + nclusters *= .25; + + Double_t likea = TMath::Exp(-sumda * fRecoPars->GetPhiSlope()); 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 < fRecoPars->GetChi2YCut()) likechi2y += TMath::Exp(-TMath::Sqrt(chi2y) * fRecoPars->GetChi2YSlope()); + Double_t likechi2z = TMath::Exp(-chi2z * fRecoPars->GetChi2ZSlope()); + 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::RecoParam()->GetStreamLevel() >= 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; + TTreeSRedirector &cstreamer = *fReconstructor->GetDebugStream(AliTRDReconstructor::kTracker); cstreamer << "CookLikelihood" << "EventNumber=" << eventNumber << "CandidateNumber=" << candidateNumber @@ -2637,22 +2856,21 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed, Int_t planes[4] << "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"; } return like; } - - //____________________________________________________________________ void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4]) { @@ -2948,78 +3166,65 @@ Float_t AliTRDtrackerV1::CalculateReferenceX(AliTRDseedV1 *tracklets){ 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; - - // 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; - -} - -//_____________________________________________________________________________ -Float_t AliTRDtrackerV1::GetChi2Y(AliTRDseedV1 *tracklets) const -{ - // Chi2 definition on y-direction +// //_____________________________________________________________________________ +// 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; +// +// } - 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); - chi2 += distLayer * distLayer; - } - return chi2; -} //____________________________________________________________________ void AliTRDtrackerV1::ResetSeedTB() @@ -3033,18 +3238,35 @@ void AliTRDtrackerV1::ResetSeedTB() } } + +//_____________________________________________________________________________ +Float_t AliTRDtrackerV1::GetChi2Y(AliTRDseedV1 *tracklets) const +{ + // Calculates normalized chi2 in y-direction + // chi2 = Sum chi2 / n_tracklets + + Double_t chi2 = 0.; Int_t n = 0; + for(Int_t ipl = kNPlanes; ipl--;){ + if(!tracklets[ipl].IsOK()) continue; + chi2 += tracklets[ipl].GetChi2Y(); + n++; + } + return n ? chi2/n : 0.; +} + //_____________________________________________________________________________ Float_t AliTRDtrackerV1::GetChi2Z(AliTRDseedV1 *tracklets) const { - // Chi2 definition on z-direction + // Calculates normalized chi2 in z-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].GetMeanz() - tracklets[ipl].GetZref(0); - chi2 += distLayer * distLayer; + chi2 += tracklets[ipl].GetChi2Z(); + n++; } - return chi2; + return n ? chi2/n : 0.; } /////////////////////////////////////////////////////// @@ -3104,6 +3326,8 @@ void AliTRDtrackerV1::AliTRDLeastSquare::Eval(){ // 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);