From 4f1c04d342d28b759082cb6f3b25df4b80e7aa2d Mon Sep 17 00:00:00 2001 From: hristov Date: Wed, 26 Jan 2005 17:03:17 +0000 Subject: [PATCH] Removing wide clusters from the reconstruction (M.Ivanov) --- TRD/AliTRDcluster.cxx | 3 +- TRD/AliTRDcluster.h | 16 +- TRD/AliTRDtrack.cxx | 127 ++++++++-- TRD/AliTRDtrack.h | 25 +- TRD/AliTRDtracker.cxx | 569 +++++++++++++++++++++++++++--------------- TRD/AliTRDtracker.h | 1 + 6 files changed, 505 insertions(+), 236 deletions(-) diff --git a/TRD/AliTRDcluster.cxx b/TRD/AliTRDcluster.cxx index bc8b8c72a91..2a0e378176c 100644 --- a/TRD/AliTRDcluster.cxx +++ b/TRD/AliTRDcluster.cxx @@ -50,7 +50,7 @@ ClassImp(AliTRDcluster) fSigmaY2 = 0.2; fSigmaZ2 = 5.; - + fNPads =0; } //_____________________________________________________________________________ @@ -72,6 +72,7 @@ AliTRDcluster::AliTRDcluster(const AliTRDcluster &c):AliCluster() fDetector = c.GetDetector(); fTimeBin = c.GetLocalTimeBin(); fQ = c.GetQ(); + fNPads = c.fNPads; } diff --git a/TRD/AliTRDcluster.h b/TRD/AliTRDcluster.h index 451df071615..189bc20dbbe 100644 --- a/TRD/AliTRDcluster.h +++ b/TRD/AliTRDcluster.h @@ -14,7 +14,7 @@ class AliTRDcluster : public AliCluster { public: - AliTRDcluster() : AliCluster() { fQ=0; fTimeBin=0; fDetector=0; } + AliTRDcluster() : AliCluster() { fQ=0; fTimeBin=0; fDetector=0; fNPads=0; } AliTRDcluster(const AliTRDcluster &c); AliTRDcluster(const AliTRDrecPoint &p); @@ -39,12 +39,12 @@ class AliTRDcluster : public AliCluster { Int_t GetLocalTimeBin() const { return fTimeBin; } Float_t GetQ() const { return fQ; } - void Set2pad() { SetBit(k2pad); } - void Set3pad() { SetBit(k3pad); } - void Set4pad() { SetBit(k4pad); } - void Set5pad() { SetBit(k5pad); } - void SetLarge() { SetBit(kLarge); } - + void Set2pad() { SetBit(k2pad); fNPads=2; } + void Set3pad() { SetBit(k3pad); fNPads=3; } + void Set4pad() { SetBit(k4pad); fNPads=4; } + void Set5pad() { SetBit(k5pad); fNPads=5; } + void SetLarge() { SetBit(kLarge);fNPads=6; } + Int_t GetNPads() const {return fNPads;} protected: enum { @@ -58,7 +58,7 @@ class AliTRDcluster : public AliCluster { Int_t fDetector; // TRD detector number Int_t fTimeBin; // Time bin number within the detector Float_t fQ; // amplitude - + Int_t fNPads; // number of pads in cluster ClassDef(AliTRDcluster,1) // Cluster for the TRD }; diff --git a/TRD/AliTRDtrack.cxx b/TRD/AliTRDtrack.cxx index 27a88d8292c..96ec89c6cb8 100644 --- a/TRD/AliTRDtrack.cxx +++ b/TRD/AliTRDtrack.cxx @@ -63,6 +63,12 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, fNWrong = 0; fNRotate = 0; fStopped = 0; + fNCross =0; + fNLast =0; + fChi2Last=0; + fNExpected=0; + fNExpectedLast=0; + fNdedx=0; Double_t q = TMath::Abs(c->GetQ()); Double_t s = fX*fC - fE, t=fT; if(s*s < 1) q *= TMath::Sqrt((1-s*s)/(1+t*t)); @@ -75,7 +81,6 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, fIndex[i] = 0; fIndexBackup[i] = 0; //bacup indexes MI } - fNCross =0; fBackupTrack =0; } @@ -99,9 +104,17 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) : AliKalmanTrack(t) { fNWrong = t.fNWrong; fNRotate = t.fNRotate; fStopped = t.fStopped; + fNCross = t.fNCross; + fNExpected = t.fNExpected; + fNExpectedLast = t.fNExpectedLast; + fNdedx = t.fNdedx; + fNLast = t.fNLast; + fChi2Last = t.fChi2Last; + fBackupTrack =0; fAlpha=t.fAlpha; fX=t.fX; + fY=t.fY; fZ=t.fZ; fE=t.fE; fT=t.fT; fC=t.fC; fCyy=t.fCyy; @@ -124,8 +137,6 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) : AliKalmanTrack(t) { fIndex[i] = 0; fIndexBackup[i] = 0; //MI backup indexes } - fNCross =0; - fBackupTrack =0; } //_____________________________________________________________________________ @@ -150,6 +161,13 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack& t, Double_t alpha) fNWrong = 0; fNRotate = 0; fStopped = 0; + fNExpected=0; + fNExpectedLast=0; + fNdedx =0; + fNCross =0; + fNLast =0; + fChi2Last =0; + fBackupTrack =0; fAlpha = alpha; if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi(); @@ -188,8 +206,7 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack& t, Double_t alpha) fIndex[i] = 0; fIndexBackup[i] = 0; // MI backup indexes } - fNCross =0; - fBackupTrack =0; + } //_____________________________________________________________________________ AliTRDtrack::AliTRDtrack(const AliESDtrack& t) @@ -217,6 +234,13 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t) fNWrong = 0; fStopped = 0; fNRotate = 0; + fNExpected =0; + fNExpectedLast=0; + fNdedx = 0; + fNCross =0; + fNLast =0; + fChi2Last =0; + fBackupTrack =0; fAlpha = t.GetAlpha(); if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi(); @@ -259,8 +283,6 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t) fdQdl[i] = 0; // fIndex[i] = 0; //MI store indexes } - fNCross =0; - fBackupTrack =0; if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return; StartTimeIntegral(); @@ -279,6 +301,22 @@ AliTRDtrack::~AliTRDtrack() } + +Float_t AliTRDtrack::StatusForTOF() +{ + Int_t status=0; + if (GetNumberOfClusters()<20) return 0; // + if (fN>110&&fChi2/(Float_t(fN))<3) return 3; //gold + if (fNLast>30&&fChi2Last/(Float_t(fNLast))<3) return 3; //gold + if (fNLast>20&&fChi2Last/(Float_t(fNLast))<2) return 3; //gold + if (fNLast/(fNExpectedLast+3.)>0.8 && fChi2Last/Float_t(fNLast)<5&&fNLast>20) return 2; //silber + if (fNLast>5 &&((fNLast+1.)/(fNExpectedLast+1.))>0.8&&fChi2Last/(fNLast-5.)<6) return 1; + // + + return status; +} + + //____________________________________________________________________________ void AliTRDtrack::GetExternalParameters(Double_t& xr, Double_t x[5]) const { // @@ -347,7 +385,12 @@ void AliTRDtrack::CookdEdx(Double_t low, Double_t up) { //----------------------------------------------------------------- Int_t i; - Int_t nc=GetNumberOfClusters(); + //Int_t nc=GetNumberOfClusters(); + Int_t nc=fNdedx; + if (nc<10) { + SetdEdx(0); + return; + } Float_t sorted[kMAX_CLUSTERS_PER_TRACK]; for (i=0; i < nc; i++) { @@ -610,7 +653,7 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index, } //_____________________________________________________________________________ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index, Double_t h01, - Int_t plane) + Int_t /*plane*/) { // Assignes found cluster to the track and updates track information @@ -625,8 +668,15 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index } Double_t tangent = TMath::Sqrt(tangent2); if ((fC*fX-fE)<0) tangent*=-1; - Double_t correction = 0*plane; - Double_t errang = tangent2*0.04; // + // Double_t correction = 0*plane; + Double_t errang = tangent2*0.04; // + Double_t errsys =0.025*0.025*20; //systematic error part + Float_t extend =1; + if (c->GetNPads()==4) extend=2; + //if (c->GetNPads()==5) extend=3; + //if (c->GetNPads()==6) extend=3; + //if (c->GetQ()<15) return 1; + /* if (corrector!=0){ //if (0){ @@ -640,9 +690,9 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index } */ // - Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.); + // Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.); - Double_t r00=c->GetSigmaY2() +errang, r01=0., r11=c->GetSigmaZ2()*10000.; + Double_t r00=(c->GetSigmaY2() +errang+errsys)*extend, r01=0., r11=c->GetSigmaZ2()*10000.; r00+=fCyy; r01+=fCzy; r11+=fCzz; Double_t det=r00*r11 - r01*r01; Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det; @@ -670,19 +720,29 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index fC = cur; } else { + Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12); + Double_t xu_factor = 1000.; // empirical factor set by C.Xu // in the first tilt version dy=c->GetY() - fY; dz=c->GetZ() - fZ; - dy=dy+h01*dz+correction; + //dy=dy+h01*dz+correction; + + Double_t tiltdz = dz; + if (TMath::Abs(tiltdz)>padlength/2.) { + tiltdz = TMath::Sign(padlength/2,dz); + } + // dy=dy+h01*dz; + dy=dy+h01*tiltdz; + Double_t add=0; if (TMath::Abs(dz)>padlength/2.){ //Double_t dy2 = c->GetY() - fY; //Double_t sign = (dz>0) ? -1.: 1.; //dy2-=h01*sign*padlength/2.; //dy = dy2; - add =1.; + add =1; } - Double_t s00 = c->GetSigmaY2()+errang+add; // error pad + Double_t s00 = (c->GetSigmaY2()+errang)*extend+errsys+add; // error pad Double_t s11 = c->GetSigmaZ2()*xu_factor; // error pad-row // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00; @@ -838,6 +898,8 @@ Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) con chi2 = (dy*dy)/r00; } else { + Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12); + // r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2(); r00+=fCyy; r01+=fCzy; r11+=fCzz; @@ -849,7 +911,12 @@ Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) con } Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01; Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ; - dy=dy+h01*dz; + Double_t tiltdz = dz; + if (TMath::Abs(tiltdz)>padlength/2.) { + tiltdz = TMath::Sign(padlength/2,dz); + } + // dy=dy+h01*dz; + dy=dy+h01*tiltdz; chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; } @@ -911,17 +978,39 @@ void AliTRDtrack::ResetCovariance(Float_t mult) { // fCyy*=mult; - fCzy*=0.; fCzz*=mult; + fCzy*=0.; fCzz*=1.; fCey*=0.; fCez*=0.; fCee*=mult; - fCty*=0.; fCtz*=0.; fCte*=0.; fCtt*=mult; + fCty*=0.; fCtz*=0.; fCte*=0.; fCtt*=1.; fCcy*=0.; fCcz*=0.; fCce*=0.; fCct*=0.; fCcc*=mult; } + + + void AliTRDtrack::MakeBackupTrack() { // // if (fBackupTrack) delete fBackupTrack; fBackupTrack = new AliTRDtrack(*this); + +} + +Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z){ + // + // Find prolongation at given x + // return 0 if not exist + + Double_t c1=fC*fX - fE; + if (TMath::Abs(c1)>1.) return 0; + Double_t r1=TMath::Sqrt(1.- c1*c1); + Double_t c2=fC*xk - fE; + if (TMath::Abs(c2)>1.) return 0; + Double_t r2=TMath::Sqrt(1.- c2*c2); + y =fY + (xk-fX)*(c1+c2)/(r1+r2); + z =fZ + (xk-fX)*(c1+c2)/(c1*r2 + c2*r1)*fT; + + return 1; + } diff --git a/TRD/AliTRDtrack.h b/TRD/AliTRDtrack.h index 1cac1566296..69a6bb08c6b 100644 --- a/TRD/AliTRDtrack.h +++ b/TRD/AliTRDtrack.h @@ -20,7 +20,7 @@ const unsigned kMAX_CLUSTERS_PER_TRACK=210; class AliTRDtrack : public AliKalmanTrack { // Represents reconstructed TRD track - + friend class AliTRDtracker; public: AliTRDtrack():AliKalmanTrack(){fBackupTrack=0;} @@ -31,8 +31,8 @@ public: AliTRDtrack(const AliESDtrack& t); ~AliTRDtrack(); Int_t Compare(const TObject *o) const; - void CookdEdx(Double_t low=0.05, Double_t up=0.70); - + void CookdEdx(Double_t low=0.05, Double_t up=0.55); + Float_t StatusForTOF(); Double_t GetAlpha() const {return fAlpha;} Int_t GetSector() const { //if (fabs(fAlpha) < AliTRDgeometry::GetAlpha()/2) return 0; @@ -74,8 +74,6 @@ public: Double_t GetZ() const {return fZ;} UInt_t * GetBackupIndexes() {return fIndexBackup;} UInt_t * GetIndexes() {return fIndex;} - - Double_t GetYat(Double_t xk) const { //----------------------------------------------------------------- // This function calculates the Y-coordinate of a track at the plane x=xk. @@ -85,6 +83,8 @@ public: Double_t c2=fC*xk - fE, r2=TMath::Sqrt(1.- c2*c2); return fY + (xk-fX)*(c1+c2)/(r1+r2); } + Int_t GetProlongation(Double_t xk, Double_t &y, Double_t &z); + void SetStop(Bool_t stop) {fStopped=stop;} Bool_t GetStop() const {return fStopped;} @@ -100,9 +100,15 @@ public: void SetLikelihoodElectron(Float_t l) { fLhElectron = l; }; void SetSampledEdx(Float_t q, Int_t i) { + Double_t s=GetSnp(), t=GetTgl(); + q*= TMath::Sqrt((1-s*s)/(1+t*t)); + fdQdl[i]=q; + } + void SetSampledEdx(Float_t q) { Double_t s=GetSnp(), t=GetTgl(); q*= TMath::Sqrt((1-s*s)/(1+t*t)); - fdQdl[i]=q; + fdQdl[fNdedx]=q; + fNdedx++; } void SetSeedLabel(Int_t lab) { fSeedLab=lab; } @@ -118,7 +124,7 @@ public: Int_t GetNWrong() const {return fNWrong;} Int_t GetNRotate() const {return fNRotate;} Int_t GetNCross() const {return fNCross;} - void IncCross() {fNCross++;} + void IncCross() {fNCross++; if (fBackupTrack) fBackupTrack->IncCross();} AliTRDtrack * GetBackupTrack(){return fBackupTrack;} void MakeBackupTrack(); // @@ -156,6 +162,11 @@ protected: Int_t fNWrong; // number of wrong clusters Int_t fNRotate; // number of rotation Int_t fNCross; // number of the cross materials + Int_t fNExpected; //expected number of cluster + Int_t fNLast; //number of clusters in last 2 layers + Int_t fNExpectedLast; //number of expected clusters on last 2 layers + Int_t fNdedx; //number of clusters for dEdx measurment + Float_t fChi2Last; //chi2 in the last 2 layers AliTRDtrack * fBackupTrack; //! backup track ClassDef(AliTRDtrack,2) // TRD reconstructed tracks diff --git a/TRD/AliTRDtracker.cxx b/TRD/AliTRDtracker.cxx index bee93c684a3..dc221e493d5 100644 --- a/TRD/AliTRDtracker.cxx +++ b/TRD/AliTRDtracker.cxx @@ -453,10 +453,23 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { Int_t found=0; Float_t foundMin = 20; - Int_t n = event->GetNumberOfTracks(); + // + //Sort tracks + Float_t *quality =new Float_t[n]; + Int_t *index =new Int_t[n]; for (Int_t i=0; iGetTrack(i); + Double_t covariance[15]; + seed->GetExternalCovariance(covariance); + quality[i] = covariance[0]+covariance[2]; + } + TMath::Sort(n,quality,index,kFALSE); + // + for (Int_t i=0; iGetTrack(i); + AliESDtrack* seed=event->GetTrack(index[i]); + ULong_t status=seed->GetStatus(); if ( (status & AliESDtrack::kTPCout ) == 0 ) continue; if ( (status & AliESDtrack::kTRDout) != 0 ) continue; @@ -466,7 +479,7 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { track->SetSeedLabel(lbl); seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); //make backup fNseeds++; - Float_t p4 = track->GetC(); + Float_t p4 = track->GetC(); // Int_t expectedClr = FollowBackProlongation(*track); /* @@ -482,54 +495,57 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { delete track2; } */ - if (TMath::Abs(track->GetC()-p4)/TMath::Abs(p4)>0.2) { - delete track; - continue; //too big change of curvature - to be checked - } - - Int_t foundClr = track->GetNumberOfClusters(); - if (foundClr >= foundMin) { - track->CookdEdx(0.,1.); - CookdEdxTimBin(*track); - - CookLabel(track, 1-fgkLabelFraction); - if(track->GetChi2()/track->GetNumberOfClusters()<6) { // sign only gold tracks - UseClusters(track); - } - Bool_t isGold = kFALSE; - - if (track->GetChi2()/track->GetNumberOfClusters()<5) { //full gold track - seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); - isGold = kTRUE; - } - if (!isGold && track->GetNCross()==0&&track->GetChi2()/track->GetNumberOfClusters()<7){ //almost gold track - seed->UpdateTrackParams(track, 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); + if (TMath::Abs(track->GetC()-p4)/TMath::Abs(p4)<0.2 || TMath::Abs(track->GetPt())>0.8 ) { + // + //make backup for back propagation + // + Int_t foundClr = track->GetNumberOfClusters(); + if (foundClr >= foundMin) { + track->CookdEdx(); + CookLabel(track, 1-fgkLabelFraction); + if(track->GetChi2()/track->GetNumberOfClusters()<4) { // sign only gold tracks + if (seed->GetKinkIndex(0)==0&&TMath::Abs(track->GetPt())<1.5 ) UseClusters(track); + } + Bool_t isGold = kFALSE; + + if (track->GetChi2()/track->GetNumberOfClusters()<5) { //full gold track + seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); + isGold = kTRUE; + } + if (!isGold && track->GetNCross()==0&&track->GetChi2()/track->GetNumberOfClusters()<7){ //almost gold track + seed->UpdateTrackParams(track, 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; + } + } } } - else{ - delete track; - continue; - } - + // + //Propagation to the TOF (I.Belikov) if (track->GetStop()==kFALSE){ - + Double_t xtof=371.; Double_t c2=track->GetC()*xtof - track->GetEta(); - if (TMath::Abs(c2)>=0.85) { + if (TMath::Abs(c2)>=0.99) { delete track; continue; } - Double_t xTOF0 = 371. ; + Double_t xTOF0 = 365. ; PropagateToOuterPlane(*track,xTOF0); + // + //energy losses taken to the account - check one more time + c2=track->GetC()*xtof - track->GetEta(); + if (TMath::Abs(c2)>=0.99) { + delete track; + continue; + } + // Double_t ymax=xtof*TMath::Tan(0.5*AliTRDgeometry::GetAlpha()); Double_t y=track->GetYat(xtof); @@ -566,7 +582,7 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { found++; } } - + seed->SetTRDQuality(track->StatusForTOF()); delete track; //End of propagation to the TOF @@ -580,7 +596,9 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { cerr<<"Number of back propagated TRD tracks: "<Clear(); fNseeds=0; - + delete [] index; + delete [] quality; + return 0; } @@ -601,35 +619,42 @@ Int_t AliTRDtracker::RefitInward(AliESD* event) Int_t nseed = 0; Int_t found = 0; Int_t innerTB = fTrSec[0]->GetInnerTimeBin(); + AliTRDtrack seed2; Int_t n = event->GetNumberOfTracks(); for (Int_t i=0; iGetTrack(i); - AliTRDtrack* seed2 = new AliTRDtrack(*seed); - if (seed2->GetX()<270){ - seed->UpdateTrackParams(seed2, AliESDtrack::kTRDbackup); // backup TPC track - only update - delete seed2; + new(&seed2) AliTRDtrack(*seed); + if (seed2.GetX()<270){ + seed->UpdateTrackParams(&seed2, AliESDtrack::kTRDbackup); // backup TPC track - only update continue; } ULong_t status=seed->GetStatus(); if ( (status & AliESDtrack::kTRDout ) == 0 ) { - delete seed2; continue; } if ( (status & AliESDtrack::kTRDin) != 0 ) { - delete seed2; continue; } nseed++; - seed2->ResetCovariance(5.); - AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha()); + if (1/seed2.Get1Pt()>5.&& seed2.GetX()>260.) { + Double_t oldx = seed2.GetX(); + seed2.PropagateTo(500.); + seed2.ResetCovariance(1.); + seed2.PropagateTo(oldx); + } + else{ + seed2.ResetCovariance(5.); + } + + AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha()); + UInt_t * indexes2 = seed2.GetIndexes(); for (Int_t i=0;iSetPIDsignals(seed2->GetPIDsignals(i),i); - pt->SetPIDTimBin(seed2->GetPIDTimBin(i),i); + pt->SetPIDsignals(seed2.GetPIDsignals(i),i); + pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i); } - UInt_t * indexes2 = seed2->GetIndexes(); UInt_t * indexes3 = pt->GetBackupIndexes(); for (Int_t i=0;i<200;i++) { if (indexes2[i]==0) break; @@ -638,20 +663,6 @@ Int_t AliTRDtracker::RefitInward(AliESD* event) //AliTRDtrack *pt = seed2; AliTRDtrack &t=*pt; FollowProlongation(t, innerTB); - /* - if (t.GetNumberOfClusters()GetTRDclusters(indexes3)*0.5){ - // debug - why we dont go back? - AliTRDtrack *pt2 = new AliTRDtrack(*seed2,seed2->GetAlpha()); - UInt_t * indexes2 = seed2->GetIndexes(); - UInt_t * indexes3 = pt2->GetBackupIndexes(); - for (Int_t i=0;i<200;i++) { - if (indexes2[i]==0) break; - indexes3[i] = indexes2[i]; - } - FollowProlongation(*pt2, innerTB); - delete pt2; - } - */ if (t.GetNumberOfClusters() >= foundMin) { // UseClusters(&t); //CookLabel(pt, 1-fgkLabelFraction); @@ -683,8 +694,6 @@ Int_t AliTRDtracker::RefitInward(AliESD* event) } delete pt2; } - - delete seed2; delete pt; } @@ -811,7 +820,7 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf) wZwindow = TMath::Sqrt(2.25 * 12 * sz2); // Find the closest correct cluster for debugging purposes - if (timeBin) { + if (timeBin&&fVocal) { Float_t minDY = 1000000; for (Int_t i=0; iGetY() > y+road) break; if (c->IsUsed() > 0) continue; @@ -875,8 +884,8 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf) } if(!cl) { - - for (Int_t i=timeBin.Find(y-road); iGetY() > y+road) break; @@ -891,13 +900,16 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf) cl=c; index=timeBin.GetIndex(i); } - } + } + */ if (cl) { + wYclosest = cl->GetY(); wZclosest = cl->GetZ(); Double_t h01 = GetTiltFactor(cl); - t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); + if (cl->GetNPads()<5) + t.SetSampledEdx(cl->GetQ()/dx); //printf("Track position\t%f\t%f\t%f\n",t.GetX(),t.GetY(),t.GetZ()); //printf("Cluster position\t%d\t%f\t%f\n",cl->GetLocalTimeBin(),cl->GetY(),cl->GetZ()); Int_t det = cl->GetDetector(); @@ -913,32 +925,6 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf) //if (tryAgain==0) break; tryAgain--; } - - /* - if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) { - - printf(" %f", wIndex); //1 - printf(" %f", wTB); //2 - printf(" %f", wYrt); //3 - printf(" %f", wYclosest); //4 - printf(" %f", wYcorrect); //5 - printf(" %f", wYwindow); //6 - printf(" %f", wZrt); //7 - printf(" %f", wZclosest); //8 - printf(" %f", wZcorrect); //9 - printf(" %f", wZwindow); //10 - printf(" %f", wPx); //11 - printf(" %f", wPy); //12 - printf(" %f", wPz); //13 - printf(" %f", wSigmaC2*1000000); //14 - printf(" %f", wSigmaTgl2*1000); //15 - printf(" %f", wSigmaY2); //16 - // printf(" %f", wSigmaZ2); //17 - printf(" %f", wChi2); //17 - printf(" %f", wC); //18 - printf("\n"); - } - */ } } } @@ -960,8 +946,6 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) Float_t wIndex, wTB, wChi2; Float_t wYrt, wYclosest, wYcorrect, wYwindow; Float_t wZrt, wZclosest, wZcorrect, wZwindow; - Float_t wPx, wPy, wPz, wC; - Double_t px, py, pz; Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2; Int_t trackIndex = t.GetLabel(); @@ -971,6 +955,9 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) TVector2::Phi_0_2pi(alpha); Int_t s; + + Int_t clusters[1000]; + for (Int_t i=0;i<1000;i++) clusters[i]=-1; Int_t outerTB = fTrSec[0]->GetOuterTimeBin(); Double_t radLength, rho, x, dx, y, ymax = 0, z; @@ -1009,20 +996,23 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) Float_t angle = t.GetAlpha(); // MI - if rotation - we go through the material if (!AdjustSector(&t)) break; Int_t cross = kFALSE; - + Int_t crosz = kFALSE; if (TMath::Abs(angle - t.GetAlpha())>0.000001) cross = kTRUE; //better to stop track Int_t currentzone = fTrSec[s]->GetLayer(nr)->GetZone(z); - if (currentzone==-10) cross = kTRUE; // we are in the frame + if (currentzone==-10) {cross = kTRUE,crosz=kTRUE;} // we are in the frame if (currentzone>-10){ // layer knows where we are if (zone==-10) zone = currentzone; - if (zone!=currentzone) cross=kTRUE; + if (zone!=currentzone) { + cross=kTRUE; + crosz=kTRUE; + } } + if (TMath::Abs(t.GetSnp())>0.8 && t.GetBackupTrack()==0) t.MakeBackupTrack(); if (cross) { + if (t.GetNCross()==0 && t.GetBackupTrack()==0) t.MakeBackupTrack(); t.IncCross(); - if (t.GetNCross()==1) t.MakeBackupTrack(); - if (t.GetNCross()>2) break; + if (t.GetNCross()>4) break; } - // // s = t.GetSector(); @@ -1050,12 +1040,16 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) // now propagate to the middle plane of the next time bin fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster); - + if (crosz){ + rho = 1000*2.7; radLength = 24.01; //TEMPORARY - aluminium in between z - will be detected using GeoModeler in future versions + } x = fTrSec[s]->GetLayer(nr+1)->GetX(); - if(!t.PropagateTo(x,radLength,rho)) break; + if(!t.PropagateTo(x,radLength,rho)) break; if (!AdjustSector(&t)) break; s = t.GetSector(); - if(!t.PropagateTo(x,radLength,rho)) break; + if(!t.PropagateTo(x,radLength,rho)) break; + + if (TMath::Abs(t.GetSnp())>0.95) break; y = t.GetY(); z = t.GetZ(); @@ -1065,8 +1059,13 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) // trackIndex, nr+1, lookForCluster); if(lookForCluster) { - expectedNumberOfClusters++; +// if (clusters[nr]==-1) { +// FindClusters(s,nr,nr+30,&t,clusters); +// } + expectedNumberOfClusters++; + t.fNExpected++; + if (t.fX>345) t.fNExpectedLast++; wIndex = (Float_t) t.GetLabel(); wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex(); @@ -1080,11 +1079,6 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) wYrt = (Float_t) y; wZrt = (Float_t) z; wYwindow = (Float_t) road; - t.GetPxPyPz(px,py,pz); - wPx = (Float_t) px; - wPy = (Float_t) py; - wPz = (Float_t) pz; - wC = (Float_t) t.GetC(); wSigmaC2 = (Float_t) t.GetSigmaC2(); wSigmaTgl2 = (Float_t) t.GetSigmaTgl2(); wSigmaY2 = (Float_t) t.GetSigmaY2(); @@ -1121,7 +1115,7 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) wZwindow = TMath::Sqrt(2.25 * 12 * sz2); // Find the closest correct cluster for debugging purposes - if (timeBin) { + if (timeBin&&fVocal) { minDY = 1000000; for (Int_t i=0; i0) { + index = clusters[nr+1]; + cl = (AliTRDcluster*)GetCluster(index); + Double_t h01 = GetTiltFactor(cl); + maxChi2=t.GetPredictedChi2(cl,h01); + } + */ - for (Int_t i=timeBin.Find(y-road); iGetY() > y+road) break; - if (c->IsUsed() > 0) continue; - if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue; - - Double_t h01 = GetTiltFactor(c); - chi2=t.GetPredictedChi2(c,h01); - - if (chi2 > maxChi2) continue; - maxChi2=chi2; - cl=c; - index=timeBin.GetIndex(i); + if (!cl){ + Int_t maxn = timeBin; + for (Int_t i=timeBin.Find(y-road); iGetY() > y+road) break; + if (c->IsUsed() > 0) continue; + if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue; + Double_t h01 = GetTiltFactor(c); + chi2=t.GetPredictedChi2(c,h01); + + if (chi2 > maxChi2) continue; + maxChi2=chi2; + cl=c; + index=timeBin.GetIndex(i); + //check is correct - if((c->GetLabel(0) != trackIndex) && - (c->GetLabel(1) != trackIndex) && - (c->GetLabel(2) != trackIndex)) t.AddNWrong(); - } - + if((c->GetLabel(0) != trackIndex) && + (c->GetLabel(1) != trackIndex) && + (c->GetLabel(2) != trackIndex)) t.AddNWrong(); + } + } if(!cl) { - - for (Int_t i=timeBin.Find(y-road); iGetY() > y+road) break; if (c->IsUsed() > 0) continue; - if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue; - + // if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue; + if((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue; + // + // Double_t h01 = GetTiltFactor(c); chi2=t.GetPredictedChi2(c,h01); @@ -1185,55 +1191,29 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) if (cl) { wYclosest = cl->GetY(); wZclosest = cl->GetZ(); - - t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); + if (cl->GetNPads()<5) + t.SetSampledEdx(cl->GetQ()/dx); Double_t h01 = GetTiltFactor(cl); Int_t det = cl->GetDetector(); Int_t plane = fGeom->GetPlane(det); - + if (t.fX>345){ + t.fNLast++; + t.fChi2Last+=maxChi2; + } if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) { - //if(!t.Update(cl,maxChi2,index,h01)) { - if(!tryAgain--) return 0; + if(!t.Update(cl,maxChi2,index,h01)) { + if(!tryAgain--) return 0; + } } else tryAgain=fMaxGap; } else { if (tryAgain==0) break; - tryAgain--; - - //if (minDY < 1000000 && isNewLayer) - //cout << "\t" << nRefPlane << "\t" << "\t" << t.GetNRotate() << "\t" << - // road << "\t" << minDY << "\t" << chi2 << "\t" << wChi2 << "\t" << maxChi2 << endl; - + tryAgain--; } isNewLayer = kFALSE; - /* - if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) { - - printf(" %f", wIndex); //1 - printf(" %f", wTB); //2 - printf(" %f", wYrt); //3 - printf(" %f", wYclosest); //4 - printf(" %f", wYcorrect); //5 - printf(" %f", wYwindow); //6 - printf(" %f", wZrt); //7 - printf(" %f", wZclosest); //8 - printf(" %f", wZcorrect); //9 - printf(" %f", wZwindow); //10 - printf(" %f", wPx); //11 - printf(" %f", wPy); //12 - printf(" %f", wPz); //13 - printf(" %f", wSigmaC2*1000000); //14 - printf(" %f", wSigmaTgl2*1000); //15 - printf(" %f", wSigmaY2); //16 - // printf(" %f", wSigmaZ2); //17 - printf(" %f", wChi2); //17 - printf(" %f", wC); //18 - printf("\n"); - } - */ } } } @@ -1335,7 +1315,9 @@ Int_t AliTRDtracker::Refit(AliTRDtrack& t, Int_t rf) AliTRDcluster *cl=(AliTRDcluster*)GetCluster(iCluster[nr-1]); Double_t h01 = GetTiltFactor(cl); Double_t chi2=t.GetPredictedChi2(cl, h01); - t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); + if (cl->GetNPads()<5) t.SetSampledEdx(cl->GetQ()/dx); + + //t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); t.Update(cl,chi2,iCluster[nr-1],h01); } @@ -1464,7 +1446,7 @@ Int_t AliTRDtracker::LoadClusters(TTree *cTree) // Fills clusters into TRD tracking_sectors // Note that the numbering scheme for the TRD tracking_sectors // differs from that of TRD sectors - + cout<<"\n Read Sectors clusters"<GetTotBytes()/(sizeof(AliTRDcluster))); + TObjArray *clusterArray = new TObjArray(nsize+1000); TBranch *branch=ClusterTree->GetBranch("TRDcluster"); if (!branch) { @@ -1739,7 +1722,6 @@ Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) const Int_t nbytes = 0; AliTRDcluster *c = 0; // printf("\n"); - for (Int_t iEntry = 0; iEntry < nEntries; iEntry++) { // Import the tree @@ -1752,15 +1734,22 @@ Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) const // Loop through all TRD digits for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster); - AliTRDcluster *co = new AliTRDcluster(*c); + if (c->GetNPads()>3&&(iCluster%3>0)) { + delete clusterArray->RemoveAt(iCluster); + continue; + } + // AliTRDcluster *co = new AliTRDcluster(*c); //remove unnecesary coping - + clusters are together in memory + AliTRDcluster *co = c; co->SetSigmaY2(c->GetSigmaY2() * fSY2corr); Int_t ltb = co->GetLocalTimeBin(); if(ltb == 19) co->SetSigmaZ2(c->GetSigmaZ2()); else if(fNoTilt) co->SetSigmaZ2(c->GetSigmaZ2() * fSZ2corr); array->AddLast(co); - delete clusterArray->RemoveAt(iCluster); + // delete clusterArray->RemoveAt(iCluster); + clusterArray->RemoveAt(iCluster); } } + cout<<"Allocated"<GetEntriesFast()<<"\n"; delete clusterArray; @@ -2211,6 +2200,7 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I MapTimeBinLayers(); delete [] zc; delete [] zmax; + delete [] zmaxsensitive; } @@ -2428,32 +2418,46 @@ void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters( // and sensitivity in point // + Double_t alpha = AliTRDgeometry::GetAlpha(); + Double_t ymax = fX*TMath::Tan(0.5*alpha); + + dx = fdX; rho = fRho; radLength = fX0; lookForCluster = kFALSE; + Bool_t cross =kFALSE; + // + // + if ( (ymax-TMath::Abs(y))<3.){ //cross material + rho*=40.; + radLength*=40.; + cross=kTRUE; + } // // check dead regions in sensitive volume - if(fTimeBinIndex >= 0) { // - Int_t zone=-1; - for(Int_t ch = 0; ch < (Int_t) kZones; ch++) { - if (TMath::Abs(z - fZc[ch]) < fZmaxSensitive[ch]){ - zone = ch; - lookForCluster = !(fIsHole[zone]); - if(TMath::Abs(y) > fYmaxSensitive){ - lookForCluster = kFALSE; - } - if (fIsHole[zone]) { - //if hole - rho = 1.29e-3; - radLength = 36.66; - } - } - } - return; + Int_t zone=-1; + for(Int_t ch = 0; ch < (Int_t) kZones; ch++) { + if (TMath::Abs(z - fZc[ch]) > fZmax[ch]) continue; //not in given zone + // + if (TMath::Abs(z - fZc[ch]) < fZmaxSensitive[ch]){ + if (fTimeBinIndex>=0) lookForCluster = !(fIsHole[zone]); + if(TMath::Abs(y) > fYmaxSensitive){ + lookForCluster = kFALSE; + } + if (fIsHole[zone]) { + //if hole + rho = 1.29e-3; + radLength = 36.66; + } + }else{ + rho = 2.7; radLength = 24.01; //aluminium in between + } } // + if (fTimeBinIndex>=0) return; + // // // check hole if (fHole==kFALSE) return; @@ -2464,7 +2468,7 @@ void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters( //if hole rho = 1.29e-3; radLength = 36.66; - } + } } } return; @@ -2610,6 +2614,169 @@ void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack) } // end of function +Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters) +{ + // + // + // try to find nearest clusters to the track in timebins from t0 to t1 + // + // + Double_t x[1000],yt[1000],zt[10000]; + Double_t dz[1000],dy[10000]; + Int_t indexes[2][10000]; + AliTRDcluster *cl[2][10000]; + + for (Int_t it=t0;itGetX(); + Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt()); + Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl()); + Double_t road = 10.*sqrt(track->GetSigmaY2() + sy2); + Int_t nall=0; + Int_t nfound=0; + + for (Int_t it=t0;itGetLayer(it)); + if (timeBin==0) continue; // no indexes1 + Int_t maxn = timeBin; + x[it] = timeBin.GetX(); + Double_t y,z; + if (!track->GetProlongation(x[it],y,z)) continue; + yt[it]=y; + zt[it]=z; + Double_t chi2 =1000000; + nall++; + // + // find nearest cluster at given pad + for (Int_t i=timeBin.Find(y-road); iGetY() > y+road) break; + if (c->IsUsed() > 0) continue; + if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue; + chi2=track->GetPredictedChi2(c,h01); + if (chi2 > maxChi2) continue; + maxChi2=chi2; + cl[0][it]=c; + indexes[0][it] =timeBin.GetIndex(i); + } + // + // find nearest cluster also in adjacent 2 pads + // + for (Int_t i=timeBin.Find(y-road); iGetY() > y+road) break; + if (c->IsUsed() > 0) continue; + if((c->GetZ()-z)*(c->GetZ()-z) > 12 * sz2) continue; + chi2=track->GetPredictedChi2(c,h01); + if (chi2 > maxChi2) continue; + maxChi2=chi2; + cl[1][it]=c; + indexes[1][it]=timeBin.GetIndex(i); + if (!cl[0][it]) { + cl[0][it]=c; + indexes[0][it]=timeBin.GetIndex(i); + } + } + if (cl[0][it]||cl[1][it]) nfound++; + } + // + if (nfound<5) return 0; + // + // choose one of the variants + // + Int_t changes[2]={0,0}; + Float_t sigmay[2]={1000,1000}; + Float_t meany[2] ={1000,1000}; + Float_t meanz[2] ={1000,1000}; + Int_t sumall[2] ={0,0}; + Int_t ngood[2] ={0,0}; + Int_t nbad[2] ={0,0}; + // + // + for (Int_t ih=0; ih<2;ih++){ + Float_t lastz =-10000; + Float_t sumz =0; + Float_t sum =0; + Double_t sumdy = 0; + Double_t sumdy2= 0; + // + // how many changes ++ mean z ++mean y ++ sigma y + for (Int_t it=t0;itGetZ(); + if (TMath::Abs(lastz-cl[ih][it]->GetZ())>1) { + lastz = cl[ih][it]->GetZ(); + changes[ih]++; + } + sumz = cl[ih][it]->GetZ(); + sum++; + Double_t h01 = GetTiltFactor(cl[ih][it]); + dz[it] = cl[ih][it]->GetZ()- zt[it]; + dy[it] = cl[ih][it]->GetY()+ dz[it]*h01 -yt[it]; + sum++; + sumdy += dy[it]; + sumdy2+= dy[it]*dy[it]; + Int_t label = TMath::Abs(track->GetLabel()); + if (cl[ih][it]->GetLabel(0)==label || cl[ih][it]->GetLabel(1)==label||cl[ih][it]->GetLabel(2)==label){ + ngood[ih]++; + } + else{ + nbad[ih]++; + } + } + if (sumall[ih]>4){ + meanz[ih] = sumz/sum; + meany[ih] = sumdy/sum; + sigmay[ih] = TMath::Sqrt(sumdy2/sum-meany[ih]*meany[ih]); + } + } + Int_t best =0; + /* + if (sumall[0]quality1){ + best = 1; + } + + + // + for (Int_t it=t0;itGetZ()- zt[it]; + dy[it] = cl[best][it]->GetY()+ dz[it]*h01 -yt[it]; + // + if (TMath::Abs(dy[it])<2.5*sigmay[best]) + clusters[it] = indexes[best][it]; + } + + if (nbad[0]>4){ + nbad[0] = nfound; + } + return nfound; +} + + diff --git a/TRD/AliTRDtracker.h b/TRD/AliTRDtracker.h index 73ec38c094d..fe71d27c492 100644 --- a/TRD/AliTRDtracker.h +++ b/TRD/AliTRDtracker.h @@ -52,6 +52,7 @@ class AliTRDtracker : public AliTracker { Int_t CookSectorIndex(Int_t gs) const { return kTrackingSectors - 1 - gs; } AliTRDcluster * GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin); Int_t GetLastPlane(AliTRDtrack * track); //return last updated plane + Int_t FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters); Float_t GetSeedGap() const {return fgkSeedGap;} Int_t GetMaxGap() const {return fMaxGap;} -- 2.43.0