X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDtrackerV1.cxx;h=1b6b50367b47f958be776bd42674636db7b8e18f;hb=83dea92e6165d3705d3111d83924ef38692f994e;hp=26b4f0c4267a26c368808dee30586bf008a25e75;hpb=bccda319e1502ef207681dcfcb2999522f58c5c3;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDtrackerV1.cxx b/TRD/AliTRDtrackerV1.cxx index 26b4f0c4267..1b6b50367b4 100644 --- a/TRD/AliTRDtrackerV1.cxx +++ b/TRD/AliTRDtrackerV1.cxx @@ -1,18 +1,18 @@ /************************************************************************** - * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * * - * Author: The ALICE Off-line Project. * - * Contributors are mentioned in the code where appropriate. * - * * - * Permission to use, copy, modify and distribute this software and its * - * documentation strictly for non-commercial purposes is hereby granted * - * without fee, provided that the above copyright notice appears in all * - * copies and that both the copyright notice and this permission notice * - * appear in the supporting documentation. The authors make no claims * - * about the suitability of this software for any purpose. It is * - * provided "as is" without express or implied warranty. * - **************************************************************************/ +* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * +* * +* Author: The ALICE Off-line Project. * +* Contributors are mentioned in the code where appropriate. * +* * +* Permission to use, copy, modify and distribute this software and its * +* documentation strictly for non-commercial purposes is hereby granted * +* without fee, provided that the above copyright notice appears in all * +* copies and that both the copyright notice and this permission notice * +* appear in the supporting documentation. The authors make no claims * +* about the suitability of this software for any purpose. It is * +* provided "as is" without express or implied warranty. * +**************************************************************************/ /* $Id$ */ @@ -26,103 +26,105 @@ // // /////////////////////////////////////////////////////////////////////////////// -#include -#include -#include +// #include +// #include +// #include #include -#include -#include -#include -#include +#include #include -#include -#include #include #include -#include #include #include "AliLog.h" #include "AliESDEvent.h" -#include "AliAlignObj.h" +#include "AliGeomManager.h" #include "AliRieman.h" #include "AliTrackPointArray.h" -#include "AliTRDtracker.h" -#include "AliTRDtrackerV1.h" #include "AliTRDgeometry.h" #include "AliTRDpadPlane.h" -#include "AliTRDgeometry.h" -#include "AliTRDcluster.h" -#include "AliTRDtrack.h" -#include "AliTRDseed.h" #include "AliTRDcalibDB.h" -#include "AliTRDCommonParam.h" #include "AliTRDReconstructor.h" #include "AliTRDCalibraFillHisto.h" -#include "AliTRDtrackerFitter.h" -#include "AliTRDstackLayer.h" #include "AliTRDrecoParam.h" + +#include "AliTRDcluster.h" #include "AliTRDseedV1.h" #include "AliTRDtrackV1.h" -#include "Cal/AliTRDCalDet.h" +#include "AliTRDtrackerV1.h" +#include "AliTRDtrackerDebug.h" +#include "AliTRDtrackingChamber.h" +#include "AliTRDchamberTimeBin.h" + -#define DEBUG ClassImp(AliTRDtrackerV1) + + +const Float_t AliTRDtrackerV1::fgkMinClustersInTrack = 0.5; // +const Float_t AliTRDtrackerV1::fgkLabelFraction = 0.8; // +const Double_t AliTRDtrackerV1::fgkMaxChi2 = 12.0; // +const Double_t AliTRDtrackerV1::fgkMaxSnp = 0.95; // Maximum local sine of the azimuthal angle +const Double_t AliTRDtrackerV1::fgkMaxStep = 2.0; // Maximal step size in propagation Double_t AliTRDtrackerV1::fgTopologicQA[kNConfigs] = { - 0.1112, 0.1112, 0.1112, 0.0786, 0.0786, - 0.0786, 0.0786, 0.0579, 0.0579, 0.0474, - 0.0474, 0.0408, 0.0335, 0.0335, 0.0335 + 0.1112, 0.1112, 0.1112, 0.0786, 0.0786, + 0.0786, 0.0786, 0.0579, 0.0579, 0.0474, + 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(AliTRDrecoParam *p) - :AliTRDtracker() - ,fSieveSeeding(0) +AliTRDtrackerV1::AliTRDtrackerV1(AliTRDReconstructor *rec) + :AliTracker() + ,fReconstructor(0x0) + ,fGeom(new AliTRDgeometry()) + ,fClusters(0x0) ,fTracklets(0x0) - ,fRecoParam(p) - ,fFitter(0x0) -{ - // - // Default constructor. Nothing is initialized. - // - -} - -//____________________________________________________________________ -AliTRDtrackerV1::AliTRDtrackerV1(const TFile *in, AliTRDrecoParam *p) - :AliTRDtracker(in) + ,fTracks(0x0) ,fSieveSeeding(0) - ,fTracklets(0x0) - ,fRecoParam(p) - ,fFitter(0x0) { // - // Standard constructor. - // Setting of the geometry file, debug output (if enabled) - // and initilize fitter helper. - // + // Default constructor. + // + AliTRDcalibDB *trd = 0x0; + if (!(trd = AliTRDcalibDB::Instance())) { + AliFatal("Could not get calibration object"); + } - fFitter = new AliTRDtrackerFitter(); + if(!fgNTimeBins) fgNTimeBins = trd->GetNumberOfTimeBins(); -#ifdef DEBUG - fFitter->SetDebugStream(fDebugStreamer); -#endif + for (Int_t isector = 0; isector < AliTRDgeometry::kNsector; isector++) new(&fTrSec[isector]) AliTRDtrackingSector(fGeom, isector); + + for(Int_t isl =0; islDelete(); delete fTracklets;} + + if(fgDebugStreamer) delete fgDebugStreamer; + if(fgRieman) delete fgRieman; + if(fgTiltedRieman) delete fgTiltedRieman; + if(fgTiltedRiemanConstrained) delete fgTiltedRiemanConstrained; + for(Int_t isl =0; islDelete(); delete fTracks;} + if(fTracklets) {fTracklets->Delete(); delete fTracklets;} + if(fClusters) { + fClusters->Delete(); delete fClusters; + } + if(fGeom) delete fGeom; } //____________________________________________________________________ @@ -143,34 +145,95 @@ Int_t AliTRDtrackerV1::Clusters2Tracks(AliESDEvent *esd) // See AliTRDtrackerV1::Clusters2TracksSM() for details. // - if(!fRecoParam){ - AliError("Reconstruction configuration not initialized. Call first AliTRDtrackerV1::SetRecoParam()."); - return 0; - } - - //AliInfo("Start Track Finder ..."); - Int_t ntracks = 0; - for(int ism=0; ismGetRecoParam() ){ + AliError("Reconstruction configuration not initialized. Call first AliTRDReconstructor::SetRecoParam()."); + return 0; + } + + //AliInfo("Start Track Finder ..."); + Int_t ntracks = 0; + for(int ism=0; ismUncheckedAt(index); - // etc - return kTRUE; + //AliInfo(Form("Asking for tracklet %d", index)); + + 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(); + + Double_t local[3]; + local[0] = tracklet->GetX0(); + local[1] = tracklet->GetYfit(0); + local[2] = tracklet->GetZfit(0); + Double_t global[3]; + fGeom->RotateBack(idet, local, global); + p.SetXYZ(global[0],global[1],global[2]); + + + // setting volume id + AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1; + switch (fGeom->GetLayer(idet)) { + case 0: + iLayer = AliGeomManager::kTRD1; + break; + case 1: + iLayer = AliGeomManager::kTRD2; + break; + case 2: + iLayer = AliGeomManager::kTRD3; + break; + case 3: + iLayer = AliGeomManager::kTRD4; + break; + case 4: + iLayer = AliGeomManager::kTRD5; + break; + case 5: + iLayer = AliGeomManager::kTRD6; + break; + }; + Int_t modId = fGeom->GetSector(idet) * fGeom->Nstack() + fGeom->GetStack(idet); + UShort_t volid = AliGeomManager::LayerToVolUID(iLayer, modId); + p.SetVolumeID(volid); + + return kTRUE; } +//____________________________________________________________________ +TLinearFitter* AliTRDtrackerV1::GetTiltedRiemanFitter() +{ + if(!fgTiltedRieman) fgTiltedRieman = new TLinearFitter(4, "hyp4"); + return fgTiltedRieman; +} +//____________________________________________________________________ +TLinearFitter* AliTRDtrackerV1::GetTiltedRiemanFitterConstraint() +{ + if(!fgTiltedRiemanConstrained) fgTiltedRiemanConstrained = new TLinearFitter(2, "hyp2"); + return fgTiltedRiemanConstrained; +} + +//____________________________________________________________________ +AliRieman* AliTRDtrackerV1::GetRiemanFitter() +{ + if(!fgRieman) fgRieman = new AliRieman(AliTRDtrackingChamber::kNTimeBins * AliTRDgeometry::kNlayer); + return fgRieman; +} + //_____________________________________________________________________________ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) { @@ -183,242 +246,178 @@ Int_t AliTRDtrackerV1::PropagateBack(AliESDEvent *event) // by the TPC tracker. // - Int_t found = 0; // number of tracks found - Float_t foundMin = 20.0; - - AliTRDseed::SetNTimeBins(fTimeBinsPerPlane); - 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]; - } - // Sort tracks according to covariance of local Y and Z - TMath::Sort(nSeed,quality,index,kFALSE); - - // Backpropagate all seeds - for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) { - - // Get the seeds in sorted sequence - AliESDtrack *seed = event->GetTrack(index[iSeed]); - - // Check the seed status - ULong_t status = seed->GetStatus(); - if ((status & AliESDtrack::kTPCout) == 0) continue; - if ((status & AliESDtrack::kTRDout) != 0) continue; - - // Do the back prolongation - Int_t lbl = seed->GetLabel(); - AliTRDtrackV1 *track = new AliTRDtrackV1(*seed); - //track->Print(); - track->SetSeedLabel(lbl); - seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); // Make backup - fNseeds++; - Float_t p4 = track->GetC(); - Int_t expectedClr = FollowBackProlongation(*track); - //AliInfo(Form("\nTRACK %d Clusters %d [%d] in chi2 %f", index[iSeed], expectedClr, track->GetNumberOfClusters(), track->GetChi2())); - //track->Print(); - - //Double_t cov[15]; - //seed->GetExternalCovariance(cov); - //AliInfo(Form("track %d cov[%f %f] 0", index[iSeed], cov[0], cov[2])); - - if ((TMath::Abs(track->GetC() - 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()); // A.Bercuci 25.07.07 - CookLabel(track,1 - fgkLabelFraction); - if (track->GetBackupTrack()) UseClusters(track->GetBackupTrack()); - - - //seed->GetExternalCovariance(cov); - //AliInfo(Form("track %d cov[%f %f] 0 test", index[iSeed], cov[0], cov[2])); - - // Sign only gold tracks - if (track->GetChi2() / track->GetNumberOfClusters() < 4) { - if ((seed->GetKinkIndex(0) == 0) && - (track->Pt() < 1.5)) UseClusters(track); - } - Bool_t isGold = kFALSE; - - // Full gold track - if (track->GetChi2() / track->GetNumberOfClusters() < 5) { - if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup); - - isGold = kTRUE; - } - //seed->GetExternalCovariance(cov); - //AliInfo(Form("track %d cov[%f %f] 00", index[iSeed], cov[0], cov[2])); - - // Almost gold track - if ((!isGold) && (track->GetNCross() == 0) && - (track->GetChi2() / track->GetNumberOfClusters() < 7)) { - //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); - if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup); - - isGold = kTRUE; - } - //seed->GetExternalCovariance(cov); - //AliInfo(Form("track %d cov[%f %f] 01", index[iSeed], cov[0], cov[2])); - - if ((!isGold) && (track->GetBackupTrack())) { - if ((track->GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track->GetBackupTrack()->GetChi2()/(track->GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) { - seed->UpdateTrackParams(track->GetBackupTrack(),AliESDtrack::kTRDbackup); - isGold = kTRUE; - } - } - //seed->GetExternalCovariance(cov); - //AliInfo(Form("track %d cov[%f %f] 02", index[iSeed], cov[0], cov[2])); - - //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) { - //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup); - //} - } - } - /**/ - - /**/ - // Debug part of tracking -/* TTreeSRedirector &cstream = *fDebugStreamer; - Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number. - if (AliTRDReconstructor::StreamLevel() > 0) { - if (track->GetBackupTrack()) { - cstream << "Tracks" - << "EventNrInFile=" << eventNrInFile - << "ESD.=" << seed - << "trd.=" << track - << "trdback.=" << track->GetBackupTrack() - << "\n"; - } - else { - cstream << "Tracks" - << "EventNrInFile=" << eventNrInFile - << "ESD.=" << seed - << "trd.=" << track - << "trdback.=" << track - << "\n"; - } - }*/ - /**/ - - //seed->GetExternalCovariance(cov); - //AliInfo(Form("track %d cov[%f %f] 1", index[iSeed], cov[0], cov[2])); - - // Propagation to the TOF (I.Belikov) - if (track->GetStop() == kFALSE) { - //AliInfo("Track not stopped in TRD ..."); - Double_t xtof = 371.0; - Double_t xTOF0 = 370.0; - - Double_t c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX()); - if (TMath::Abs(c2) >= 0.99) { - delete track; - continue; - } - - PropagateToX(*track,xTOF0,fgkMaxStep); - - // Energy losses taken to the account - check one more time - c2 = track->GetSnp() + track->GetC() * (xtof - track->GetX()); - if (TMath::Abs(c2) >= 0.99) { - delete track; - continue; - } - - //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) { - // fHBackfit->Fill(7); - //delete track; - // continue; - //} - - Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha()); - Double_t y; - track->GetYAt(xtof,GetBz(),y); - if (y > ymax) { - if (!track->Rotate( AliTRDgeometry::GetAlpha())) { - delete track; - continue; - } - }else if (y < -ymax) { - if (!track->Rotate(-AliTRDgeometry::GetAlpha())) { - delete track; - continue; - } - } - - if (track->PropagateTo(xtof)) { - //AliInfo("set kTRDout"); - seed->UpdateTrackParams(track,AliESDtrack::kTRDout); - - for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) { - for (Int_t j = 0; j < AliESDtrack::kNSlice; j++) { - seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j); - } - seed->SetTRDTimBin(track->GetPIDTimBin(i),i); - } - //seed->SetTRDtrack(new AliTRDtrack(*track)); - if (track->GetNumberOfClusters() > foundMin) found++; - } - } else { - //AliInfo("Track stopped in TRD ..."); - - if ((track->GetNumberOfClusters() > 15) && - (track->GetNumberOfClusters() > 0.5*expectedClr)) { - seed->UpdateTrackParams(track,AliESDtrack::kTRDout); - - //seed->SetStatus(AliESDtrack::kTRDStop); - for (Int_t i = 0; i < AliESDtrack::kNPlane; i++) { - for (Int_t j = 0; j SetTRDsignals(track->GetPIDsignals(i,j),i,j); - } - seed->SetTRDTimBin(track->GetPIDTimBin(i),i); - } - //seed->SetTRDtrack(new AliTRDtrack(*track)); - found++; - } - } - - //if (((t->GetStatus()&AliESDtrack::kTRDout)!=0 ) - - seed->SetTRDQuality(track->StatusForTOF()); - seed->SetTRDBudget(track->GetBudget(0)); - -// if ((seed->GetStatus()&AliESDtrack::kTRDin)!=0 ) printf("TRDin "); -// if ((seed->GetStatus()&AliESDtrack::kTRDbackup)!=0 ) printf("TRDbackup "); -// if ((seed->GetStatus()&AliESDtrack::kTRDStop)!=0 ) printf("TRDstop "); -// if ((seed->GetStatus()&AliESDtrack::kTRDout)!=0 ) printf("TRDout "); -// printf("\n"); - delete track; - - //seed->GetExternalCovariance(cov); - //AliInfo(Form("track %d cov[%f %f] 2", index[iSeed], cov[0], cov[2])); - } - + // Calibration monitor + AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance(); + if (!calibra) AliInfo("Could not get Calibra instance\n"); + + 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){ + 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); + } + + // Backpropagate all seeds + Int_t expectedClr; + AliTRDtrackV1 track; + for (Int_t iSeed = 0; iSeed < nSeed; iSeed++) { + + // Get the seeds in sorted sequence + AliESDtrack *seed = event->GetTrack(index[iSeed]); + + // Check the seed status + ULong_t status = seed->GetStatus(); + if ((status & AliESDtrack::kTPCout) == 0) continue; + if ((status & AliESDtrack::kTRDout) != 0) continue; + + // Do the back prolongation + new(&track) AliTRDtrackV1(*seed); + track.SetReconstructor(fReconstructor); + + //Int_t lbl = seed->GetLabel(); + //track.SetSeedLabel(lbl); + + // Make backup and mark entrance in the TRD + seed->UpdateTrackParams(&track, AliESDtrack::kTRDin); + seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup); + Float_t p4 = track.GetC(); + 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); + } + } - AliInfo(Form("Number of seeds: %d",fNseeds)); - AliInfo(Form("Number of back propagated TRD tracks: %d",found)); - - //fSeeds->Clear(); - fNseeds = 0; - - delete [] index; - delete [] quality; - + if ((TMath::Abs(track.GetC() - 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()); + + // Sign only gold tracks + if (track.GetChi2() / track.GetNumberOfClusters() < 4) { + if ((seed->GetKinkIndex(0) == 0) && (track.Pt() < 1.5)){ + //UseClusters(&track); + } + } + Bool_t isGold = kFALSE; + + // Full gold track + if (track.GetChi2() / track.GetNumberOfClusters() < 5) { + if (track.GetBackupTrack()) seed->UpdateTrackParams(track.GetBackupTrack(),AliESDtrack::kTRDbackup); + + isGold = kTRUE; + } + + // Almost gold track + if ((!isGold) && (track.GetNCross() == 0) && (track.GetChi2() / track.GetNumberOfClusters() < 7)) { + //seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); + if (track.GetBackupTrack()) seed->UpdateTrackParams(track.GetBackupTrack(),AliESDtrack::kTRDbackup); + + isGold = kTRUE; + } + + if ((!isGold) && (track.GetBackupTrack())) { + if ((track.GetBackupTrack()->GetNumberOfClusters() > foundMin) && ((track.GetBackupTrack()->GetChi2()/(track.GetBackupTrack()->GetNumberOfClusters()+1)) < 7)) { + seed->UpdateTrackParams(track.GetBackupTrack(),AliESDtrack::kTRDbackup); + isGold = kTRUE; + } + } + + //if ((track->StatusForTOF() > 0) && (track->GetNCross() == 0) && (Float_t(track->GetNumberOfClusters()) / Float_t(track->GetNExpected()) > 0.4)) { + //seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup); + //} + } + } + + // Propagation to the TOF (I.Belikov) + if (track.IsStopped() == kFALSE) { + Double_t xtof = 371.0; + Double_t xTOF0 = 370.0; + + Double_t c2 = track.GetSnp() + track.GetC() * (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()); + if (TMath::Abs(c2) >= 0.99) continue; + + //if (!PropagateToX(*track,xTOF0,fgkMaxStep)) { + // fHBackfit->Fill(7); + //delete track; + // continue; + //} + + Double_t ymax = xtof * TMath::Tan(0.5 * AliTRDgeometry::GetAlpha()); + Double_t y; + track.GetYAt(xtof,GetBz(),y); + if (y > ymax) { + if (!track.Rotate( AliTRDgeometry::GetAlpha())) continue; + }else if (y < -ymax) { + if (!track.Rotate(-AliTRDgeometry::GetAlpha())) continue; + } + + if (track.PropagateTo(xtof)) { + seed->UpdateTrackParams(&track, AliESDtrack::kTRDout); + track.UpdateESDtrack(seed); + } + } else { + if ((track.GetNumberOfClusters() > 15) && (track.GetNumberOfClusters() > 0.5*expectedClr)) { + seed->UpdateTrackParams(&track, AliESDtrack::kTRDout); + + track.UpdateESDtrack(seed); + } + } + + 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)); + + // run stand alone tracking + if (fReconstructor->IsSeeding()) Clusters2Tracks(event); + return 0; } @@ -437,55 +436,41 @@ Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event) Int_t nseed = 0; // contor for loaded seeds Int_t found = 0; // contor for updated TRD tracks - // Calibration monitor - AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance(); - if (!calibra) AliInfo("Could not get Calibra instance\n"); - AliTRDtrackV1 track; for (Int_t itrack = 0; itrack < event->GetNumberOfTracks(); itrack++) { AliESDtrack *seed = event->GetTrack(itrack); - new(&track) AliTRDtrackV1(*seed); + new(&track) AliTRDtrackV1(*seed); if (track.GetX() < 270.0) { seed->UpdateTrackParams(&track, AliESDtrack::kTRDbackup); - //AliInfo(Form("Remove for X = %7.3f [270.]\n", track.GetX())); - continue; + continue; } 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); - // do the propagation and processing + // do the propagation and processing Bool_t kUPDATE = kFALSE; - Double_t xTPC = 250.0; - if(FollowProlongation(track)){ - // computes PID for track - track.CookPID(); - // update calibration references using this track - if(calibra->GetHisto2d()) calibra->UpdateHistogramsV1(&track); - - // Prolongate to TPC - if (PropagateToX(track, xTPC, fgkMaxStep)) { // -with update - seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit); - 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++; - kUPDATE = kTRUE; - } - } - - // Prolongate to TPC without update - if(!kUPDATE) { + Double_t xTPC = 250.0; + if(FollowProlongation(track)){ + // Prolongate to TPC + if (PropagateToX(track, xTPC, fgkMaxStep)) { // -with update + seed->UpdateTrackParams(&track, AliESDtrack::kTRDrefit); + found++; + kUPDATE = kTRUE; + } + } + + // Prolongate to TPC without update + if(!kUPDATE) { AliTRDtrackV1 tt(*seed); if (PropagateToX(tt, xTPC, fgkMaxStep)) seed->UpdateTrackParams(&tt, AliESDtrack::kTRDrefit); } @@ -493,110 +478,98 @@ Int_t AliTRDtrackerV1::RefitInward(AliESDEvent *event) AliInfo(Form("Number of loaded seeds: %d",nseed)); AliInfo(Form("Number of found tracks from loaded seeds: %d",found)); - return 0; + return 0; } - //____________________________________________________________________ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) { -// Extrapolates the TRD track in the TPC direction. -// -// Parameters -// t : the TRD track which has to be extrapolated -// -// Output -// number of clusters attached to the track -// -// Detailed description -// -// Starting from current radial position of track this function -// extrapolates the track through the 6 TRD layers. The following steps -// are being performed for each plane: -// 1. prepare track: -// a. get plane limits in the local x direction -// b. check crossing sectors -// c. check track inclination -// 2. search tracklet in the tracker list (see GetTracklet() for details) -// 3. evaluate material budget using the geo manager -// 4. propagate and update track using the tracklet information. -// -// Debug level 2 -// - - //AliInfo(""); - Int_t nClustersExpected = 0; - Float_t clength = AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick(); - Int_t lastplane = 5; //GetLastPlane(&t); - for (Int_t iplane = lastplane; iplane >= 0; iplane--) { - //AliInfo(Form("plane %d", iplane)); - Int_t row1 = GetGlobalTimeBin(0, iplane, 0); // to be modified to the true time bin in the geometrical acceptance - //AliInfo(Form("row1 %d", row1)); - - // Propagate track close to the plane if neccessary - AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row1); - Double_t currentx = layer->GetX(); - if (currentx < (-fgkMaxStep + t.GetX())) - if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break; + // Extrapolates the TRD track in the TPC direction. + // + // Parameters + // t : the TRD track which has to be extrapolated + // + // Output + // number of clusters attached to the track + // + // Detailed description + // + // Starting from current radial position of track this function + // extrapolates the track through the 6 TRD layers. The following steps + // are being performed for each plane: + // 1. prepare track: + // a. get plane limits in the local x direction + // b. check crossing sectors + // c. check track inclination + // 2. search tracklet in the tracker list (see GetTracklet() for details) + // 3. evaluate material budget using the geo manager + // 4. propagate and update track using the tracklet information. + // + // Debug level 2 + // + + Int_t nClustersExpected = 0; + Int_t lastplane = 5; //GetLastPlane(&t); + for (Int_t iplane = lastplane; iplane >= 0; 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(); + // reject tracklets which are not considered for inward refit + if(x > t.GetX()+fgkMaxStep) continue; + // append tracklet to track + t.SetTracklet(tracklet, index); + + if (x < (t.GetX()-fgkMaxStep) && !PropagateToX(t, x+fgkMaxStep, fgkMaxStep)) break; if (!AdjustSector(&t)) break; - - Int_t row0 = GetGlobalTimeBin(0,iplane,GetTimeBinsPerPlane()-1); - //AliInfo(Form("row0 %d", row0)); - + // Start global position Double_t xyz0[3]; t.GetXYZ(xyz0); - // End global position - Double_t x = fTrSec[0]->GetLayer(row0)->GetX(), y, z; + // End global position + Double_t alpha = t.GetAlpha(), y, z; if (!t.GetProlongation(x,y,z)) break; Double_t xyz1[3]; - xyz1[0] = x * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); - xyz1[1] = +x * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha()); + xyz1[0] = x * TMath::Cos(alpha) - y * TMath::Sin(alpha); + 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); + 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 - //Int_t sector = t.GetSector(); - Int_t index = 0; - //AliInfo(Form("sector %d", sector)); - AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index); - //AliInfo(Form("tracklet %p @ %d", tracklet, index)); - if(!tracklet) continue; - //AliInfo(Form("reco %p", tracklet->GetRecoParam())); - t.SetTracklet(tracklet, iplane, index); - - t.PropagateTo(tracklet->GetX0() - clength, xx0, xrho); - if (!AdjustSector(&t)) break; - + // 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(); - } - } - -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() > 1){ - Int_t index; - for(int iplane=0; iplane<6; iplane++){ - AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index); - if(!tracklet) continue; - t.SetTracklet(tracklet, iplane, index); - } - - TTreeSRedirector &cstreamer = *fDebugStreamer; - cstreamer << "FollowProlongation" - << "ncl=" << nClustersExpected - << "track.=" << &t - << "\n"; - } -#endif + if (maxChi2 < 1e+10 && t.Update(tracklet, maxChi2)){ + nClustersExpected += tracklet->GetN(); + } + } + + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ + Int_t index; + for(int iplane=0; iplane<6; iplane++){ + AliTRDseedV1 *tracklet = GetTracklet(&t, iplane, index); + if(!tracklet) continue; + t.SetTracklet(tracklet, index); + } + + Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); + TTreeSRedirector &cstreamer = *fgDebugStreamer; + cstreamer << "FollowProlongation" + << "EventNumber=" << eventNumber + << "ncl=" << nClustersExpected + //<< "track.=" << &t + << "\n"; + } return nClustersExpected; @@ -605,196 +578,866 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t) //_____________________________________________________________________________ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t) { -// Extrapolates the TRD track in the TOF direction. -// -// Parameters -// t : the TRD track which has to be extrapolated -// -// Output -// number of clusters attached to the track -// -// Detailed description -// -// Starting from current radial position of track this function -// extrapolates the track through the 6 TRD layers. The following steps -// are being performed for each plane: -// 1. prepare track: -// a. get plane limits in the local x direction -// b. check crossing sectors -// c. check track inclination -// 2. build tracklet (see AliTRDseed::AttachClusters() for details) -// 3. evaluate material budget using the geo manager -// 4. propagate and update track using the tracklet information. -// -// Debug level 2 -// - - Int_t nClustersExpected = 0; - // Loop through the TRD planes - for (Int_t iplane = 0; iplane < AliTRDgeometry::Nplan(); iplane++) { - //AliInfo(Form("Processing plane %d ...", iplane)); - // Get the global time bin for the first local time bin of the given plane - Int_t row0 = GetGlobalTimeBin(0, iplane, fTimeBinsPerPlane-1); - - // Retrive first propagation layer in the chamber - AliTRDpropagationLayer *layer = fTrSec[0]->GetLayer(row0); - - // Get the X coordinates of the propagation layer for the first time bin - Double_t currentx = layer->GetX(); // what if X is not defined ??? - if (currentx < t.GetX()) continue; + // Extrapolates the TRD track in the TOF direction. + // + // Parameters + // t : the TRD track which has to be extrapolated + // + // Output + // number of clusters attached to the track + // + // Detailed description + // + // Starting from current radial position of track this function + // extrapolates the track through the 6 TRD layers. The following steps + // are being performed for each plane: + // 1. prepare track: + // a. get plane limits in the local x direction + // b. check crossing sectors + // c. check track inclination + // 2. build tracklet (see AliTRDseed::AttachClusters() for details) + // 3. evaluate material budget using the geo manager + // 4. propagate and update track using the tracklet information. + // + // Debug level 2 + // + + Int_t nClustersExpected = 0; + Double_t clength = 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); + } + + // 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 = tracklets[ilayer]; + if(!ptrTracklet){ + ptrTracklet = new(&tracklet) AliTRDseedV1(ilayer); + ptrTracklet->SetReconstructor(fReconstructor); + alpha = t.GetAlpha(); + Int_t sector = Int_t(alpha/AliTRDgeometry::GetAlpha() + (alpha>0. ? 0 : AliTRDgeometry::kNsector)); + + if(!fTrSec[sector].GetNChambers()) continue; + + if((x = fTrSec[sector].GetX(ilayer)) < 1.) continue; + + if (!t.GetProlongation(x, y, z)) return -1/*nClustersExpected*/; + Int_t stack = fGeom->GetStack(z, ilayer); + Int_t nCandidates = stack >= 0 ? 1 : 2; + z -= stack >= 0 ? 0. : 4.; + + for(int icham=0; ichamGetStack(z, ilayer)) < 0) continue; + + if(!(chamber = fTrSec[sector].GetChamber(stack, ilayer))) continue; + + if(chamber->GetNClusters() < fgNTimeBins*fReconstructor->GetRecoParam() ->GetFindableClusters()) continue; + + x = chamber->GetX(); + + AliTRDpadPlane *pp = fGeom->GetPadPlane(ilayer, stack); + tracklet.SetTilt(TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle())); + tracklet.SetPadLength(pp->GetLengthIPad()); + tracklet.SetDetector(chamber->GetDetector()); + tracklet.SetX0(x); + if(!tracklet.Init(&t)){ + t.SetStopped(kTRUE); + return nClustersExpected; + } + if(!tracklet.AttachClustersIter(chamber, 1000./*, kTRUE*/)) continue; + tracklet.Init(&t); + + 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)) return -1/*nClustersExpected*/; + if(!AdjustSector(&t)) return -1/*nClustersExpected*/; + if(TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/; + continue; + } - // Get the global time bin for the last local time bin of the given plane - Int_t row1 = GetGlobalTimeBin(0, iplane, 0); + // Propagate closer to the current chamber if neccessary + x -= clength; + if (x > (fgkMaxStep + t.GetX()) && !PropagateToX(t, x-fgkMaxStep, fgkMaxStep)) return -1/*nClustersExpected*/; + if (!AdjustSector(&t)) return -1/*nClustersExpected*/; + if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) return -1/*nClustersExpected*/; - // Propagate closer to the current chamber if neccessary - if (currentx > (fgkMaxStep + t.GetX()) && !PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break; - - // Rotate track to adjacent sector if neccessary - if (!AdjustSector(&t)) break; - Int_t sector = Int_t(TMath::Abs(t.GetAlpha()/AliTRDgeometry::GetAlpha())); - if(t.GetAlpha() < 0) sector = AliTRDgeometry::Nsect() - sector-1; - - //AliInfo(Form("sector %d [%f]", sector, t.GetAlpha())); - - // Check whether azimuthal angle is getting too large - if (TMath::Abs(t.GetSnp()) > fgkMaxSnp) break; - + // load tracklet to the tracker and the track + ptrTracklet = SetTracklet(ptrTracklet); + t.SetTracklet(ptrTracklet, fTracklets->GetEntriesFast()-1); + + + // Calculate the mean material budget along the path inside the chamber //Calculate global entry and exit positions of the track in chamber (only track prolongation) Double_t xyz0[3]; // entry point - t.GetXYZ(xyz0); - //printf("Entry global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz0[0], xyz0[1], xyz0[2]); - - // Get local Y and Z at the X-position of the end of the chamber - Double_t x0 = fTrSec[sector]->GetLayer(row1)->GetX(), y, z; - if (!t.GetProlongation(x0, y, z)) break; - //printf("Exit local x[%7.3f] y[%7.3f] z[%7.3f]\n", x0, y, z); - - Double_t xyz1[3]; // exit point - xyz1[0] = x0 * TMath::Cos(t.GetAlpha()) - y * TMath::Sin(t.GetAlpha()); - xyz1[1] = +x0 * TMath::Sin(t.GetAlpha()) + y * TMath::Cos(t.GetAlpha()); + t.GetXYZ(xyz0); + alpha = t.GetAlpha(); + x = ptrTracklet->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; - - //printf("Exit global x[%7.3f] y[%7.3f] z[%7.3f]\n", xyz1[0], xyz1[1], xyz1[2]); - // Find tracklet along the path inside the chamber - AliTRDseedV1 tracklet(*t.GetTracklet(iplane)); - // if the track is not already build (e.g. stand alone tracker) we build it now. - if(!tracklet.GetN()){ // a better check has to be implemented TODO!!!!!!! - - //AliInfo(Form("Building tracklet for plane %d ...", iplane)); - // check if we are inside detection volume - Int_t ichmb = fGeom->GetChamber(xyz0[2], iplane); - if(ichmb<0) ichmb = fGeom->GetChamber(xyz1[2], iplane); - if(ichmb<0){ - // here we should decide what to do with the track. The space between the pads in 2 chambers is 4cm+. Is it making sense to continue building the tracklet here TODO???? - AliWarning(Form("Track prolongated in the interspace between TRD detectors in plane %d. Skip plane. To be fixed !", iplane)); - continue; - } - - // temporary until the functionalities of AliTRDpropagationLayer and AliTRDstackLayer are merged TODO - AliTRDpadPlane *pp = fGeom->GetPadPlane(iplane, ichmb); - Int_t nrows = pp->GetNrows(); - Double_t stacklength = pp->GetRow0ROC() - pp->GetRowEndROC();/*(nrows - 2) * pp->GetLengthIPad() + 2 * pp->GetLengthOPad() + (nrows - 1) * pp->GetRowSpacing();*/ - Double_t z0 = fGeom->GetRow0(iplane, ichmb, 0); - - Int_t nClustersChmb = 0; - AliTRDstackLayer stackLayer[35]; - for(int itb=0; itbGetLayer(row1 - itb))); - stackLayer[itb] = ksmLayer; -#ifdef DEBUG - stackLayer[itb].SetDebugStream(fDebugStreamer); -#endif - stackLayer[itb].SetRange(z0 - stacklength, stacklength); - stackLayer[itb].SetSector(sector); - stackLayer[itb].SetStackNr(ichmb); - stackLayer[itb].SetNRows(nrows); - stackLayer[itb].SetRecoParam(fRecoParam); - stackLayer[itb].BuildIndices(); - nClustersChmb += stackLayer[itb].GetNClusters(); - } - if(!nClustersChmb) continue; - //AliInfo(Form("Detector p[%d] c[%d]. Building tracklet from %d clusters ... ", iplane, ichmb, nClustersChmb)); - - tracklet.SetRecoParam(fRecoParam); - tracklet.SetTilt(TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle())); - tracklet.SetPadLength(pp->GetLengthIPad()); - tracklet.SetPlane(iplane); - //Int_t tbRange = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det)); - //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]); - //tracklet.SetNTimeBinsRange(tbRange); - tracklet.SetX0(x0); - tracklet.Init(&t); - if(!tracklet.AttachClustersIter(stackLayer, 1000.)) continue; - tracklet.Init(&t); - - //if(!tracklet.AttachClusters(stackLayer, kTRUE)) continue; - //if(!tracklet.Fit()) continue; - } - Int_t ncl = tracklet.GetN(); - //AliInfo(Form("N clusters %d", ncl)); - - // Discard tracklet if bad quality. - //Check if this condition is not already checked during building of the tracklet - if(ncl < fTimeBinsPerPlane * fRecoParam->GetFindableClusters()){ - //AliInfo(Form("Discard tracklet for %d nclusters", ncl)); - continue; - } - - // load tracklet to the tracker and the track - Int_t index = SetTracklet(&tracklet); - t.SetTracklet(&tracklet, iplane, index); - - // Calculate the mean material budget along the path inside the chamber 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(tracklet.GetX0(), xx0, xrho); - if (!AdjustSector(&t)) break; - Double_t maxChi2 = t.GetPredictedChi2(&tracklet); - if (maxChi2<1e+10 && t.Update(&tracklet, maxChi2)){ - nClustersExpected += ncl; - } - // Reset material budget if 2 consecutive gold - if(iplane>0 && ncl + t.GetTracklet(iplane-1)->GetN() > 20) t.SetBudget(2, 0.); - - // Make backup of the track until is gold - // TO DO update quality check of the track. - // consider comparison with fTimeBinsRange - Float_t ratio0 = ncl / Float_t(fTimeBinsPerPlane); - //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1); + + // Propagate and update track + if (!t.PropagateTo(x, xx0, xrho)) return -1/*nClustersExpected*/; + if (!AdjustSector(&t)) return -1/*nClustersExpected*/; + Double_t maxChi2 = t.GetPredictedChi2(ptrTracklet); + if (!t.Update(ptrTracklet, maxChi2)) return -1/*nClustersExpected*/; + if (maxChi2<1e+10) { + nClustersExpected += ptrTracklet->GetN(); + //t.SetTracklet(&tracklet, index); + } + // Reset material budget if 2 consecutive gold + if(ilayer>0 && t.GetTracklet(ilayer-1) && ptrTracklet->GetN() + t.GetTracklet(ilayer-1)->GetN() > 20) t.SetBudget(2, 0.); + + // Make backup of the track until is gold + // TO DO update quality check of the track. + // consider comparison with fTimeBinsRange + Float_t ratio0 = ptrTracklet->GetN() / Float_t(fgNTimeBins); + //Float_t ratio1 = Float_t(t.GetNumberOfClusters()+1) / Float_t(t.GetNExpected()+1); //printf("tracklet.GetChi2() %f [< 18.0]\n", tracklet.GetChi2()); - //printf("ratio0 %f [> 0.8]\n", ratio0); - //printf("ratio1 %f [> 0.6]\n", ratio1); - //printf("ratio0+ratio1 %f [> 1.5]\n", ratio0+ratio1); - //printf("t.GetNCross() %d [== 0]\n", t.GetNCross()); - //printf("TMath::Abs(t.GetSnp()) %f [< 0.85]\n", TMath::Abs(t.GetSnp())); - //printf("t.GetNumberOfClusters() %d [> 20]\n", t.GetNumberOfClusters()); + //printf("ratio0 %f [> 0.8]\n", ratio0); + //printf("ratio1 %f [> 0.6]\n", ratio1); + //printf("ratio0+ratio1 %f [> 1.5]\n", ratio0+ratio1); + //printf("t.GetNCross() %d [== 0]\n", t.GetNCross()); + //printf("TMath::Abs(t.GetSnp()) %f [< 0.85]\n", TMath::Abs(t.GetSnp())); + //printf("t.GetNumberOfClusters() %d [> 20]\n", t.GetNumberOfClusters()); - if (//(tracklet.GetChi2() < 18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update + if (//(tracklet.GetChi2() < 18.0) && TO DO check with FindClusters and move it to AliTRDseed::Update (ratio0 > 0.8) && //(ratio1 > 0.6) && //(ratio0+ratio1 > 1.5) && (t.GetNCross() == 0) && (TMath::Abs(t.GetSnp()) < 0.85) && (t.GetNumberOfClusters() > 20)) t.MakeBackupTrack(); - - } // end planes loop + + } // end layers loop + + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ + TTreeSRedirector &cstreamer = *fgDebugStreamer; + Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); + //AliTRDtrackV1 *debugTrack = new AliTRDtrackV1(t); + //debugTrack->SetOwner(); + cstreamer << "FollowBackProlongation" + << "EventNumber=" << eventNumber + << "ncl=" << nClustersExpected + //<< "track.=" << debugTrack + << "\n"; + } + + return nClustersExpected; +} + +//_________________________________________________________________________ +Float_t AliTRDtrackerV1::FitRieman(AliTRDseedV1 *tracklets, Double_t *chi2, Int_t *planes){ + // + // Fits a Riemann-circle to the given points without tilting pad correction. + // The fit is performed using an instance of the class AliRieman (equations + // and transformations see documentation of this class) + // Afterwards all the tracklets are Updated + // + // Parameters: - Array of tracklets (AliTRDseedV1) + // - Storage for the chi2 values (beginning with direction z) + // - Seeding configuration + // Output: - The curvature + // + AliRieman *fitter = AliTRDtrackerV1::GetRiemanFitter(); + fitter->Reset(); + Int_t allplanes[] = {0, 1, 2, 3, 4, 5}; + Int_t *ppl = &allplanes[0]; + Int_t maxLayers = 6; + if(planes){ + maxLayers = 4; + ppl = planes; + } + 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->Update(); + // Set the reference position of the fit and calculate the chi2 values + memset(chi2, 0, sizeof(Double_t) * 2); + for(Int_t il = 0; il < maxLayers; il++){ + // Reference positions + tracklets[ppl[il]].Init(fitter); + + // chi2 + if((!tracklets[ppl[il]].IsOK()) && (!planes)) continue; + chi2[0] += tracklets[ppl[il]].GetChi2Y(); + chi2[1] += tracklets[ppl[il]].GetChi2Z(); + } + return fitter->GetC(); +} + +//_________________________________________________________________________ +void AliTRDtrackerV1::FitRieman(AliTRDcluster **seedcl, Double_t chi2[2]) +{ + // + // Performs a Riemann helix fit using the seedclusters as spacepoints + // Afterwards the chi2 values are calculated and the seeds are updated + // + // Parameters: - The four seedclusters + // - The tracklet array (AliTRDseedV1) + // - The seeding configuration + // - Chi2 array + // + // debug level 2 + // + AliRieman *fitter = AliTRDtrackerV1::GetRiemanFitter(); + fitter->Reset(); + for(Int_t i = 0; i < 4; i++) + fitter->AddPoint(seedcl[i]->GetX(), seedcl[i]->GetY(), seedcl[i]->GetZ(), 1, 10); + fitter->Update(); + + + // Update the seed and calculated the chi2 value + chi2[0] = 0; chi2[1] = 0; + for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){ + // chi2 + chi2[0] += (seedcl[ipl]->GetZ() - fitter->GetZat(seedcl[ipl]->GetX())) * (seedcl[ipl]->GetZ() - fitter->GetZat(seedcl[ipl]->GetX())); + chi2[1] += (seedcl[ipl]->GetY() - fitter->GetYat(seedcl[ipl]->GetX())) * (seedcl[ipl]->GetY() - fitter->GetYat(seedcl[ipl]->GetX())); + } +} + + +//_________________________________________________________________________ +Float_t AliTRDtrackerV1::FitTiltedRiemanConstraint(AliTRDseedV1 *tracklets, Double_t zVertex) +{ + // + // Fits a helix to the clusters. Pad tilting is considered. As constraint it is + // assumed that the vertex position is set to 0. + // This method is very usefull for high-pt particles + // Basis for the fit: (x - x0)^2 + (y - y0)^2 - R^2 = 0 + // x0, y0: Center of the circle + // Measured y-position: ymeas = y - tan(phiT)(zc - zt) + // zc: center of the pad row + // Equation which has to be fitted (after transformation): + // a + b * u + e * v + 2*(ymeas + tan(phiT)(z - zVertex))*t = 0 + // Transformation: + // t = 1/(x^2 + y^2) + // u = 2 * x * t + // v = 2 * x * tan(phiT) * t + // Parameters in the equation: + // a = -1/y0, b = x0/y0, e = dz/dx + // + // The Curvature is calculated by the following equation: + // - curv = a/Sqrt(b^2 + 1) = 1/R + // Parameters: - the 6 tracklets + // - the Vertex constraint + // Output: - the Chi2 value of the track + // + // debug level 5 + // + + TLinearFitter *fitter = GetTiltedRiemanFitterConstraint(); + fitter->StoreData(kTRUE); + fitter->ClearPoints(); + AliTRDcluster *cl = 0x0; + + Float_t x, y, z, w, t, error, tilt; + Double_t uvt[2]; + 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++){ + if(!tracklets[ilr].IsUsable(itb)) continue; + cl = tracklets[ilr].GetClusters(itb); + x = cl->GetX(); + y = cl->GetY(); + z = cl->GetZ(); + tilt = tracklets[ilr].GetTilt(); + // Transformation + t = 1./(x * x + y * y); + uvt[0] = 2. * x * t; + uvt[1] = 2. * x * t * tilt ; + w = 2. * (y + tilt * (z - zVertex)) * t; + error = 2. * 0.2 * t; + fitter->AddPoint(uvt, w, error); + nPoints++; + } + } + fitter->Eval(); + + // Calculate curvature + Double_t a = fitter->GetParameter(0); + Double_t b = fitter->GetParameter(1); + Double_t curvature = a/TMath::Sqrt(b*b + 1); + + Float_t chi2track = fitter->GetChisquare()/Double_t(nPoints); + for(Int_t ip = 0; ip < AliTRDtrackerV1::kNPlanes; ip++) + tracklets[ip].SetCC(curvature); + +/* 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); + Double_t zref = slope * xref; + Float_t chi2Z = CalculateChi2Z(tracklets, zref, slope, xref); + Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); + Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); + TTreeSRedirector &treeStreamer = *fgDebugStreamer; + treeStreamer << "FitTiltedRiemanConstraint" + << "EventNumber=" << eventNumber + << "CandidateNumber=" << candidateNumber + << "Curvature=" << curvature + << "Chi2Track=" << chi2track + << "Chi2Z=" << chi2Z + << "zref=" << zref + << "\n"; + }*/ + return chi2track; +} + +//_________________________________________________________________________ +Float_t AliTRDtrackerV1::FitTiltedRieman(AliTRDseedV1 *tracklets, Bool_t sigError) +{ + // + // Performs a Riemann fit taking tilting pad correction into account + // The equation of a Riemann circle, where the y position is substituted by the + // measured y-position taking pad tilting into account, has to be transformed + // into a 4-dimensional hyperplane equation + // Riemann circle: (x-x0)^2 + (y-y0)^2 -R^2 = 0 + // Measured y-Position: ymeas = y - tan(phiT)(zc - zt) + // zc: center of the pad row + // zt: z-position of the track + // The z-position of the track is assumed to be linear dependent on the x-position + // Transformed equation: a + b * u + c * t + d * v + e * w - 2 * (ymeas + tan(phiT) * zc) * t = 0 + // Transformation: u = 2 * x * t + // v = 2 * tan(phiT) * t + // w = 2 * tan(phiT) * (x - xref) * t + // t = 1 / (x^2 + ymeas^2) + // Parameters: a = -1/y0 + // b = x0/y0 + // c = (R^2 -x0^2 - y0^2)/y0 + // d = offset + // e = dz/dx + // If the offset respectively the slope in z-position is impossible, the parameters are fixed using + // results from the simple riemann fit. Afterwards the fit is redone. + // The curvature is calculated according to the formula: + // curv = a/(1 + b^2 + c*a) = 1/R + // + // Paramters: - Array of tracklets (connected to the track candidate) + // - Flag selecting the error definition + // Output: - Chi2 values of the track (in Parameter list) + // + TLinearFitter *fitter = GetTiltedRiemanFitter(); + fitter->StoreData(kTRUE); + fitter->ClearPoints(); + AliTRDLeastSquare zfitter; + AliTRDcluster *cl = 0x0; + + Double_t xref = CalculateReferenceX(tracklets); + Double_t x, y, z, t, tilt, dx, w, we; + Double_t uvt[4]; + Int_t nPoints = 0; + // Containers for Least-square fitter + for(Int_t ipl = 0; ipl < kNPlanes; ipl++){ + if(!tracklets[ipl].IsOK()) continue; + for(Int_t itb = 0; itb < 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; + // Transformation + t = 1./(x*x + y*y); + uvt[0] = 2. * x * t; + uvt[1] = t; + uvt[2] = 2. * tilt * t; + uvt[3] = 2. * tilt * dx * t; + w = 2. * (y + tilt*z) * t; + // error definition changes for the different calls + we = 2. * t; + we *= sigError ? tracklets[ipl].GetSigmaY() : 0.2; + fitter->AddPoint(uvt, w, we); + zfitter.AddPoint(&x, z, static_cast(TMath::Sqrt(cl->GetSigmaZ2()))); + nPoints++; + } + } + fitter->Eval(); + zfitter.Eval(); + + Double_t offset = fitter->GetParameter(3); + Double_t slope = fitter->GetParameter(4); + + // Linear fitter - not possible to make boundaries + // Do not accept non possible z and dzdx combinations + Bool_t acceptablez = kTRUE; + Double_t zref = 0.0; + 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) + acceptablez = kFALSE; + } + if (!acceptablez) { + Double_t dzmf = zfitter.GetFunctionParameter(1); + Double_t zmf = zfitter.GetFunctionValue(&xref); + fgTiltedRieman->FixParameter(3, zmf); + fgTiltedRieman->FixParameter(4, dzmf); + fitter->Eval(); + fitter->ReleaseParameter(3); + fitter->ReleaseParameter(4); + offset = fitter->GetParameter(3); + slope = fitter->GetParameter(4); + } + + // Calculate Curvarture + Double_t a = fitter->GetParameter(0); + Double_t b = fitter->GetParameter(1); + Double_t c = fitter->GetParameter(2); + Double_t curvature = 1.0 + b*b - c*a; + if (curvature > 0.0) + curvature = a / TMath::Sqrt(curvature); + + Double_t chi2track = fitter->GetChisquare()/Double_t(nPoints); + + // Update the tracklets + Double_t dy, dz; + for(Int_t iLayer = 0; iLayer < AliTRDtrackerV1::kNPlanes; iLayer++) { + + x = tracklets[iLayer].GetX0(); + y = 0; + z = 0; + dy = 0; + dz = 0; + + // y: R^2 = (x - x0)^2 + (y - y0)^2 + // => y = y0 +/- Sqrt(R^2 - (x - x0)^2) + // R = Sqrt() = 1/Curvature + // => y = y0 +/- Sqrt(1/Curvature^2 - (x - x0)^2) + Double_t res = (x * a + b); // = (x - x0)/y0 + res *= res; + res = 1.0 - c * a + b * b - res; // = (R^2 - (x - x0)^2)/y0^2 + if (res >= 0) { + res = TMath::Sqrt(res); + y = (1.0 - res) / a; + } + + // dy: R^2 = (x - x0)^2 + (y - y0)^2 + // => y = +/- Sqrt(R^2 - (x - x0)^2) + y0 + // => dy/dx = (x - x0)/Sqrt(R^2 - (x - x0)^2) + // Curvature: cr = 1/R = a/Sqrt(1 + b^2 - c*a) + // => dy/dx = (x - x0)/(1/(cr^2) - (x - x0)^2) + Double_t x0 = -b / a; + if (-c * a + b * b + 1 > 0) { + if (1.0/(curvature * curvature) - (x - x0) * (x - x0) > 0.0) { + Double_t yderiv = (x - x0) / TMath::Sqrt(1.0/(curvature * curvature) - (x - x0) * (x - x0)); + if (a < 0) yderiv *= -1.0; + dy = yderiv; + } + } + z = offset + slope * (x - xref); + dz = slope; + tracklets[iLayer].SetYref(0, y); + tracklets[iLayer].SetYref(1, dy); + tracklets[iLayer].SetZref(0, z); + tracklets[iLayer].SetZref(1, dz); + tracklets[iLayer].SetC(curvature); + tracklets[iLayer].SetChi2(chi2track); + } + +/* if(fReconstructor->GetStreamLevel() >=5){ + TTreeSRedirector &cstreamer = *fgDebugStreamer; + Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); + Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); + Double_t chi2z = CalculateChi2Z(tracklets, offset, slope, xref); + cstreamer << "FitTiltedRieman0" + << "EventNumber=" << eventNumber + << "CandidateNumber=" << candidateNumber + << "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(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; ipStoreData(kTRUE); + fitter->ClearPoints(); + AliTRDLeastSquare zfitter; + AliTRDcluster *cl = 0x0; + + AliTRDseedV1 work[kNPlanes], *tracklet = 0x0; + if(!tracklets){ + for(Int_t ipl = 0; ipl < kNPlanes; ipl++){ + if(!(tracklet = track->GetTracklet(ipl))) continue; + if(!tracklet->IsOK()) continue; + new(&work[ipl]) AliTRDseedV1(*tracklet); + } + tracklets = &work[0]; + } + + Double_t xref = CalculateReferenceX(tracklets); + Double_t x, y, z, t, tilt, dx, w, we; + Double_t uvt[4]; + Int_t nPoints = 0; + // Containers for Least-square fitter + for(Int_t ipl = 0; ipl < kNPlanes; ipl++){ + if(!tracklets[ipl].IsOK()) continue; + for(Int_t itb = 0; itb < 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; + // Transformation + t = 1./(x*x + y*y); + uvt[0] = 2. * x * t; + uvt[1] = t; + uvt[2] = 2. * tilt * t; + uvt[3] = 2. * tilt * dx * t; + w = 2. * (y + tilt*z) * t; + // error definition changes for the different calls + we = 2. * t; + we *= sigError ? tracklets[ipl].GetSigmaY() : 0.2; + fitter->AddPoint(uvt, w, we); + zfitter.AddPoint(&x, z, static_cast(TMath::Sqrt(cl->GetSigmaZ2()))); + nPoints++; + } + } + if(fitter->Eval()) return 1.E10; + + Double_t z0 = fitter->GetParameter(3); + Double_t dzdx = fitter->GetParameter(4); + + + // Linear fitter - not possible to make boundaries + // Do not accept non possible z and dzdx combinations + Bool_t accept = kTRUE; + Double_t zref = 0.0; + for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) { + if(!tracklets[iLayer].IsOK()) continue; + zref = z0 + dzdx * (tracklets[iLayer].GetX0() - xref); + if (TMath::Abs(tracklets[iLayer].GetZProb() - zref) > tracklets[iLayer].GetPadLength() * 0.5 + 1.0) + accept = kFALSE; + } + if (!accept) { + zfitter.Eval(); + Double_t dzmf = zfitter.GetFunctionParameter(1); + Double_t zmf = zfitter.GetFunctionValue(&xref); + fitter->FixParameter(3, zmf); + fitter->FixParameter(4, dzmf); + fitter->Eval(); + fitter->ReleaseParameter(3); + fitter->ReleaseParameter(4); + z0 = fitter->GetParameter(3); // = zmf ? + dzdx = fitter->GetParameter(4); // = dzmf ? + } + + // Calculate Curvature + Double_t a = fitter->GetParameter(0); + Double_t b = fitter->GetParameter(1); + Double_t c = fitter->GetParameter(2); + Double_t y0 = 1. / a; + Double_t x0 = -b * y0; + Double_t tmp = y0*y0 + x0*x0 - c*y0; + if(tmp<=0.) return 1.E10; + Double_t R = TMath::Sqrt(tmp); + Double_t C = 1.0 + b*b - c*a; + if (C > 0.0) C = a / TMath::Sqrt(C); + + // Calculate chi2 of the fit + Double_t chi2 = fitter->GetChisquare()/Double_t(nPoints); + + // Update the tracklets + if(!track){ + for(Int_t ip = 0; ip < kNPlanes; ip++) { + x = tracklets[ip].GetX0(); + tmp = R*R-(x-x0)*(x-x0); + if(tmp <= 0.) continue; + tmp = TMath::Sqrt(tmp); + + // y: R^2 = (x - x0)^2 + (y - y0)^2 + // => y = y0 +/- Sqrt(R^2 - (x - x0)^2) + tracklets[ip].SetYref(0, y0 - (y0>0.?1.:-1)*tmp); + // => dy/dx = (x - x0)/Sqrt(R^2 - (x - x0)^2) + tracklets[ip].SetYref(1, (x - x0) / tmp); + tracklets[ip].SetZref(0, z0 + dzdx * (x - xref)); + tracklets[ip].SetZref(1, dzdx); + tracklets[ip].SetC(C); + tracklets[ip].SetChi2(chi2); + } + } + + //update track points array + if(np && points){ + Float_t xyz[3]; + for(int ip=0; ip0.?1.:-1.)*TMath::Sqrt(R*R-(xyz[0]-x0)*(xyz[0]-x0)); + xyz[2] = z0 + dzdx * (xyz[0] - xref); + points[ip].SetXYZ(xyz); + } + } + + return chi2; +} + + +//____________________________________________________________________ +Double_t AliTRDtrackerV1::FitKalman(AliTRDtrackV1 *track, AliTRDseedV1 *tracklets, Bool_t up, Int_t np, AliTrackPoint *points) +{ +// Kalman filter implementation for the TRD. +// It returns the positions of the fit in the array "points" +// +// Author : A.Bercuci@gsi.de -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() > 1){ - TTreeSRedirector &cstreamer = *fDebugStreamer; - cstreamer << "FollowBackProlongation" - << "ncl=" << nClustersExpected - << "track.=" << &t - << "\n"; - } -#endif + //printf("Start track @ x[%f]\n", track->GetX()); - return nClustersExpected; + //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; + + //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) +{ + // + // Calculates the chi2-value of the track in z-Direction including tilting pad correction. + // A linear dependence on the x-value serves as a model. + // The parameters are related to the tilted Riemann fit. + // Parameters: - Array of tracklets (AliTRDseedV1) related to the track candidate + // - the offset for the reference x + // - the slope + // - the reference x position + // Output: - The Chi2 value of the track in z-Direction + // + Float_t chi2Z = 0, nLayers = 0; + for (Int_t iLayer = 0; iLayer < AliTRDgeometry::kNlayer; iLayer++) { + if(!tracklets[iLayer].IsOK()) continue; + Double_t z = offset + slope * (tracklets[iLayer].GetX0() - xref); + chi2Z += TMath::Abs(tracklets[iLayer].GetMeanz() - z); + nLayers++; + } + chi2Z /= TMath::Max((nLayers - 3.0),1.0); + return chi2Z; } //_____________________________________________________________________________ @@ -844,7 +1487,7 @@ 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])) { @@ -863,108 +1506,263 @@ Int_t AliTRDtrackerV1::PropagateToX(AliTRDtrackV1 &t, Double_t xToGo, Double_t m } -//____________________________________________________________________ -void AliTRDtrackerV1::UnloadClusters() -{ + +//_____________________________________________________________________________ +Int_t AliTRDtrackerV1::ReadClusters(TClonesArray* &array, TTree *clusterTree) const +{ // - // Clears the arrays of clusters and tracks. Resets sectors and timebins + // Reads AliTRDclusters from the file. + // The names of the cluster tree and branches + // should match the ones used in AliTRDclusterizer::WriteClusters() // - Int_t i; - Int_t nentr; + Int_t nsize = Int_t(clusterTree->GetTotBytes() / (sizeof(AliTRDcluster))); + TObjArray *clusterArray = new TObjArray(nsize+1000); + + TBranch *branch = clusterTree->GetBranch("TRDcluster"); + if (!branch) { + AliError("Can't get the branch !"); + return 1; + } + branch->SetAddress(&clusterArray); + + if(!fClusters){ + Float_t nclusters = fReconstructor->GetRecoParam()->GetNClusters(); + if(fReconstructor->IsHLT()) nclusters /= AliTRDgeometry::kNsector; + array = new TClonesArray("AliTRDcluster", Int_t(nclusters)); + array->SetOwner(kTRUE); + } + + // Loop through all entries in the tree + Int_t nEntries = (Int_t) clusterTree->GetEntries(); + Int_t nbytes = 0; + Int_t ncl = 0; + AliTRDcluster *c = 0x0; + for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { + // Import the tree + nbytes += clusterTree->GetEvent(iEntry); + + // Get the number of points in the detector + Int_t nCluster = clusterArray->GetEntriesFast(); + for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { + if(!(c = (AliTRDcluster *) clusterArray->UncheckedAt(iCluster))) continue; + c->SetInChamber(); + new((*fClusters)[ncl++]) AliTRDcluster(*c); + delete (clusterArray->RemoveAt(iCluster)); + } + + } + delete clusterArray; + + return 0; +} + +//_____________________________________________________________________________ +Int_t AliTRDtrackerV1::LoadClusters(TTree *cTree) +{ + // + // Fills clusters into TRD tracking sectors + // + + if(!fReconstructor->IsWritingClusters()){ + fClusters = AliTRDReconstructor::GetClusters(); + } else { + if (ReadClusters(fClusters, cTree)) { + AliError("Problem with reading the clusters !"); + return 1; + } + } + SetClustersOwner(); - nentr = fClusters->GetEntriesFast(); - for (i = 0; i < nentr; i++) { - delete fClusters->RemoveAt(i); + if(!fClusters || !fClusters->GetEntriesFast()){ + AliInfo("No TRD clusters"); + return 1; } - fNclusters = 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(fTracklets){ - for (i = 0; i < fTracklets->GetEntriesFast(); i++) delete fTracklets->RemoveAt(i); - } + if(!clusters || !clusters->GetEntriesFast()){ + AliInfo("No TRD clusters"); + return 1; + } + + 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)); - nentr = fSeeds->GetEntriesFast(); - for (i = 0; i < nentr; i++) { - delete fSeeds->RemoveAt(i); + 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->GetStack(detector); + Int_t layer = fGeom->GetLayer(detector); + + fTrSec[sector].GetChamber(stack, layer, kTRUE)->InsertCluster(c, icl); } - nentr = fTracks->GetEntriesFast(); - for (i = 0; i < nentr; i++) { - delete fTracks->RemoveAt(i); + const AliTRDCalDet *cal = AliTRDcalibDB::Instance()->GetT0Det(); + for(int isector =0; isectorGetNumberOfLayers(); pl++) { - fTrSec[i]->GetLayer(pl)->Clear(); - } + return nin; +} + + + +//____________________________________________________________________ +void AliTRDtrackerV1::UnloadClusters() +{ + // + // Clears the arrays of clusters and tracks. Resets sectors and timebins + // + + if(fTracks) fTracks->Delete(); + if(fTracklets) fTracklets->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(); + + // Increment the Event Number + AliTRDtrackerDebug::SetEventNumber(AliTRDtrackerDebug::GetEventNumber() + 1); +} + +//_____________________________________________________________________________ +Bool_t AliTRDtrackerV1::AdjustSector(AliTRDtrackV1 *track) +{ + // + // Rotates the track when necessary + // + + 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; + } + } + else if (y < -ymax) { + if (!track->Rotate(-alpha)) { + return kFALSE; + } + } + + return kTRUE; + } + //____________________________________________________________________ AliTRDseedV1* AliTRDtrackerV1::GetTracklet(AliTRDtrackV1 *track, Int_t p, Int_t &idx) { -// Find tracklet for TRD track -// Parameters -// - track -// - sector -// - plane -// - index -// Output -// tracklet -// index -// Detailed description -// - idx = track->GetTrackletIndex(p); - //AliInfo(Form("looking for tracklet in plane %d idx %d [%d]", p, idx, track->GetTrackletIndex(p))); - AliTRDseedV1 *tracklet = idx<0 ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx); - //AliInfo(Form("found 0x%x @ %d", tracklet, idx)); - -// Int_t *index = track->GetTrackletIndexes(); -// for (UInt_t i = 0; i < 6; i++) AliInfo(Form("index[%d] = %d", i, index[i])); -// -// for (UInt_t i = 0; i < 6/*kMaxTimeBinIndex*/; i++) { -// if (index[i] < 0) continue; -// -// tracklet = (AliTRDseedV1*)fTracklets->UncheckedAt(index[i]); -// if(!tracklet) break; -// -// if(tracklet->GetPlane() != p) continue; -// -// idx = index[i]; -// } + // Find tracklet for TRD track + // Parameters + // - track + // - sector + // - plane + // - index + // Output + // tracklet + // index + // Detailed description + // + idx = track->GetTrackletIndex(p); + AliTRDseedV1 *tracklet = (idx==0xffff) ? 0x0 : (AliTRDseedV1*)fTracklets->UncheckedAt(idx); - return tracklet; + return tracklet; } //____________________________________________________________________ -Int_t AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet) +AliTRDseedV1* AliTRDtrackerV1::SetTracklet(AliTRDseedV1 *tracklet) { -// Add this tracklet to the list of tracklets stored in the tracker -// -// Parameters -// - tracklet : pointer to the tracklet to be added to the list -// -// Output -// - the index of the new tracklet in the tracker tracklets list -// -// Detailed description -// Build the tracklets list if it is not yet created (late initialization) -// and adds the new tracklet to the list. -// - if(!fTracklets){ - fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsect()*kMaxTracksStack); - fTracklets->SetOwner(kTRUE); - } - Int_t nentries = fTracklets->GetEntriesFast(); - new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet); - return nentries; + // Add this tracklet to the list of tracklets stored in the tracker + // + // Parameters + // - tracklet : pointer to the tracklet to be added to the list + // + // Output + // - the index of the new tracklet in the tracker tracklets list + // + // Detailed description + // Build the tracklets list if it is not yet created (late initialization) + // and adds the new tracklet to the list. + // + if(!fTracklets){ + fTracklets = new TClonesArray("AliTRDseedV1", AliTRDgeometry::Nsector()*kMaxTracksStack); + fTracklets->SetOwner(kTRUE); + } + Int_t nentries = fTracklets->GetEntriesFast(); + return new ((*fTracklets)[nentries]) AliTRDseedV1(*tracklet); } //____________________________________________________________________ -Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *sector - , AliESDEvent *esd) +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) { // // Steer tracking for one SM. @@ -982,60 +1780,42 @@ Int_t AliTRDtrackerV1::Clusters2TracksSM(AliTRDtracker::AliTRDtrackingSector *se // 1. Unpack AliTRDpropagationLayers objects for each stack. // 2. Launch stack tracking. // See AliTRDtrackerV1::Clusters2TracksStack() for details. - // 3. Pack results in the ESD event. - // - - AliTRDpadPlane *pp = 0x0; - - // allocate space for esd tracks in this SM - TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack); - esdTrackList.SetOwner(); - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - Int_t nTimeBins = cal->GetNumberOfTimeBins(); - const Int_t kFindable = Int_t(fRecoParam->GetFindableClusters()*6.*nTimeBins); - - Int_t ntracks = 0; - Int_t nClStack = 0; - for(int istack = 0; istackGetPadPlane((Int_t)(ilayer/nTimeBins), istack); - Double_t stacklength = (pp->GetNrows() - 2) * pp->GetLengthIPad() - + 2 * pp->GetLengthOPad() + 2 * pp->GetLengthRim(); - //Debug - Double_t z0 = fGeom->GetRow0((Int_t)(ilayer/nTimeBins),istack,0); - const AliTRDpropagationLayer ksmLayer(*(sector->GetLayer(ilayer))); - stackLayer[ilayer] = ksmLayer; -#ifdef DEBUG - stackLayer[ilayer].SetDebugStream(fDebugStreamer); -#endif - stackLayer[ilayer].SetRange(z0 - stacklength, stacklength); - stackLayer[ilayer].SetSector(sector->GetSector()); - stackLayer[ilayer].SetStackNr(istack); - stackLayer[ilayer].SetNRows(pp->GetNrows()); - stackLayer[ilayer].SetRecoParam(fRecoParam); - stackLayer[ilayer].BuildIndices(); - nClStack += stackLayer[ilayer].GetNClusters(); - } - //AliInfo(Form("Finish building stack propagation layers. nClusters %d.", nClStack)); - if(nClStack < kFindable) continue; - ntracks += Clusters2TracksStack(&stackLayer[0], &esdTrackList); - } - //AliInfo(Form("Found %d tracks in SM", ntracks)); - - for(int itrack=0; itrackAddTrack((AliESDtrack*)esdTrackList[itrack]); + // 3. Pack results in the ESD event. + // + + // allocate space for esd tracks in this SM + TClonesArray esdTrackList("AliESDtrack", 2*kMaxTracksStack); + esdTrackList.SetOwner(); + + Int_t nTracks = 0; + Int_t nChambers = 0; + AliTRDtrackingChamber **stack = 0x0, *chamber = 0x0; + for(int istack = 0; istackGetNClusters() < fgNTimeBins * fReconstructor->GetRecoParam() ->GetFindableClusters()) continue; + nChambers++; + //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)); + nTracks += Clusters2TracksStack(stack, &esdTrackList); + } + //AliInfo(Form("Found %d tracks in SM %d [%d]\n", nTracks, sector, esd->GetNumberOfTracks())); + + for(int itrack=0; itrackAddTrack((AliESDtrack*)esdTrackList[itrack]); - return ntracks; + // Reset Track and Candidate Number + AliTRDtrackerDebug::SetCandidateNumber(0); + AliTRDtrackerDebug::SetTrackNumber(0); + return nTracks; } //____________________________________________________________________ -Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer - , TClonesArray *esdTrackList) +Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDtrackingChamber **stack, TClonesArray *esdTrackList) { // // Make tracks in one TRD stack. @@ -1063,290 +1843,310 @@ Int_t AliTRDtrackerV1::Clusters2TracksStack(AliTRDstackLayer *layer // 8. Build ESD track and register it to the output list // - AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized - Int_t pars[4]; // MakeSeeds parameters + const AliTRDCalDet *cal = AliTRDcalibDB::Instance()->GetT0Det(); + AliTRDtrackingChamber *chamber = 0x0; + AliTRDseedV1 sseed[kMaxTracksStack*6]; // to be initialized + Int_t pars[4]; // MakeSeeds parameters - //Double_t alpha = AliTRDgeometry::GetAlpha(); - //Double_t shift = .5 * alpha; - Int_t configs[kNConfigs]; - - // Build initial seeding configurations - Double_t quality = BuildSeedingConfigs(layer, configs); -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() > 1) - AliInfo(Form("Plane config %d %d %d Quality %f" - , configs[0], configs[1], configs[2], quality)); -#endif - - // 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; - do{ - // Loop over seeding configurations - ntracks = 0; ntracks1 = 0; - for (Int_t iconf = 0; iconf<3; iconf++) { - pars[0] = configs[iconf]; - pars[1] = layer->GetStackNr(); - pars[2] = ntracks; - ntracks = MakeSeeds(layer, &sseed[6*ntracks], pars); - if(ntracks == kMaxTracksStack) break; - } -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Candidate TRD tracks %d in stack %d iteration %d.", ntracks, pars[1], fSieveSeeding)); -#endif - if(!ntracks) break; - - // Sort the seeds according to their quality - Int_t sort[kMaxTracksStack]; - TMath::Sort(ntracks, fTrackQuality, sort, kTRUE); - - // Initialize number of tracks so far and logic switches - Int_t ntracks0 = esdTrackList->GetEntriesFast(); - Bool_t signedTrack[kMaxTracksStack]; - Bool_t fakeTrack[kMaxTracksStack]; - for (Int_t i=0; iGetLabel(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); - - - // 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; - - - // Build track parameters - AliTRDseedV1 *lseed =&sseed[trackIndex*6]; - Int_t idx = 0; - while(idx<3 && !lseed->IsOK()) { - idx++; - lseed++; - } - Double_t cR = lseed->GetC(); - trackParams[1] = lseed->GetYref(0); - trackParams[2] = 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; - trackParams[0] = lseed->GetX0(); - trackParams[6] = layer[0].GetSector();/* *alpha+shift; // Supermodule*/ - -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() > 1) printf("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1]); - - if(AliTRDReconstructor::StreamLevel() >= 1){ - Int_t sector = layer[0].GetSector(); - 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(); - //for(int ic=0; ic<30; ic++) if(sseed[trackIndex*6+is].GetClusters(ic)) printf("l[%d] tb[%d] cptr[%p]\n", is, ic, sseed[trackIndex*6+is].GetClusters(ic)); - } - //Int_t eventNrInFile = esd->GetEventNumberInFile(); - //AliInfo(Form("Number of clusters %d.", nclusters)); - TTreeSRedirector &cstreamer = *fDebugStreamer; - cstreamer << "Clusters2TracksStack" - << "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] - << "Sector=" << sector - << "Stack=" << pars[1] - << "Label=" << label - << "Label1=" << label1 - << "Label2=" << label2 - << "FakeRatio=" << fakeratio - << "Freq=" << frequency - << "Ncl=" << ncl - << "NLayers=" << nlayers - << "Findable=" << findable - << "NUsed=" << nused - << "\n"; - //???for(int is=0; is<6; is++) delete dseed[is]; - } -#endif - - AliTRDtrackV1 *track = AliTRDtrackerV1::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()); - new ((*esdTrackList)[ntracks0++]) AliESDtrack(esdTrack); - ntracks1++; - } - - jSieve++; - } while(jSieve<5 && candidates); // end track candidates sieve - if(!ntracks1) break; - - // increment counters - ntracks2 += ntracks1; - fSieveSeeding++; - - // Rebuild plane configurations and indices taking only unused clusters into account - quality = BuildSeedingConfigs(layer, configs); - //if(quality < fRecoParam->GetPlaneQualityThreshold()) break; - - for(Int_t il = 0; il < kNPlanes * fTimeBinsPerPlane; il++) layer[il].BuildIndices(fSieveSeeding); - -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() > 1) AliInfo(Form("Sieve level %d Plane config %d %d %d Quality %f", fSieveSeeding, configs[0], configs[1], configs[2], quality)); -#endif - } while(fSieveSeeding<10); // end stack clusters sieve - + //Double_t alpha = AliTRDgeometry::GetAlpha(); + //Double_t shift = .5 * alpha; + Int_t configs[kNConfigs]; + + // Build initial seeding configurations + Double_t quality = BuildSeedingConfigs(stack, configs); + 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(icGetStack((*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(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1) AliInfo(Form("Candidate TRD tracks %d in iteration %d.", ntracks, fSieveSeeding)); + + if(!ntracks) break; + + // Sort the seeds according to their quality + Int_t sort[kMaxTracksStack]; + TMath::Sort(ntracks, fTrackQuality, sort, kTRUE); + + // Initialize number of tracks so far and logic switches + Int_t ntracks0 = esdTrackList->GetEntriesFast(); + Bool_t signedTrack[kMaxTracksStack]; + Bool_t fakeTrack[kMaxTracksStack]; + for (Int_t i=0; i 30){ + //printf("Skip %d nused %d\n", trackIndex, nused); + fakeTrack[trackIndex] = kTRUE; + continue; + } + if (Float_t(nused)/ncl >= .25){ + //printf("Skip %d nused/ncl >= .25\n", trackIndex); + fakeTrack[trackIndex] = kTRUE; + continue; + } + + // Classify tracks + Bool_t skip = kFALSE; + switch(jSieve){ + case 0: + if(nlayers < 6) {skip = kTRUE; break;} + if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;} + break; + + case 1: + if(nlayers < findable){skip = kTRUE; break;} + if(TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -4.){skip = kTRUE; break;} + break; + + case 2: + if ((nlayers == findable) || (nlayers == 6)) { skip = kTRUE; break;} + if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -6.0){skip = kTRUE; break;} + break; + + case 3: + if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) < -5.){skip = kTRUE; break;} + break; + + case 4: + if (nlayers == 3){skip = kTRUE; break;} + //if (TMath::Log(1.E-9+fTrackQuality[trackIndex]) - nused/(nlayers-3.0) < -15.0){skip = kTRUE; break;} + break; + } + if(skip){ + candidates++; + //printf("REJECTED : %d [%d] nlayers %d trackQuality = %e nused %d\n", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused); + continue; + } + signedTrack[trackIndex] = kTRUE; + + + // Sign clusters + AliTRDcluster *cl = 0x0; Int_t clusterIndex = -1; + for (Int_t jLayer = 0; jLayer < kNPlanes; jLayer++) { + Int_t jseed = kNPlanes*trackIndex+jLayer; + if(!sseed[jseed].IsOK()) continue; + if(TMath::Abs(sseed[jseed].GetYfit(1) - sseed[jseed].GetYfit(1)) >= .2) continue; // check this condition with Marian + sseed[jseed].UseClusters(); + if(!cl){ + ic = 0; + while(!(cl = sseed[jseed].GetClusters(ic))) ic++; + clusterIndex = sseed[jseed].GetIndexes(ic); + } + } + if(!cl) continue; + + + // Build track parameters + AliTRDseedV1 *lseed =&sseed[trackIndex*6]; +/* Int_t idx = 0; + while(idx<3 && !lseed->IsOK()) { + idx++; + lseed++; + }*/ + Double_t x = lseed->GetX0();// - 3.5; + trackParams[0] = x; //NEW AB + trackParams[1] = lseed->GetYref(0); // lseed->GetYat(x); + trackParams[2] = lseed->GetZref(0); // lseed->GetZat(x); + trackParams[3] = TMath::Sin(TMath::ATan(lseed->GetYref(1))); + trackParams[4] = lseed->GetZref(1) / TMath::Sqrt(1. + lseed->GetYref(1) * lseed->GetYref(1)); + trackParams[5] = lseed->GetC(); + Int_t ich = 0; while(!(chamber = stack[ich])) ich++; + trackParams[6] = fGeom->GetSector(chamber->GetDetector());/* *alpha+shift; // Supermodule*/ + + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 1){ + AliInfo(Form("Track %d [%d] nlayers %d trackQuality = %e nused %d, yref = %3.3f", itrack, trackIndex, nlayers, fTrackQuality[trackIndex], nused, trackParams[1])); + + Int_t nclusters = 0; + AliTRDseedV1 *dseed[6]; + + // Build track label - what happens if measured data ??? + Int_t labels[1000]; + Int_t outlab[1000]; + Int_t nlab = 0; + + Int_t labelsall[1000]; + Int_t nlabelsall = 0; + Int_t naccepted = 0; + + for (Int_t iLayer = 0; iLayer < kNPlanes; iLayer++) { + Int_t jseed = kNPlanes*trackIndex+iLayer; + dseed[iLayer] = new AliTRDseedV1(sseed[jseed]); + dseed[iLayer]->SetOwner(); + nclusters += sseed[jseed].GetN2(); + if(!sseed[jseed].IsOK()) continue; + for(int ilab=0; ilab<2; ilab++){ + if(sseed[jseed].GetLabels(ilab) < 0) continue; + labels[nlab] = sseed[jseed].GetLabels(ilab); + nlab++; + } + + // Cooking label + for (Int_t itime = 0; itime < fgNTimeBins; itime++) { + if(!sseed[jseed].IsUsable(itime)) continue; + naccepted++; + Int_t tindex = 0, ilab = 0; + while(ilab<3 && (tindex = sseed[jseed].GetClusters(itime)->GetLabel(ilab)) >= 0){ + labelsall[nlabelsall++] = tindex; + ilab++; + } + } + } + Freq(nlab,labels,outlab,kFALSE); + Int_t label = outlab[0]; + Int_t frequency = outlab[1]; + Freq(nlabelsall,labelsall,outlab,kFALSE); + Int_t label1 = outlab[0]; + Int_t label2 = outlab[2]; + Float_t fakeratio = (naccepted - outlab[1]) / Float_t(naccepted); + + //Int_t eventNrInFile = esd->GetEventNumberInFile(); + //AliInfo(Form("Number of clusters %d.", nclusters)); + Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); + Int_t trackNumber = AliTRDtrackerDebug::GetTrackNumber(); + Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); + TTreeSRedirector &cstreamer = *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"; + } + + AliTRDtrackV1 *track = MakeTrack(&sseed[trackIndex*kNPlanes], trackParams); + if(!track){ + AliWarning("Fail to build a TRD Track."); + continue; + } + + //AliInfo("End of MakeTrack()"); + AliESDtrack *esdTrack = new ((*esdTrackList)[ntracks0++]) AliESDtrack(); + esdTrack->UpdateTrackParams(track, AliESDtrack::kTRDout); + esdTrack->SetLabel(track->GetLabel()); + track->UpdateESDtrack(esdTrack); + // write ESD-friends if neccessary + if (fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) > 0){ + AliTRDtrackV1 *calibTrack = new AliTRDtrackV1(*track); + calibTrack->SetOwner(); + esdTrack->AddCalibObject(calibTrack); + } + ntracks1++; + AliTRDtrackerDebug::SetTrackNumber(AliTRDtrackerDebug::GetTrackNumber() + 1); + } + + jSieve++; + } while(jSieve<5 && candidates); // end track candidates sieve + if(!ntracks1) break; + + // 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; //fReconstructor->GetRecoParam() ->GetPlaneQualityThreshold()) break; + + for(Int_t ip = 0; ip < kNPlanes; ip++){ + if(!(chamber = stack[ip])) continue; + chamber->Build(fGeom, cal);//Indices(fSieveSeeding); + } + + 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 + - //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1])); + //AliInfo(Form("Registered TRD tracks %d in stack %d.", ntracks2, pars[1])); - return ntracks2; + return ntracks2; } //___________________________________________________________________ -Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers - , Int_t *configs) +Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDtrackingChamber **stack, Int_t *configs) { // // Assign probabilities to chambers according to their @@ -1372,28 +2172,31 @@ Double_t AliTRDtrackerV1::BuildSeedingConfigs(AliTRDstackLayer *layers // The overall chamber quality is given by the product of this 2 contributions. // - Double_t chamberQA[kNPlanes]; - for(int iplane=0; iplaneGetQuality() : 0.; + } + + Double_t tconfig[kNConfigs]; + Int_t planes[4]; + for(int iconf=0; iconf 1) AliInfo(Form("Making seeds Stack[%d] Config[%d] Tracks[%d]...", istack, config, ntracks)); -#endif - - // Init chambers geometry - Int_t det/*, tbRange[6]*/; // time bins inside the detector geometry - Double_t hL[kNPlanes]; // Tilting angle - Float_t padlength[kNPlanes]; // pad lenghts - AliTRDpadPlane *pp; - for(int il=0; ilGetPadPlane(il, istack); // istack has to be imported - hL[il] = TMath::Tan(-TMath::DegToRad()*pp->GetTiltingAngle()); - padlength[il] = pp->GetLengthIPad(); - det = il; // to be fixed !!!!! - //tbRange[il] = fTimeBinsPerPlane; //Int_t(AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght() * AliTRDCommonParam::Instance()->GetSamplingFrequency()/AliTRDcalibDB::Instance()->GetVdriftDet()->GetValue(det)); - //printf("%d hl[%f] pl[%f] tb[%d]\n", il, hL[il], padlength[il], tbRange[il]); - } - - Double_t cond0[4], cond1[4], cond2[4]; - // make seeding layers (to be moved in Clusters2TracksStack) - AliTRDstackLayer *layer[] = {0x0, 0x0, 0x0, 0x0}; - for(int isl=0; islBuildCond(c[3], cond0, 0); - layer[0]->GetClusters(cond0, index, ncl); - Int_t jcl = 0; - while(jclGetX() - 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); - - Int_t kcl = 0; - while(kclBuildCond(c[1], cond2, 2, theta, phi); - c[2] = layer[2]->GetNearestCluster(cond2); - if(!c[2]) continue; - - //AliInfo("Seeding clusters found. Building seeds ..."); - //for(Int_t i = 0; i < kNSeedPlanes; i++) printf("%i. coordinates: x = %3.3f, y = %3.3f, z = %3.3f\n", i, c[i]->GetX(), c[i]->GetY(), c[i]->GetZ()); - for (Int_t il = 0; il < 6; il++) cseed[il].Reset(); - - fFitter->Reset(); - - fFitter->FitRieman(c, kNSeedPlanes); - - chi2[0] = 0.; chi2[1] = 0.; - AliTRDseedV1 *tseed = 0x0; - for(int iLayer=0; iLayerSetRecoParam(fRecoParam); - tseed->SetPlane(jLayer); - tseed->SetTilt(hL[jLayer]); - tseed->SetPadLength(padlength[jLayer]); - //tseed->SetNTimeBinsRange(tbRange[jLayer]); - tseed->SetX0(layer[iLayer]->GetX());//layers[jLayer*fTimeBinsPerPlane].GetX()); - - tseed->Init(fFitter->GetRiemanFitter()); - // temporary until new AttachClusters() - tseed->SetX0(layers[(jLayer+1)*fTimeBinsPerPlane-1].GetX()); - chi2[0] += tseed->GetChi2Z(c[iLayer]->GetZ()); - chi2[1] += tseed->GetChi2Y(c[iLayer]->GetY()); - } - - Bool_t isFake = kFALSE; - 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; -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() >= 2){ - Float_t yref[4], ycluster[4]; - for(int il=0; il<4; il++){ - tseed = &cseed[planes[il]]; - yref[il] = tseed->GetYref(0); - ycluster[il] = c[il]->GetY(); - } - Float_t threshold = .5;//1./(3. - sLayer); - Int_t ll = c[3]->GetLabel(0); - TTreeSRedirector &cs0 = *fDebugStreamer; - cs0 << "MakeSeeds0" - <<"isFake=" << isFake - <<"label=" << ll - <<"threshold=" << threshold - <<"chi2=" << chi2[1] - <<"yref0="< fRecoParam->GetChi2Z()/*7./(3. - sLayer)*//*iter*/){ - //AliInfo(Form("Failed chi2 filter on chi2Z [%f].", chi2[0])); - continue; - } - if(chi2[1] > fRecoParam->GetChi2Y()/*1./(3. - sLayer)*//*iter*/){ - //AliInfo(Form("Failed chi2 filter on chi2Y [%f].", chi2[1])); - continue; - } - //AliInfo("Passed chi2 filter."); - -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() >= 2){ - Float_t minmax[2] = { -100.0, 100.0 }; - for (Int_t iLayer = 0; iLayer < 4; iLayer++) { - Float_t max = c[iLayer]->GetZ() + cseed[planes[iLayer]].GetPadLength() * 0.5 + 1.0 - cseed[planes[iLayer]].GetZref(0); - if (max < minmax[1]) minmax[1] = max; - Float_t min = c[iLayer]->GetZ()-cseed[planes[iLayer]].GetPadLength() * 0.5 - 1.0 - cseed[planes[iLayer]].GetZref(0); - if (min > minmax[0]) minmax[0] = min; - } - Double_t xpos[4]; - for(Int_t l = 0; l < kNSeedPlanes; l++) xpos[l] = layer[l]->GetX(); - TTreeSRedirector &cstreamer = *fDebugStreamer; - cstreamer << "MakeSeeds1" - << "isFake=" << isFake - << "config=" << config - << "Cl0.=" << c[0] - << "Cl1.=" << c[1] - << "Cl2.=" << c[2] - << "Cl3.=" << c[3] - << "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() - << "Y2exp=" << cond2[0] - << "Z2exp=" << cond2[1] - << "Chi2R=" << chi2[0] - << "Chi2Z=" << chi2[1] - << "Seed0.=" << &cseed[planes[0]] - << "Seed1.=" << &cseed[planes[1]] - << "Seed2.=" << &cseed[planes[2]] - << "Seed3.=" << &cseed[planes[3]] - << "Zmin=" << minmax[0] - << "Zmax=" << minmax[1] - << "\n" ; - } -#endif - // try attaching clusters to tracklets - Int_t nUsedCl = 0; - Int_t nlayers = 0; - for(int iLayer=0; iLayer 25) break; - nlayers++; - } - if(nlayers < kNSeedPlanes){ - //AliInfo("Failed updating all seeds."); - continue; - } - // fit tracklets and cook likelihood - chi2[0] = 0.; chi2[1] = 0.; - fFitter->FitRieman(&cseed[0], &planes[0]); - AliRieman *rim = fFitter->GetRiemanFitter(); - for(int iLayer=0; iLayer<4; iLayer++){ - cseed[planes[iLayer]].Init(rim); - chi2[0] += (Float_t)cseed[planes[iLayer]].GetChi2Z(); - chi2[1] += cseed[planes[iLayer]].GetChi2Y(); - } - Double_t chi2r = chi2[1], chi2z = chi2[0]; - Double_t like = CookLikelihood(&cseed[0], planes, chi2); // to be checked - if (TMath::Log(1.E-9 + like) < fRecoParam->GetTrackLikelihood()){ - //AliInfo(Form("Failed likelihood %f[%e].", TMath::Log(1.E-9 + like), like)); - 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 lextrap[2]; - GetExtrapolationConfig(config, lextrap); - Int_t nusedf = 0; // debug value - for(int iLayer=0; iLayer<2; iLayer++){ - Int_t jLayer = lextrap[iLayer]; - - // prepare extrapolated seed - cseed[jLayer].Reset(); - cseed[jLayer].SetRecoParam(fRecoParam); - cseed[jLayer].SetPlane(jLayer); - cseed[jLayer].SetTilt(hL[jLayer]); - cseed[jLayer].SetX0(layers[(jLayer +1) * fTimeBinsPerPlane-1].GetX()); - cseed[jLayer].SetPadLength(padlength[jLayer]); - //cseed[jLayer].SetNTimeBinsRange(tbRange[jLayer]); - cseed[jLayer].Init(rim); -// AliTRDcluster *cd = FindSeedingCluster(&layers[jLayer*fTimeBinsPerPlane], &cseed[jLayer]); -// if(cd == 0x0) continue; - - // fit extrapolated seed - AliTRDseedV1::FitRiemanTilt(cseed, kTRUE); - if ((jLayer == 0) && !(cseed[1].IsOK())) continue; - if ((jLayer == 5) && !(cseed[4].IsOK())) continue; - AliTRDseedV1 tseed = cseed[jLayer]; - if(!tseed.AttachClustersIter(&layers[jLayer*fTimeBinsPerPlane], 1000.)) continue; - cseed[jLayer] = tseed; - nusedf += cseed[jLayer].GetNUsed(); // debug value - AliTRDseedV1::FitRiemanTilt(cseed, kTRUE); - } - //AliInfo("Extrapolation done."); - - ImproveSeedQuality(layers, cseed); - //AliInfo("Improve seed quality done."); - - nlayers = 0; - Int_t nclusters = 0; - Int_t findable = 0; - for (Int_t iLayer = 0; iLayer < 6; iLayer++) { - if (TMath::Abs(cseed[iLayer].GetYref(0) / cseed[iLayer].GetX0()) < 0.15) findable++; - if (!cseed[iLayer].IsOK()) continue; - nclusters += cseed[iLayer].GetN2(); - nlayers++; - } - if (nlayers < 3){ - //AliInfo("Failed quality check on seeds."); - continue; - } - - // fit full track and cook likelihoods - fFitter->FitRieman(&cseed[0]); - Double_t chi2ZF = 0., chi2RF = 0.; - for(int ilayer=0; ilayer<6; ilayer++){ - cseed[ilayer].Init(fFitter->GetRiemanFitter()); - if (!cseed[ilayer].IsOK()) continue; - //tchi2 = cseed[ilayer].GetChi2Z(); - //printf("layer %d chi2 %e\n", ilayer, tchi2); - chi2ZF += cseed[ilayer].GetChi2Z(); - chi2RF += cseed[ilayer].GetChi2Y(); - } - chi2ZF /= TMath::Max((nlayers - 3.), 1.); - chi2RF /= TMath::Max((nlayers - 3.), 1.); - - // do the final track fitting - fFitter->SetLayers(nlayers); -#ifdef DEBUG - fFitter->SetDebugStream(fDebugStreamer); -#endif - fTrackQuality[ntracks] = fFitter->FitHyperplane(&cseed[0], chi2ZF, GetZ()); - Double_t param[3]; - Double_t chi2[2]; - fFitter->GetHyperplaneFitResults(param); - fFitter->GetHyperplaneFitChi2(chi2); - //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].SetC(param[1]/*cR*/); - cseed[iLayer].SetCC(param[0]/*cC*/); - cseed[iLayer].SetChi2(chi2[0]); - cseed[iLayer].SetChi2Z(chi2ZF); - } - -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() >= 2){ - Double_t curv = (fFitter->GetRiemanFitter())->GetC(); - TTreeSRedirector &cstreamer = *fDebugStreamer; - cstreamer << "MakeSeeds2" - << "C=" << curv - << "Chi2R=" << chi2r - << "Chi2Z=" << chi2z - << "Chi2TR=" << chi2[0] - << "Chi2TC=" << chi2[1] - << "Chi2RF=" << chi2RF - << "Chi2ZF=" << chi2ZF - << "Ncl=" << nclusters - << "Nlayers=" << nlayers - << "NUsedS=" << nUsedCl - << "NUsed=" << nusedf - << "Findable" << findable - << "Like=" << like - << "S0.=" << &cseed[0] - << "S1.=" << &cseed[1] - << "S2.=" << &cseed[2] - << "S3.=" << &cseed[3] - << "S4.=" << &cseed[4] - << "S5.=" << &cseed[5] - << "Label=" << label - << "Freq=" << frequency - << "\n"; - } -#endif - - ntracks++; - 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; + // 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 + Double_t hL[kNPlanes]; // Tilting angle + Float_t padlength[kNPlanes]; // pad lenghts + AliTRDpadPlane *pp = 0x0; + for(int iplane=0; iplaneGetPadPlane(iplane, istack); + hL[iplane] = TMath::Tan(TMath::DegToRad()*pp->GetTiltingAngle()); + padlength[iplane] = pp->GetLengthIPad(); + } + + 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, fReconstructor)) continue; + nlayers++; + } + if(nlayers < 4) return 0; + + + // Start finding seeds + Double_t cond0[4], cond1[4], cond2[4]; + Int_t icl = 0; + while((c[3] = (*fSeedTB[3])[icl++])){ + if(!c[3]) continue; + 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(jclGetX() - c[0]->GetX(); + Double_t theta = (c[3]->GetZ() - c[0]->GetZ())/dx; + Double_t phi = (c[3]->GetY() - c[0]->GetY())/dx; + fSeedTB[1]->BuildCond(c[0], cond1, 1, theta, phi); + fSeedTB[1]->GetClusters(cond1, jndex, mcl); + //printf("Found c[0] candidates 1 %d\n", mcl); + + Int_t kcl = 0; + while(kclBuildCond(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; iLayerSetDetector((*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 25) break; + 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; iLayerGetRecoParam() ->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; iLayerGetStreamLevel(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; + } + } + } + + return ntracks; } //_____________________________________________________________________________ @@ -1815,6 +2597,7 @@ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params) // To be discussed with Marian !! // + Double_t alpha = AliTRDgeometry::GetAlpha(); Double_t shift = AliTRDgeometry::GetAlpha()/2.0; Double_t c[15]; @@ -1825,32 +2608,64 @@ AliTRDtrackV1* AliTRDtrackerV1::MakeTrack(AliTRDseedV1 *seeds, Double_t *params) 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); - AliInfo(Form("N clusters for track %d", nc)); - if (nc < 30) { - delete track; - track = 0x0; - } else { -// track->CookdEdx(); -// track->CookdEdxTimBin(-1); -// CookLabel(track, 0.9); + 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; ipGetEntriesFast()-1); + } + return SetTrack(&track); } - return 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; -//____________________________________________________________________ -void AliTRDtrackerV1::CookLabel(AliKalmanTrack */*pt*/, Float_t /*wrong*/) const -{ - // to be implemented, preferably at the level of TRD tracklet. !!!!!!! + 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; } + //____________________________________________________________________ -void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers - , AliTRDseedV1 *cseed) +Int_t AliTRDtrackerV1::ImproveSeedQuality(AliTRDtrackingChamber **stack, AliTRDseedV1 *cseed) { // // Sort tracklets according to "quality" and try to "improve" the first 4 worst @@ -1868,109 +2683,123 @@ void AliTRDtrackerV1::ImproveSeedQuality(AliTRDstackLayer *layers // tracklet seed such that the seed quality (see AliTRDseed::GetQuality()) // can be maximized. If some optimization is found the old seeds are replaced. // - - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - Int_t nTimeBins = cal->GetNumberOfTimeBins(); - - // make a local working copy - AliTRDseedV1 bseed[6]; - for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer]; - - - Float_t lastquality = 10000.0; - Float_t lastchi2 = 10000.0; - Float_t chi2 = 1000.0; - - for (Int_t iter = 0; iter < 4; iter++) { - Float_t sumquality = 0.0; - Float_t squality[6]; - Int_t sortindexes[6]; - - for (Int_t jLayer = 0; jLayer < 6; jLayer++) { - squality[jLayer] = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : -1.; - sumquality +=squality[jLayer]; - } - if ((sumquality >= lastquality) || (chi2 > lastchi2)) break; - - - lastquality = sumquality; - lastchi2 = chi2; - if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer]; - - - TMath::Sort(6, squality, sortindexes, kFALSE); - for (Int_t jLayer = 5; jLayer > 1; jLayer--) { - Int_t bLayer = sortindexes[jLayer]; - bseed[bLayer].AttachClustersIter(&layers[bLayer*nTimeBins], squality[bLayer], kTRUE); - } - - chi2 = AliTRDseedV1::FitRiemanTilt(bseed,kTRUE); - } // Loop: iter + // debug level: 7 + // + + // make a local working copy + AliTRDtrackingChamber *chamber = 0x0; + AliTRDseedV1 bseed[6]; + Int_t nLayers = 0; + for (Int_t jLayer = 0; jLayer < 6; jLayer++) bseed[jLayer] = cseed[jLayer]; + + Float_t lastquality = 10000.0; + Float_t lastchi2 = 10000.0; + Float_t chi2 = 1000.0; + + for (Int_t iter = 0; iter < 4; iter++) { + Float_t sumquality = 0.0; + Float_t squality[6]; + Int_t sortindexes[6]; + + for (Int_t jLayer = 0; jLayer < 6; jLayer++) { + squality[jLayer] = bseed[jLayer].IsOK() ? bseed[jLayer].GetQuality(kTRUE) : 1000.; + sumquality += squality[jLayer]; + } + if ((sumquality >= lastquality) || (chi2 > lastchi2)) break; + + nLayers = 0; + lastquality = sumquality; + lastchi2 = chi2; + if (iter > 0) for (Int_t jLayer = 0; jLayer < 6; jLayer++) cseed[jLayer] = bseed[jLayer]; + + TMath::Sort(6, squality, sortindexes, kFALSE); + 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); + if(bseed[bLayer].IsOK()) nLayers++; + } + + chi2 = FitTiltedRieman(bseed, kTRUE); + if(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 7){ + Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); + Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); + TLinearFitter *tiltedRieman = GetTiltedRiemanFitter(); + TTreeSRedirector &cstreamer = *fgDebugStreamer; + cstreamer << "ImproveSeedQuality" + << "EventNumber=" << eventNumber + << "CandidateNumber=" << candidateNumber + << "Iteration=" << iter + << "S0.=" << &bseed[0] + << "S1.=" << &bseed[1] + << "S2.=" << &bseed[2] + << "S3.=" << &bseed[3] + << "S4.=" << &bseed[4] + << "S5.=" << &bseed[5] + << "FitterT.=" << tiltedRieman + << "\n"; + } + } // Loop: iter + + // we are sure that at least 2 tracklets are OK ! + return nLayers+2; } -//____________________________________________________________________ -Double_t AliTRDtrackerV1::MakeSeedingPlanes(AliTRDstackLayer *layers) -{ +//_________________________________________________________________________ +Double_t AliTRDtrackerV1::CalculateTrackLikelihood(AliTRDseedV1 *tracklets, Double_t *chi2){ // - // Calculate plane quality for seeding. - // + // Calculates the Track Likelihood value. This parameter serves as main quality criterion for + // the track selection + // The likelihood value containes: + // - The chi2 values from the both fitters and the chi2 values in z-direction from a linear fit + // - The Sum of the Parameter |slope_ref - slope_fit|/Sigma of the tracklets + // For all Parameters an exponential dependency is used // - // Parameters : - // layers : Array of propagation layers for this plane. + // Parameters: - Array of tracklets (AliTRDseedV1) related to the track candidate + // - Array of chi2 values: + // * Non-Constrained Tilted Riemann fit + // * Vertex-Constrained Tilted Riemann fit + // * z-Direction from Linear fit + // Output: - The calculated track likelihood // - // Output : - // plane quality factor for seeding - // - // Detailed description + // debug level 2 // - // The quality of the plane for seeding is higher if: - // 1. the average timebin population is closer to an integer number - // 2. the distribution of clusters/timebin is closer to a uniform distribution. - // - the slope of the first derivative of a parabolic fit is small or - // - the slope of a linear fit is small - // - - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - Int_t nTimeBins = cal->GetNumberOfTimeBins(); - -// Double_t x; -// TLinearFitter fitter(1, "pol1"); -// fitter.ClearPoints(); - Int_t ncl = 0; - Int_t nused = 0; - Int_t nClLayer; - for(int itb=0; itbIsUsed()) nused++; - } - - // calculate the deviation of the mean number of clusters from the - // closest integer values - Float_t nclMed = float(ncl-nused)/nTimeBins; - Int_t ncli = Int_t(nclMed); - Float_t nclDev = TMath::Abs(nclMed - TMath::Max(ncli, 1)); - nclDev -= (nclDev>.5) && ncli ? 0. : 1.; - /*Double_t quality = */ return TMath::Exp(2.*nclDev); - -// // get slope of the derivative -// if(!fitter.Eval()) return quality; -// fitter.PrintResults(3); -// Double_t a = fitter.GetParameter(1); -// -// printf("ncl_dev(%f) a(%f)\n", ncl_dev, a); -// return quality*TMath::Exp(-a); + + Double_t sumdaf = 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()); + nLayers++; + } + sumdaf /= Float_t (nLayers - 2.0); + + Double_t likeChi2Z = TMath::Exp(-chi2[2] * 0.14); // Chi2Z + 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(fReconstructor->GetStreamLevel(AliTRDReconstructor::kTracker) >= 2){ + Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber(); + Int_t candidateNumber = AliTRDtrackerDebug::GetCandidateNumber(); + TTreeSRedirector &cstreamer = *fgDebugStreamer; + cstreamer << "CalculateTrackLikelihood0" + << "EventNumber=" << eventNumber + << "CandidateNumber=" << candidateNumber + << "LikeChi2Z=" << likeChi2Z + << "LikeChi2TR=" << likeChi2TR + << "LikeChi2TC=" << likeChi2TC + << "LikeAF=" << likeAF + << "TrackLikelihood=" << trackLikelihood + << "\n"; + } + + return trackLikelihood; } //____________________________________________________________________ -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. @@ -1995,435 +2824,68 @@ Double_t AliTRDtrackerV1::CookLikelihood(AliTRDseedV1 *cseed // The distributions for each type of probabilities are given below as of // (date). They have to be checked to assure consistency of estimation. // - - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - Int_t nTimeBins = cal->GetNumberOfTimeBins(); - // ratio of the total number of clusters/track which are expected to be found by the tracker. - Float_t fgFindable = fRecoParam->GetFindableClusters(); - - - Int_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); - Double_t likechi2y = 0.0000000001; - if (chi2[1] < 0.5) likechi2y += TMath::Exp(-TMath::Sqrt(chi2[1]) * 7.73); - Double_t likechi2z = TMath::Exp(-chi2[0] * 0.088) / TMath::Exp(-chi2[0] * 0.019); - Int_t enc = Int_t(fgFindable*4.*nTimeBins); // Expected Number Of Clusters, normally 72 - Double_t likeN = TMath::Exp(-(enc - nclusters) * 0.19); - - Double_t like = likea * likechi2y * likechi2z * likeN; - -#ifdef DEBUG - //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){ - TTreeSRedirector &cstreamer = *fDebugStreamer; - cstreamer << "CookLikelihood" - << "sumda=" << sumda - << "chi0=" << chi2[0] - << "chi1=" << chi2[1] - << "likea=" << likea - << "likechi2y=" << likechi2y - << "likechi2z=" << likechi2z - << "nclusters=" << nclusters - << "likeN=" << likeN - << "like=" << like - << "\n"; - } -#endif - - return like; -} - -//___________________________________________________________________ -void AliTRDtrackerV1::GetMeanCLStack(AliTRDstackLayer *layers - , Int_t *planes - , Double_t *params) -{ - // - // Determines the Mean number of clusters per layer. - // Needed to determine good Seeding Layers - // - // Parameters: - // - Array of AliTRDstackLayers - // - Container for the params - // - // Detailed description - // - // Two Iterations: - // In the first Iteration the mean is calculted using all layers. - // After this, all layers outside the 1-sigma-region are rejected. - // Then the mean value and the standard-deviation are calculted a second - // time in order to select all layers in the 1-sigma-region as good-candidates. - // - - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - Int_t nTimeBins = cal->GetNumberOfTimeBins(); - - Float_t mean = 0, stdev = 0; - Double_t ncl[kNTimeBins*kNSeedPlanes], mcl[kNTimeBins*kNSeedPlanes]; - Int_t position = 0; - memset(ncl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes); - memset(mcl, 0, sizeof(Int_t)*kNTimeBins*kNSeedPlanes); - Int_t nused = 0; - for(Int_t ipl = 0; ipl < kNSeedPlanes; ipl++){ - for(Int_t ils = 0; ils < nTimeBins; ils++){ - position = planes[ipl]*nTimeBins + ils; - ncl[ipl * nTimeBins + ils] = layers[position].GetNClusters(); - nused = 0; - for(Int_t icl = 0; icl < ncl[ipl * nTimeBins + ils]; icl++) - if((layers[position].GetCluster(icl))->IsUsed()) nused++; - ncl[ipl * nTimeBins + ils] -= nused; - } - } - // Declaration of quartils: - //Double_t qvals[3] = {0.0, 0.0, 0.0}; - //Double_t qprop[3] = {0.16667, 0.5, 0.83333}; - // Iterations - Int_t counter; - Double_t *array; - Int_t *limit; - Int_t nLayers = nTimeBins * kNSeedPlanes; - for(Int_t iter = 0; iter < 2; iter++){ - array = (iter == 0) ? &ncl[0] : &mcl[0]; - limit = (iter == 0) ? &nLayers : &counter; - counter = 0; - if(iter == 1){ - for(Int_t i = 0; i < nTimeBins *kNSeedPlanes; i++){ - if((ncl[i] > mean + stdev) || (ncl[i] < mean - stdev)) continue; // Outside 1-sigma region -// if((ncl[i] > qvals[2]) || (ncl[i] < qvals[0])) continue; // Outside 1-sigma region - if(ncl[i] == 0) continue; // 0-Layers also rejected - mcl[counter] = ncl[i]; - counter++; - } - } - if(*limit == 0) break; - printf("Limit = %d\n", *limit); - //using quartils instead of mean and RMS -// TMath::Quantiles(*limit,3,array,qvals,qprop,kFALSE); - mean = TMath::Median(*limit, array, 0x0); - stdev = TMath::RMS(*limit, array); - } -// printf("Quantiles: 0.16667 = %3.3f, 0.5 = %3.3f, 0.83333 = %3.3f\n", qvals[0],qvals[1],qvals[2]); -// memcpy(params,qvals,sizeof(Double_t)*3); - params[1] = (Double_t)TMath::Nint(mean); - params[0] = (Double_t)TMath::Nint(mean - stdev); - params[2] = (Double_t)TMath::Nint(mean + stdev); -} + // ratio of the total number of clusters/track which are expected to be found by the tracker. + const AliTRDrecoParam *fRecoPars = fReconstructor->GetRecoParam(); + + 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)); + } + nclusters *= .25; + + Double_t likea = TMath::Exp(-sumda * fRecoPars->GetPhiCut()); + Double_t likechi2y = 0.0000000001; + 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(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" + << "EventNumber=" << eventNumber + << "CandidateNumber=" << candidateNumber + << "tracklet0.=" << &cseed[0] + << "tracklet1.=" << &cseed[1] + << "tracklet2.=" << &cseed[2] + << "tracklet3.=" << &cseed[3] + << "tracklet4.=" << &cseed[4] + << "tracklet5.=" << &cseed[5] + << "sumda=" << sumda + << "chi2y=" << chi2y + << "chi2z=" << chi2z + << "likea=" << likea + << "likechi2y=" << likechi2y + << "likechi2z=" << likechi2z + << "nclusters=" << nclusters + << "likeN=" << likeN + << "like=" << like + << "meanncls=" << mean_ncls + << "\n"; + } -//___________________________________________________________________ -Int_t AliTRDtrackerV1::GetSeedingLayers(AliTRDstackLayer *layers - , Double_t *params) -{ - // - // Algorithm to find optimal seeding layer - // Layers inside one sigma region (given by Quantiles) are sorted - // according to their difference. - // All layers outside are sorted according t - // - // Parameters: - // - Array of AliTRDstackLayers (in the current plane !!!) - // - Container for the Indices of the seeding Layer candidates - // - // Output: - // - Number of Layers inside the 1-sigma-region - // - // The optimal seeding layer should contain the mean number of - // custers in the layers in one chamber. - // - - //printf("Params: %3.3f, %3.3f, %3.3f\n", params[0], params[1], params[2]); - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - const Int_t kMaxClustersLayer = AliTRDstackLayer::kMaxClustersLayer; - Int_t nTimeBins = cal->GetNumberOfTimeBins(); - Int_t ncl[kNTimeBins], indices[kNTimeBins], bins[kMaxClustersLayer]; - memset(ncl, 0, sizeof(Int_t)*kNTimeBins); - memset(indices, 0, sizeof(Int_t)*kNTimeBins); - memset(bins, 0, sizeof(Int_t)*kMaxClustersLayer); - Int_t nused = 0; - for(Int_t ils = 0; ils < nTimeBins; ils++){ - ncl[ils] = layers[ils].GetNClusters(); - nused = 0; - for(Int_t icl = 0; icl < ncl[ils]; icl++) - if((layers[ils].GetCluster(icl))->IsUsed()) nused++; - ncl[ils] -= nused; - } - - Float_t mean = params[1]; - for(Int_t ils = 0; ils < nTimeBins; ils++){ - memmove(indices + bins[ncl[ils]+1] + 1, indices + bins[ncl[ils]+1], sizeof(Int_t)*(nTimeBins - ils)); - indices[bins[ncl[ils]+1]] = ils; - for(Int_t i = ncl[ils]+1; i < kMaxClustersLayer; i++) - bins[i]++; - } - - //for(Int_t i = 0; i < nTimeBins; i++) printf("Bin %d = %d\n", i, bins[i]); - Int_t sbin = -1; - Int_t nElements; - Int_t position = 0; - TRandom *r = new TRandom(); - Int_t iter = 0; - while(1){ - while(sbin < (Int_t)params[0] || sbin > (Int_t)params[2]){ - // Randomly selecting one bin - sbin = (Int_t)r->Poisson(mean); - } - printf("Bin = %d\n",sbin); - //Randomly selecting one Layer in the bin - nElements = bins[sbin + 1] - bins[sbin]; - printf("nElements = %d\n", nElements); - if(iter == 5){ - position = (Int_t)(gRandom->Rndm()*(nTimeBins-1)); - break; - } - else if(nElements==0){ - iter++; - continue; - } - position = (Int_t)(gRandom->Rndm()*(nElements-1)) + bins[sbin]; - break; - } - delete r; - return indices[position]; + return like; } -//____________________________________________________________________ -AliTRDcluster *AliTRDtrackerV1::FindSeedingCluster(AliTRDstackLayer *layers - , AliTRDseedV1/*AliRieman*/ *reference) -{ - // - // Finds a seeding Cluster for the extrapolation chamber. - // - // The seeding cluster should be as close as possible to the assumed - // track which is represented by a Rieman fit. - // Therefore the selecting criterion is the minimum distance between - // the best fitting cluster and the Reference which is derived from - // the AliTRDseed. Because all layers are assumed to be equally good - // a linear search is performed. - // - // Imput parameters: - layers: array of AliTRDstackLayers (in one chamber!!!) - // - sfit: the reference - // - // Output: - the best seeding cluster - // - - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - Int_t nTimeBins = cal->GetNumberOfTimeBins(); - - // distances as squared distances - Int_t index = 0; - Float_t ypos = 0.0, zpos = 0.0, distance = 0.0, nearestDistance =100000.0; - ypos = reference->GetYref(0); - zpos = reference->GetZref(0); - AliTRDcluster *currentBest = 0x0, *temp = 0x0; - for(Int_t ils = 0; ils < nTimeBins; ils++){ - // Reference positions -// ypos = reference->GetYat(layers[ils].GetX()); -// zpos = reference->GetZat(layers[ils].GetX()); - index = layers[ils].SearchNearestCluster(ypos, zpos, fRecoParam->GetRoad2y(), fRecoParam->GetRoad2z()); - if(index == -1) continue; - temp = layers[ils].GetCluster(index); - if(!temp) continue; - distance = (temp->GetY() - ypos) * (temp->GetY() - ypos) + (temp->GetZ() - zpos) * (temp->GetZ() - zpos); - if(distance < nearestDistance){ - nearestDistance = distance; - currentBest = temp; - } - } - return currentBest; -} -//____________________________________________________________________ -AliTRDstackLayer *AliTRDtrackerV1::MakeSeedingLayer(AliTRDstackLayer *layers - , Int_t plane) -{ - // - // Creates a seeding layer - // - - // constants - const Int_t kMaxRows = 16; - const Int_t kMaxCols = 144; - const Int_t kMaxPads = 2304; - - // Get the calculation - AliTRDcalibDB *cal = AliTRDcalibDB::Instance(); - Int_t nTimeBins = cal->GetNumberOfTimeBins(); - - // Get the geometrical data of the chamber - AliTRDpadPlane *pp = fGeom->GetPadPlane(plane, layers[0].GetStackNr()); - Int_t nCols = pp->GetNcols(); - Float_t ymin = TMath::Min(pp->GetCol0(), pp->GetColEnd()); - Float_t ymax = TMath::Max(pp->GetCol0(), pp->GetColEnd()); - Float_t zmin = TMath::Min(pp->GetRow0(), pp->GetRowEnd()); - Float_t zmax = TMath::Max(pp->GetRow0(), pp->GetRowEnd()); - Int_t nRows = pp->GetNrows(); - Float_t binlength = (ymax - ymin)/nCols; - //AliInfo(Form("ymin(%f) ymax(%f) zmin(%f) zmax(%f) nRows(%d) binlength(%f)", ymin, ymax, zmin, zmax, nRows, binlength)); - - // Fill the histogram - Int_t arrpos; - Float_t ypos; - Int_t irow, nClusters; - Int_t *histogram[kMaxRows]; // 2D-Histogram - Int_t hvals[kMaxPads]; memset(hvals, 0, sizeof(Int_t)*kMaxPads); - Float_t *sigmas[kMaxRows]; - Float_t svals[kMaxPads]; memset(svals, 0, sizeof(Float_t)*kMaxPads); - AliTRDcluster *c = 0x0; - for(Int_t irs = 0; irs < kMaxRows; irs++){ - histogram[irs] = &hvals[irs*kMaxCols]; - sigmas[irs] = &svals[irs*kMaxCols]; - } - for(Int_t iTime = 0; iTime < nTimeBins; iTime++){ - nClusters = layers[iTime].GetNClusters(); - for(Int_t incl = 0; incl < nClusters; incl++){ - c = layers[iTime].GetCluster(incl); - ypos = c->GetY(); - if(ypos > ymax && ypos < ymin) continue; - irow = pp->GetPadRowNumber(c->GetZ()); // Zbin - if(irow < 0)continue; - arrpos = static_cast((ypos - ymin)/binlength); - if(ypos == ymax) arrpos = nCols - 1; - histogram[irow][arrpos]++; - sigmas[irow][arrpos] += c->GetSigmaZ2(); - } - } - -// Now I have everything in the histogram, do the selection -// printf("Starting the analysis\n"); - //Int_t nPads = nCols * nRows; - // This is what we are interested in: The center of gravity of the best candidates - Float_t cogyvals[kMaxPads]; memset(cogyvals, 0, sizeof(Float_t)*kMaxPads); - Float_t cogzvals[kMaxPads]; memset(cogzvals, 0, sizeof(Float_t)*kMaxPads); - Float_t *cogy[kMaxRows]; - Float_t *cogz[kMaxRows]; - // Lookup-Table storing coordinates according ti the bins - Float_t yLengths[kMaxCols]; - Float_t zLengths[kMaxRows]; - for(Int_t icnt = 0; icnt < nCols; icnt++){ - yLengths[icnt] = pp->GetColPos(nCols - 1 - icnt) + binlength/2; - } - for(Int_t icnt = 0; icnt < nRows; icnt++){ - zLengths[icnt] = pp->GetRowPos(icnt) - pp->GetRowSize(icnt)/2; - } - - // A bitfield is used to mask the pads as usable - Short_t mask[kMaxCols]; memset(mask, 0 ,sizeof(Short_t) * kMaxCols);//bool mvals[kMaxPads]; - for(UChar_t icount = 0; icount < nRows; icount++){ - cogy[icount] = &cogyvals[icount*kMaxCols]; - cogz[icount] = &cogzvals[icount*kMaxCols]; - } - // In this array the array position of the best candidates will be stored - Int_t cand[kMaxTracksStack]; - Float_t sigcands[kMaxTracksStack]; - - // helper variables - Int_t indices[kMaxPads]; memset(indices, 0, sizeof(Int_t)*kMaxPads); - Int_t nCandidates = 0; - Float_t norm, cogv; - // histogram filled -> Select best bins - TMath::Sort(kMaxPads, hvals, indices); // bins storing a 0 should not matter - // Set Threshold - Int_t maximum = hvals[indices[0]]; // best - Int_t threshold = static_cast(maximum * fRecoParam->GetFindableClusters()); - Int_t col, row, lower, lower1, upper, upper1; - for(Int_t ib = 0; ib < kMaxPads; ib++){ - if(nCandidates >= kMaxTracksStack){ - AliWarning(Form("Number of seed candidates %d exceeded maximum allowed per stack %d", nCandidates, kMaxTracksStack)); - break; - } - // Positions - row = indices[ib]/nCols; - col = indices[ib]%nCols; - // here will be the threshold condition: - if((mask[col] & (1 << row)) != 0) continue; // Pad is masked: continue - if(histogram[row][col] < TMath::Max(threshold, 1)){ // of course at least one cluster is needed - break; // number of clusters below threshold: break; - } - // passing: Mark the neighbors - lower = TMath::Max(col - 1, 0); upper = TMath::Min(col + 2, nCols); - lower1 = TMath::Max(row - 1, 0); upper1 = TMath::Min(row + 2, nCols); - for(Int_t ic = lower; ic < upper; ++ic) - for(Int_t ir = lower1; ir < upper1; ++ir){ - if(ic == col && ir == row) continue; - mask[ic] |= (1 << ir); - } - // Storing the position in an array - // testing for neigboring - cogv = 0; - norm = 0; - lower = TMath::Max(col - 1,0); - upper = TMath::Min(col + 2, nCols); - for(Int_t inb = lower; inb < upper; ++inb){ - cogv += yLengths[inb] * histogram[row][inb]; - norm += histogram[row][inb]; - } - cogy[row][col] = cogv / norm; - cogv = 0; norm = 0; - lower = TMath::Max(row - 1, 0); - upper = TMath::Min(row + 2, nRows); - for(Int_t inb = lower; inb < upper; ++inb){ - cogv += zLengths[inb] * histogram[inb][col]; - norm += histogram[inb][col]; - } - cogz[row][col] = cogv / norm; - // passed the filter - cand[nCandidates] = row*kMaxCols + col; // store the position of a passig candidate into an Array - sigcands[nCandidates] = sigmas[row][col] / histogram[row][col]; // never be a floating point exeption - // Analysis output - nCandidates++; - } - AliTRDstackLayer *fakeLayer = new AliTRDstackLayer(layers[0].GetZ0(), layers[0].GetDZ0(), layers[0].GetStackNr()); - fakeLayer->SetX((TMath::Abs(layers[nTimeBins-1].GetX() + layers[0].GetX()))/2); - fakeLayer->SetSector(layers[0].GetSector()); - AliTRDcluster **fakeClusters = 0x0; - UInt_t *fakeIndices = 0x0; - if(nCandidates){ - fakeClusters = new AliTRDcluster*[nCandidates]; - fakeIndices = new UInt_t[nCandidates]; - UInt_t fakeIndex = 0; - for(Int_t ican = 0; ican < nCandidates; ican++){ - fakeClusters[ican] = new AliTRDcluster(); - fakeClusters[ican]->SetX(fakeLayer->GetX()); - fakeClusters[ican]->SetY(cogyvals[cand[ican]]); - fakeClusters[ican]->SetZ(cogzvals[cand[ican]]); - fakeClusters[ican]->SetSigmaZ2(sigcands[ican]); - fakeIndices[ican] = fakeIndex++;// fantasy number - } - } - fakeLayer->SetRecoParam(fRecoParam); - fakeLayer->SetClustersArray(fakeClusters, nCandidates); - fakeLayer->SetIndexArray(fakeIndices); - fakeLayer->SetNRows(nRows); - fakeLayer->BuildIndices(); - //fakeLayer->PrintClusters(); - -#ifdef DEBUG - if(AliTRDReconstructor::StreamLevel() >= 3){ - TMatrixD hist(nRows, nCols); - for(Int_t i = 0; i < nRows; i++) - for(Int_t j = 0; j < nCols; j++) - hist(i,j) = histogram[i][j]; - TTreeSRedirector &cstreamer = *fDebugStreamer; - cstreamer << "MakeSeedingLayer" - << "Iteration=" << fSieveSeeding - << "plane=" << plane - << "ymin=" << ymin - << "ymax=" << ymax - << "zmin=" << zmin - << "zmax=" << zmax - << "L.=" << fakeLayer - << "Histogram.=" << &hist - << "\n"; - } -#endif - return fakeLayer; -} //____________________________________________________________________ void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4]) @@ -2472,98 +2934,98 @@ void AliTRDtrackerV1::GetSeedingConfig(Int_t iconfig, Int_t planes[4]) //End_Html // - switch(iconfig){ - case 0: // 5432 TQ 0 - planes[0] = 2; - planes[1] = 3; - planes[2] = 4; - planes[3] = 5; - break; - case 1: // 4321 TQ 0 - planes[0] = 1; - planes[1] = 2; - planes[2] = 3; - planes[3] = 4; - break; - case 2: // 3210 TQ 0 - planes[0] = 0; - planes[1] = 1; - planes[2] = 2; - planes[3] = 3; - break; - case 3: // 5321 TQ 1 - planes[0] = 1; - planes[1] = 2; - planes[2] = 3; - planes[3] = 5; - break; - case 4: // 4210 TQ 1 - planes[0] = 0; - planes[1] = 1; - planes[2] = 2; - planes[3] = 4; - break; - case 5: // 5431 TQ 1 - planes[0] = 1; - planes[1] = 3; - planes[2] = 4; - planes[3] = 5; - break; - case 6: // 4320 TQ 1 - planes[0] = 0; - planes[1] = 2; - planes[2] = 3; - planes[3] = 4; - break; - case 7: // 5430 TQ 2 - planes[0] = 0; - planes[1] = 3; - planes[2] = 4; - planes[3] = 5; - break; - case 8: // 5210 TQ 2 - planes[0] = 0; - planes[1] = 1; - planes[2] = 2; - planes[3] = 5; - break; - case 9: // 5421 TQ 3 - planes[0] = 1; - planes[1] = 2; - planes[2] = 4; - planes[3] = 5; - break; - case 10: // 4310 TQ 3 - planes[0] = 0; - planes[1] = 1; - planes[2] = 3; - planes[3] = 4; - break; - case 11: // 5410 TQ 4 - planes[0] = 0; - planes[1] = 1; - planes[2] = 4; - planes[3] = 5; - break; - case 12: // 5420 TQ 5 - planes[0] = 0; - planes[1] = 2; - planes[2] = 4; - planes[3] = 5; - break; - case 13: // 5320 TQ 5 - planes[0] = 0; - planes[1] = 2; - planes[2] = 3; - planes[3] = 5; - break; - case 14: // 5310 TQ 5 - planes[0] = 0; - planes[1] = 1; - planes[2] = 3; - planes[3] = 5; - break; - } + switch(iconfig){ + case 0: // 5432 TQ 0 + planes[0] = 2; + planes[1] = 3; + planes[2] = 4; + planes[3] = 5; + break; + case 1: // 4321 TQ 0 + planes[0] = 1; + planes[1] = 2; + planes[2] = 3; + planes[3] = 4; + break; + case 2: // 3210 TQ 0 + planes[0] = 0; + planes[1] = 1; + planes[2] = 2; + planes[3] = 3; + break; + case 3: // 5321 TQ 1 + planes[0] = 1; + planes[1] = 2; + planes[2] = 3; + planes[3] = 5; + break; + case 4: // 4210 TQ 1 + planes[0] = 0; + planes[1] = 1; + planes[2] = 2; + planes[3] = 4; + break; + case 5: // 5431 TQ 1 + planes[0] = 1; + planes[1] = 3; + planes[2] = 4; + planes[3] = 5; + break; + case 6: // 4320 TQ 1 + planes[0] = 0; + planes[1] = 2; + planes[2] = 3; + planes[3] = 4; + break; + case 7: // 5430 TQ 2 + planes[0] = 0; + planes[1] = 3; + planes[2] = 4; + planes[3] = 5; + break; + case 8: // 5210 TQ 2 + planes[0] = 0; + planes[1] = 1; + planes[2] = 2; + planes[3] = 5; + break; + case 9: // 5421 TQ 3 + planes[0] = 1; + planes[1] = 2; + planes[2] = 4; + planes[3] = 5; + break; + case 10: // 4310 TQ 3 + planes[0] = 0; + planes[1] = 1; + planes[2] = 3; + planes[3] = 4; + break; + case 11: // 5410 TQ 4 + planes[0] = 0; + planes[1] = 1; + planes[2] = 4; + planes[3] = 5; + break; + case 12: // 5420 TQ 5 + planes[0] = 0; + planes[1] = 2; + planes[2] = 4; + planes[3] = 5; + break; + case 13: // 5320 TQ 5 + planes[0] = 0; + planes[1] = 2; + planes[2] = 3; + planes[3] = 5; + break; + case 14: // 5310 TQ 5 + planes[0] = 0; + planes[1] = 1; + planes[2] = 3; + planes[3] = 5; + break; + } } //____________________________________________________________________ @@ -2582,66 +3044,344 @@ void AliTRDtrackerV1::GetExtrapolationConfig(Int_t iconfig, Int_t planes[2]) // Detailed description // - switch(iconfig){ - case 0: // 5432 TQ 0 - planes[0] = 1; - planes[1] = 0; - break; - case 1: // 4321 TQ 0 - planes[0] = 5; - planes[1] = 0; - break; - case 2: // 3210 TQ 0 - planes[0] = 4; - planes[1] = 5; - break; - case 3: // 5321 TQ 1 - planes[0] = 4; - planes[1] = 0; - break; - case 4: // 4210 TQ 1 - planes[0] = 5; - planes[1] = 3; - break; - case 5: // 5431 TQ 1 - planes[0] = 2; - planes[1] = 0; - break; - case 6: // 4320 TQ 1 - planes[0] = 5; - planes[1] = 1; - break; - case 7: // 5430 TQ 2 - planes[0] = 2; - planes[1] = 1; - break; - case 8: // 5210 TQ 2 - planes[0] = 4; - planes[1] = 3; - break; - case 9: // 5421 TQ 3 - planes[0] = 3; - planes[1] = 0; - break; - case 10: // 4310 TQ 3 - planes[0] = 5; - planes[1] = 2; - break; - case 11: // 5410 TQ 4 - planes[0] = 3; - planes[1] = 2; - break; - case 12: // 5420 TQ 5 - planes[0] = 3; - planes[1] = 1; - break; - case 13: // 5320 TQ 5 - planes[0] = 4; - planes[1] = 1; - break; - case 14: // 5310 TQ 5 - planes[0] = 4; - planes[1] = 2; - break; - } + switch(iconfig){ + case 0: // 5432 TQ 0 + planes[0] = 1; + planes[1] = 0; + break; + case 1: // 4321 TQ 0 + planes[0] = 5; + planes[1] = 0; + break; + case 2: // 3210 TQ 0 + planes[0] = 4; + planes[1] = 5; + break; + case 3: // 5321 TQ 1 + planes[0] = 4; + planes[1] = 0; + break; + case 4: // 4210 TQ 1 + planes[0] = 5; + planes[1] = 3; + break; + case 5: // 5431 TQ 1 + planes[0] = 2; + planes[1] = 0; + break; + case 6: // 4320 TQ 1 + planes[0] = 5; + planes[1] = 1; + break; + case 7: // 5430 TQ 2 + planes[0] = 2; + planes[1] = 1; + break; + case 8: // 5210 TQ 2 + planes[0] = 4; + planes[1] = 3; + break; + case 9: // 5421 TQ 3 + planes[0] = 3; + planes[1] = 0; + break; + case 10: // 4310 TQ 3 + planes[0] = 5; + planes[1] = 2; + break; + case 11: // 5410 TQ 4 + planes[0] = 3; + planes[1] = 2; + break; + case 12: // 5420 TQ 5 + planes[0] = 3; + planes[1] = 1; + break; + case 13: // 5320 TQ 5 + planes[0] = 4; + planes[1] = 1; + break; + case 14: // 5310 TQ 5 + planes[0] = 4; + planes[1] = 2; + break; + } +} + +//____________________________________________________________________ +AliCluster* AliTRDtrackerV1::GetCluster(Int_t idx) const +{ + Int_t ncls = fClusters->GetEntriesFast(); + 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; +} + +//____________________________________________________________________ +AliKalmanTrack* AliTRDtrackerV1::GetTrack(Int_t idx) const +{ + Int_t ntrk = fTracks->GetEntriesFast(); + return idx >= 0 && idx < ntrk ? (AliKalmanTrack*)fTracks->UncheckedAt(idx) : 0x0; +} + +//____________________________________________________________________ +Float_t AliTRDtrackerV1::CalculateReferenceX(AliTRDseedV1 *tracklets){ + // + // Calculates the reference x-position for the tilted Rieman fit defined as middle + // of the stack (middle between layers 2 and 3). For the calculation all the tracklets + // are taken into account + // + // Parameters: - Array of tracklets(AliTRDseedV1) + // + // Output: - The reference x-position(Float_t) + // + Int_t nDistances = 0; + Float_t meanDistance = 0.; + Int_t startIndex = 5; + for(Int_t il =5; il > 0; il--){ + if(tracklets[il].IsOK() && tracklets[il -1].IsOK()){ + Float_t xdiff = tracklets[il].GetX0() - tracklets[il -1].GetX0(); + meanDistance += xdiff; + nDistances++; + } + if(tracklets[il].IsOK()) startIndex = il; + } + if(tracklets[0].IsOK()) startIndex = 0; + if(!nDistances){ + // We should normally never get here + Float_t xpos[2]; memset(xpos, 0, sizeof(Float_t) * 2); + Int_t iok = 0, idiff = 0; + // This attempt is worse and should be avoided: + // check for two chambers which are OK and repeat this without taking the mean value + // Strategy avoids a division by 0; + for(Int_t il = 5; il >= 0; il--){ + if(tracklets[il].IsOK()){ + xpos[iok] = tracklets[il].GetX0(); + iok++; + startIndex = il; + } + if(iok) idiff++; // to get the right difference; + if(iok > 1) break; + } + if(iok > 1){ + meanDistance = (xpos[0] - xpos[1])/idiff; + } + else{ + // we have do not even have 2 layers which are OK? The we do not need to fit at all + return 331.; + } + } + else{ + meanDistance /= nDistances; + } + return tracklets[startIndex].GetX0() + (2.5 - startIndex) * meanDistance - 0.5 * (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick()); +} + +//_____________________________________________________________________________ +Int_t AliTRDtrackerV1::Freq(Int_t n, const Int_t *inlist + , Int_t *outlist, Bool_t down) +{ + // + // Sort eleements according occurancy + // The size of output array has is 2*n + // + + if (n <= 0) { + return 0; + } + + Int_t *sindexS = new Int_t[n]; // Temporary array for sorting + Int_t *sindexF = new Int_t[2*n]; + for (Int_t i = 0; i < n; i++) { + sindexF[i] = 0; + } + + TMath::Sort(n,inlist,sindexS,down); + + Int_t last = inlist[sindexS[0]]; + Int_t val = last; + sindexF[0] = 1; + sindexF[0+n] = last; + Int_t countPos = 0; + + // 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; + +} + + +//____________________________________________________________________ +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 +{ + // Chi2 definition on y-direction + + 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));// /tracklets[ipl].GetSigmaY(); + chi2 += distLayer * distLayer; + } + return chi2; +} + +//____________________________________________________________________ +void AliTRDtrackerV1::ResetSeedTB() +{ +// reset buffer for seeding time bin layers. If the time bin +// layers are not allocated this function allocates them + + for(Int_t isl=0; islClear(); + } +} + +//_____________________________________________________________________________ +Float_t AliTRDtrackerV1::GetChi2Z(AliTRDseedV1 *tracklets) const +{ + // 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)); // /(tracklets[ipl].GetPadLength()/2); + chi2 += distLayer * distLayer; + } + return chi2; +} + +/////////////////////////////////////////////////////// +// // +// Resources of class AliTRDLeastSquare // +// // +/////////////////////////////////////////////////////// + +//_____________________________________________________________________________ +AliTRDtrackerV1::AliTRDLeastSquare::AliTRDLeastSquare(){ + // + // Constructor of the nested class AliTRDtrackFitterLeastSquare + // + memset(fParams, 0, sizeof(Double_t) * 2); + memset(fSums, 0, sizeof(Double_t) * 5); + memset(fCovarianceMatrix, 0, sizeof(Double_t) * 3); + +} + +//_____________________________________________________________________________ +void AliTRDtrackerV1::AliTRDLeastSquare::AddPoint(Double_t *x, Double_t y, Double_t sigmaY){ + // + // Adding Point to the fitter + // + Double_t weight = 1/(sigmaY * sigmaY); + Double_t &xpt = *x; + // printf("Adding point x = %f, y = %f, sigma = %f\n", xpt, y, sigmaY); + fSums[0] += weight; + fSums[1] += weight * xpt; + fSums[2] += weight * y; + fSums[3] += weight * xpt * y; + fSums[4] += weight * xpt * xpt; + fSums[5] += weight * y * y; +} + +//_____________________________________________________________________________ +void AliTRDtrackerV1::AliTRDLeastSquare::RemovePoint(Double_t *x, Double_t y, Double_t sigmaY){ + // + // Remove Point from the sample + // + Double_t weight = 1/(sigmaY * sigmaY); + Double_t &xpt = *x; + fSums[0] -= weight; + fSums[1] -= weight * xpt; + fSums[2] -= weight * y; + fSums[3] -= weight * xpt * y; + fSums[4] -= weight * xpt * xpt; + fSums[5] -= weight * y * y; +} + +//_____________________________________________________________________________ +void AliTRDtrackerV1::AliTRDLeastSquare::Eval(){ + // + // Evaluation of the fit: + // Calculation of the parameters + // Calculation of the covariance matrix + // + + 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); + fParams[0] = (fSums[2] * fSums[4] - fSums[1] * fSums[3])/ denominator; + fParams[1] = (fSums[0] * fSums[3] - fSums[1] * fSums[2]) / denominator; + // printf("fParams[0] = %f, fParams[1] = %f\n", fParams[0], fParams[1]); + + // Covariance matrix + fCovarianceMatrix[0] = fSums[4] - fSums[1] * fSums[1] / fSums[0]; + fCovarianceMatrix[1] = fSums[5] - fSums[2] * fSums[2] / fSums[0]; + fCovarianceMatrix[2] = fSums[3] - fSums[1] * fSums[2] / fSums[0]; +} + +//_____________________________________________________________________________ +Double_t AliTRDtrackerV1::AliTRDLeastSquare::GetFunctionValue(Double_t *xpos) const { + // + // Returns the Function value of the fitted function at a given x-position + // + return fParams[0] + fParams[1] * (*xpos); +} + +//_____________________________________________________________________________ +void AliTRDtrackerV1::AliTRDLeastSquare::GetCovarianceMatrix(Double_t *storage) const { + // + // Copies the values of the covariance matrix into the storage + // + memcpy(storage, fCovarianceMatrix, sizeof(Double_t) * 3); } +