From 7ad193388829ce599187d8ac16cbd0fb3891205f Mon Sep 17 00:00:00 2001 From: hristov Date: Tue, 21 Jun 2005 15:26:01 +0000 Subject: [PATCH] Updated version of the TRD reconstruction (M.Ivanov) --- TRD/AliTRDclusterizerV1.cxx | 137 ++-- TRD/AliTRDclusterizerV1.h | 3 +- TRD/AliTRDtrack.cxx | 270 +++++--- TRD/AliTRDtrack.h | 143 +++-- TRD/AliTRDtracker.cxx | 1168 +++++++++++++++++++++++------------ TRD/AliTRDtracker.h | 29 +- TRD/TRDrecLinkDef.h | 1 + 7 files changed, 1116 insertions(+), 635 deletions(-) diff --git a/TRD/AliTRDclusterizerV1.cxx b/TRD/AliTRDclusterizerV1.cxx index 24991b29fc9..c41b2551101 100644 --- a/TRD/AliTRDclusterizerV1.cxx +++ b/TRD/AliTRDclusterizerV1.cxx @@ -164,10 +164,10 @@ Bool_t AliTRDclusterizerV1::MakeClusters() } fPar->Init(); - Float_t timeBinSize = fPar->GetDriftVelocity() - / fPar->GetSamplingFrequency(); + //Float_t timeBinSize = fPar->GetDriftVelocity() + // / fPar->GetSamplingFrequency(); // Half of ampl.region - const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; + // const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; Float_t omegaTau = fPar->GetOmegaTau(); if (fVerbose > 0) { @@ -186,7 +186,6 @@ Bool_t AliTRDclusterizerV1::MakeClusters() Int_t maxThresh = fPar->GetClusMaxThresh(); // Threshold value for the digit signal Int_t sigThresh = fPar->GetClusSigThresh(); - // Iteration limit for unfolding procedure const Float_t kEpsilon = 0.01; @@ -195,11 +194,11 @@ Bool_t AliTRDclusterizerV1::MakeClusters() const Int_t kNtrack = 3 * kNclus; Int_t iType = 0; - Int_t iUnfold = 0; - + Int_t iUnfold = 0; Double_t ratioLeft = 1.0; Double_t ratioRight = 1.0; + // Double_t padSignal[kNsig]; Double_t clusterSignal[kNclus]; Double_t clusterPads[kNclus]; @@ -260,12 +259,19 @@ Bool_t AliTRDclusterizerV1::MakeClusters() Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time)); Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time)); +// // Look for the maximum +// if (signalM >= maxThresh) { +// if (((signalL >= sigThresh) && +// (signalL < signalM)) || +// ((signalR >= sigThresh) && +// (signalR < signalM))) { +// // Maximum found, mark the position by a negative signal +// digits->SetDataUnchecked(row,col-1,time,-signalM); +// } +// } // Look for the maximum if (signalM >= maxThresh) { - if (((signalL >= sigThresh) && - (signalL < signalM)) || - ((signalR >= sigThresh) && - (signalR < signalM))) { + if ( (signalL<=signalM) && (signalR<=signalM) && (signalL+signalR)>sigThresh ) { // Maximum found, mark the position by a negative signal digits->SetDataUnchecked(row,col-1,time,-signalM); } @@ -336,10 +342,7 @@ Bool_t AliTRDclusterizerV1::MakeClusters() break; }; - // Don't analyze large clusters - //if (iType == 4) continue; - - // Look for 5 pad cluster with minimum in the middle + // Look for 5 pad cluster with minimum in the middle Bool_t fivePadCluster = kFALSE; if (col < nColMax-3) { if (digits->GetDataUnchecked(row,col+2,time) < 0) { @@ -362,7 +365,7 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // of the cluster which remains from a previous unfolding if (iUnfold) { clusterSignal[0] *= ratioLeft; - iType = 3; + iType = 5; iUnfold = 0; } @@ -378,10 +381,11 @@ Bool_t AliTRDclusterizerV1::MakeClusters() ratioRight = Unfold(kEpsilon,iplan,padSignal); ratioLeft = 1.0 - ratioRight; clusterSignal[2] *= ratioRight; - iType = 3; + iType = 5; iUnfold = 1; } + Double_t clusterCharge = clusterSignal[0] + clusterSignal[1] + clusterSignal[2]; @@ -392,29 +396,31 @@ Bool_t AliTRDclusterizerV1::MakeClusters() clusterPads[2] = time - nTimeBefore + 0.5; if (fPar->LUTOn()) { - // Calculate the position of the cluster by using the // lookup table method -// clusterPads[1] = col + 0.5 -// + fPar->LUTposition(iplan,clusterSignal[0] -// ,clusterSignal[1] -// ,clusterSignal[2]); - clusterPads[1] = 0.5 - + fPar->LUTposition(iplan,clusterSignal[0] + clusterPads[1] = + fPar->LUTposition(iplan,clusterSignal[0] ,clusterSignal[1] ,clusterSignal[2]); - } else { - // Calculate the position of the cluster by using the // center of gravity method -// clusterPads[1] = col + 0.5 -// + (clusterSignal[2] - clusterSignal[0]) -// / clusterCharge; - clusterPads[1] = 0.5 - + (clusterSignal[2] - clusterSignal[0]) - / clusterCharge; + for (Int_t i=0;i<5;i++) padSignal[i]=0; + padSignal[2] = TMath::Abs(digits->GetDataUnchecked(row,col,time)); // central pad + padSignal[1] = TMath::Abs(digits->GetDataUnchecked(row,col-1,time)); // left pad + padSignal[3] = TMath::Abs(digits->GetDataUnchecked(row,col+1,time)); // right pad + if (col>2 &&TMath::Abs(digits->GetDataUnchecked(row,col-2,time)GetDataUnchecked(row,col-2,time)); + } + if (colGetDataUnchecked(row,col+2,time)GetDataUnchecked(row,col+2,time)); + } + clusterPads[1] = GetCOG(padSignal); + Double_t check = fPar->LUTposition(iplan,clusterSignal[0] + ,clusterSignal[1] + ,clusterSignal[2]); + // clusterPads[1] = check; } @@ -424,49 +430,16 @@ Bool_t AliTRDclusterizerV1::MakeClusters() Double_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) / (clusterCharge*clusterCharge); - if (fVerbose > 1) { - printf("-----------------------------------------------------------\n"); - printf("Create cluster no. %d\n",nClusters); - printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0] - ,clusterPads[1] - ,clusterPads[2]); - printf("Indices: %d, %d, %d\n",clusterDigit[0] - ,clusterDigit[1] - ,clusterDigit[2]); - printf("Total charge = %f\n",clusterCharge); - printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0] - ,clusterTracks[1] - ,clusterTracks[2]); - printf(" pad1 %d, %d, %d\n",clusterTracks[3] - ,clusterTracks[4] - ,clusterTracks[5]); - printf(" pad2 %d, %d, %d\n",clusterTracks[6] - ,clusterTracks[7] - ,clusterTracks[8]); - printf("Type = %d, Number of pads = %d\n",iType,nPadCount); - } - // Calculate the position and the error + Double_t colSize = padPlane->GetColSize(col); + Double_t rowSize = padPlane->GetRowSize(row); Double_t clusterPos[3]; -// clusterPos[0] = clusterPads[1] * colSize + col0; -// clusterPos[1] = clusterPads[0] * rowSize + row0; - clusterPos[0] = padPlane->GetColPos(col) - clusterPads[1]; - clusterPos[1] = padPlane->GetRowPos(row) - clusterPads[0]; + clusterPos[0] = padPlane->GetColPos(col) + (clusterPads[1]-0.5)*colSize; // MI change + clusterPos[1] = padPlane->GetRowPos(row) -0.5*rowSize; //MI change clusterPos[2] = clusterPads[2]; Double_t clusterSig[2]; - Double_t colSize = padPlane->GetColSize(col); - Double_t rowSize = padPlane->GetRowSize(row); clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize; - clusterSig[1] = rowSize * rowSize / 12.; - - // Correct for ExB displacement - if (fPar->ExBOn()) { - Int_t local_time_bin = (Int_t) clusterPads[2]; - Double_t driftLength = local_time_bin * timeBinSize + kAmWidth; - Double_t deltaY = omegaTau * driftLength; - clusterPos[1] = clusterPos[1] - deltaY; - } - + clusterSig[1] = rowSize * rowSize / 12.; // Add the cluster to the output array AddCluster(clusterPos ,idet @@ -483,24 +456,12 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // Compress the arrays digits->Compress(1,0); track0->Compress(1,0); - track1->Compress(1,0); + track1->Compress(1,0); track2->Compress(1,0); // Write the cluster and reset the array WriteClusters(idet); ResetRecPoints(); - - if (fVerbose > 0) { - printf(" "); - printf("Found %d clusters in total.\n" - ,nClusters); - printf(" 2pad: %d\n",nClusters2pad); - printf(" 3pad: %d\n",nClusters3pad); - printf(" 4pad: %d\n",nClusters4pad); - printf(" 5pad: %d\n",nClusters5pad); - printf(" Large: %d\n",nClustersLarge); - } - } } } @@ -514,6 +475,18 @@ Bool_t AliTRDclusterizerV1::MakeClusters() } +Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5]) +{ + // + // get COG position + // used for clusters with more than 3 pads - where LUT not applicable + Double_t sum = signal[0]+signal[1]+signal[2]+signal[3]+signal[4]; + Double_t res = (0.0*(-signal[0]+signal[4])+(-signal[1]+signal[3]))/sum; + return res; +} + + + //_____________________________________________________________________________ Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSignal) { diff --git a/TRD/AliTRDclusterizerV1.h b/TRD/AliTRDclusterizerV1.h index 5d5e461784a..52c7a14c197 100644 --- a/TRD/AliTRDclusterizerV1.h +++ b/TRD/AliTRDclusterizerV1.h @@ -33,9 +33,8 @@ class AliTRDclusterizerV1 : public AliTRDclusterizer { AliTRDdigitsManager *fDigitsManager; //! TRD digits manager private: - virtual Double_t Unfold(Double_t eps, Int_t plane, Double_t *padSignal); - + Double_t GetCOG(Double_t signal[5]); // get COG position ClassDef(AliTRDclusterizerV1,5) // TRD-Cluster finder, slow simulator }; diff --git a/TRD/AliTRDtrack.cxx b/TRD/AliTRDtrack.cxx index 214f814a932..aaba24d6bd0 100644 --- a/TRD/AliTRDtrack.cxx +++ b/TRD/AliTRDtrack.cxx @@ -23,14 +23,18 @@ #include "AliTRDtrack.h" #include "AliTRDclusterCorrection.h" +ClassImp(AliTRDtracklet) ClassImp(AliTRDtrack) + + AliTRDtracklet::AliTRDtracklet():fY(0),fX(0),fAlpha(0),fSigma2(0),fP0(0),fP1(0),fNFound(0),fNCross(0),fPlane(0),fExpectedSigma2(0),fChi2(0){ +} + //_____________________________________________________________________________ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, const Double_t xx[5], const Double_t cc[15], - Double_t xref, Double_t alpha) : AliKalmanTrack() -{ + Double_t xref, Double_t alpha) : AliKalmanTrack() { //----------------------------------------------------------------- // This is the main track constructor. //----------------------------------------------------------------- @@ -77,7 +81,7 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, fdQdl[0] = q; // initialisation [SR, GSI 18.02.2003] (i startd for 1) - for(UInt_t i=1; i2) return -1; // sure it's stopped 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; - // + + //comp->fTree->SetAlias("nlast2","track.fTracklets[5].fNFound+track.fTracklets[4].fNFound"); + //comp->fTree->SetAlias("goldtrack","abs((track.fTracklets[5].fP1+track.fTracklets[4].fP1))<0.5&&nlast2>14"); + Int_t nlast2 = fTracklets[5].fNFound+fTracklets[4].fNFound; + if (abs((fTracklets[5].fP1+fTracklets[4].fP1))<0.3 &&nlast2>20) return 3; + if (abs((fTracklets[5].fP1+fTracklets[4].fP1))<0.3 &&nlast2>14) return 2; + if (abs((fTracklets[5].fP1+fTracklets[4].fP1))<0.5 &&nlast2>14) return 1; return status; } + //____________________________________________________________________________ void AliTRDtrack::GetExternalParameters(Double_t& xr, Double_t x[5]) const { // @@ -356,12 +359,9 @@ void AliTRDtrack::GetExternalCovariance(Double_t cc[15]) const { } + //_____________________________________________________________________________ -void AliTRDtrack::GetCovariance(Double_t cc[15]) const -{ - // - // Returns the covariance matrix - // +void AliTRDtrack::GetCovariance(Double_t cc[15]) const { cc[0]=fCyy; cc[1]=fCzy; cc[2]=fCzz; @@ -402,33 +402,22 @@ void AliTRDtrack::CookdEdx(Double_t low, Double_t up) { return; } - Float_t sorted[kMAXCLUSTERSPERTRACK]; + Float_t sorted[kMAX_CLUSTERS_PER_TRACK]; for (i=0; i < nc; i++) { - sorted[i]=fdQdl[i]; + sorted[i]=TMath::Abs(fdQdl[i]); } - /* - Int_t swap; - - do { - swap=0; - for (i=0; i2.) beta2*= 0.4; // magic const - theta2 for relativistic particles + // - not valid for electrons Double_t ey=fC*fX - fE, ez=fT; Double_t xz=fC*ez, zz1=ez*ez+1, xy=fE+ey; @@ -522,12 +513,18 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho) //Energy losses************************ if((5940*beta2/(1-beta2+1e-10) - beta2) < 0) return 0; - Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho; + Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho; + dE *= 1.2; // magic const - to normalize pools - waiting for geo manager if (x1 < x2) dE=-dE; cc=fC; fC*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE); fE+=fX*(fC-cc); - + // + Double_t deltac = fC-cc; + fCcc += 4*deltac*deltac; // fluctuation of energy losses + fCee += 4*fX*fX*deltac*deltac; // local angle unchanged + fCce += 4*fX*deltac*deltac; // correlation 1 + // // track time measurement [SR, GSI 17.02.2002] if (x1 < x2) if (IsStartedTimeIntegral()) { @@ -538,6 +535,7 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho) return 1; } + //_____________________________________________________________________________ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index, Double_t h01) { @@ -581,7 +579,7 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index, fC = cur; } else { - Double_t xuFactor = 100.; // empirical factor set by C.Xu + Double_t xu_factor = 100.; // empirical factor set by C.Xu // in the first tilt version dy=c->GetY() - fY; dz=c->GetZ() - fZ; dy=dy+h01*dz; @@ -596,7 +594,7 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index, - r00=c->GetSigmaY2()+errang+add, r01=0., r11=c->GetSigmaZ2()*xuFactor; + r00=c->GetSigmaY2()+errang+add, r01=0., r11=c->GetSigmaZ2()*xu_factor; r00+=(fCyy+2.0*h01*fCzy+h01*h01*fCzz); r01+=(fCzy+h01*fCzz); @@ -659,7 +657,6 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index, return 1; } - //_____________________________________________________________________________ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index, Double_t h01, Int_t /*plane*/) @@ -677,14 +674,7 @@ 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 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; + Double_t errsys =(0.025*0.025*20)*(1+tangent2); //systematic error part /* if (corrector!=0){ @@ -699,9 +689,7 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index } */ // - // Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.); - - Double_t r00=(c->GetSigmaY2() +errang+errsys)*extend, r01=0., r11=c->GetSigmaZ2()*10000.; + Double_t r00=(c->GetSigmaY2()+errsys), 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; @@ -729,30 +717,17 @@ 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 padlength = TMath::Sqrt(c->GetSigmaZ2()*12); - Double_t xuFactor = 1000.; // empirical factor set by C.Xu + 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; 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; - } - Double_t s00 = (c->GetSigmaY2()+errang)*extend+errsys+add; // error pad - Double_t s11 = c->GetSigmaZ2()*xuFactor; // error pad-row + Double_t s00 = c->GetSigmaY2()+errsys; // error pad + Double_t s11 = c->GetSigmaZ2()*xu_factor; // error pad-row // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00; r01 = fCzy + fCzz*h01; @@ -828,6 +803,116 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index return 1; } + + + + +//_____________________________________________________________________________ +Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet) +{ + // + // Assignes found tracklet to the track and updates track information + // + // + Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 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; + // + + Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ; + + + Double_t s00 = tracklet.GetTrackletSigma2(); // error pad + Double_t s11 = 100000; // error pad-row + Float_t h01 = tracklet.GetTilt(); + // + // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00; + r00 = fCyy + fCzz*h01*h01+s00; + // r01 = fCzy + fCzz*h01; + r01 = fCzy ; + r11 = fCzz + s11; + det = r00*r11 - r01*r01; + // inverse matrix + tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det; + + Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11; + Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11; + Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11; + Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11; + Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11; + + // K matrix +// k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01); +// k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01); +// k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01); +// k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01); +// k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01); + // + //Update measurement + Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz; + // cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz; + if (TMath::Abs(cur*fX-eta) >= 0.90000) { + //Int_t n=GetNumberOfClusters(); + // if (n>4) cerr<SetPHOShole(); fGeom->SetRICHhole(); } @@ -135,12 +155,14 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker() if (!fPar) { printf("AliTRDtracker::AliTRDtracker(): can't find TRD parameter!\n"); printf("The DEFAULT TRD parameter will be used\n"); - fPar = new AliTRDparameter(); + fPar = new AliTRDparameter("Pica","Vyjebana"); } + fPar = new AliTRDparameter("Pica","Vyjebana"); fPar->Init(); savedir->cd(); + // fGeom->SetT0(fTzero); fNclusters = 0; @@ -157,9 +179,9 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker() fHoles[icham][trS]=fGeom->IsHole(0,icham,geomS); } } - AliTRDpadPlane *padPlane = fPar->GetPadPlane(0,0); - Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle()); + Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle()); + // Float_t tiltAngle = TMath::Abs(fPar->GetTiltingAngle()); if(tiltAngle < 0.1) { fNoTilt = kTRUE; } @@ -175,8 +197,9 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker() Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region - Double_t dx = (Double_t) fPar->GetDriftVelocity() + Double_t dx = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity() / fPar->GetSamplingFrequency(); + Int_t tbAmp = fPar->GetTimeBefore(); Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx); if(kTRUE) maxAmp = 0; // intentional until we change the parameter class @@ -190,6 +213,8 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker() fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fgkSkipDepth); fVocal = kFALSE; + + fDebugStreamer = new TTreeSRedirector("TRDdebug.root"); savedir->cd(); } @@ -219,7 +244,10 @@ AliTRDtracker::~AliTRDtracker() for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) { delete fTrSec[geomS]; } - + if (fDebugStreamer) { + //fDebugStreamer->Close(); + delete fDebugStreamer; + } } //_____________________________________________________________________ @@ -300,7 +328,7 @@ inline Double_t f3trd(Double_t x1,Double_t y1, } -AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin){ +AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin, UInt_t &index){ // //try to find cluster in the backup list // @@ -314,6 +342,7 @@ AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_ Int_t iplane = fGeom->GetPlane(cli->GetDetector()); if (iplane==plane) { cl = cli; + index = indexes[i]; break; } } @@ -360,7 +389,7 @@ Int_t AliTRDtracker::Clusters2Tracks(AliESD* event) if ( (status & AliESDtrack::kTRDout ) == 0 ) continue; if ( (status & AliESDtrack::kTRDin) != 0 ) continue; nseed++; - + AliTRDtrack* seed2 = new AliTRDtrack(*seed); //seed2->ResetCovariance(); AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha()); @@ -512,11 +541,13 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { Bool_t isGold = kFALSE; if (track->GetChi2()/track->GetNumberOfClusters()<5) { //full gold track - seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); + // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); + if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup); isGold = kTRUE; } if (!isGold && track->GetNCross()==0&&track->GetChi2()/track->GetNumberOfClusters()<7){ //almost gold track - seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); + // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); + if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup); isGold = kTRUE; } if (!isGold && track->GetBackupTrack()){ @@ -526,11 +557,14 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { isGold = kTRUE; } } + if (track->StatusForTOF()>0 &&track->fNCross==0 && Float_t(track->fN)/Float_t(track->fNExpected)>0.4){ + seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup); + } } } // //Propagation to the TOF (I.Belikov) - CookdEdxTimBin(*track); + if (track->GetStop()==kFALSE){ Double_t xtof=371.; @@ -570,7 +604,7 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { seed->SetTRDsignals(track->GetPIDsignals(i),i); seed->SetTRDTimBin(track->GetPIDTimBin(i),i); } - seed->SetTRDtrack(new AliTRDtrack(*track)); + // seed->SetTRDtrack(new AliTRDtrack(*track)); if (track->GetNumberOfClusters()>foundMin) found++; } }else{ @@ -581,13 +615,32 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { seed->SetTRDsignals(track->GetPIDsignals(i),i); seed->SetTRDTimBin(track->GetPIDTimBin(i),i); } - seed->SetTRDtrack(new AliTRDtrack(*track)); + //seed->SetTRDtrack(new AliTRDtrack(*track)); found++; } } - seed->SetTRDQuality(track->StatusForTOF()); + seed->SetTRDQuality(track->StatusForTOF()); + // + // Debug part of tracking + TTreeSRedirector& cstream = *fDebugStreamer; + Int_t eventNr = event->GetEventNumber(); + if (track->GetBackupTrack()){ + cstream<<"Tracks"<< + "EventNr="<foundMin) // seed->UpdateTrackParams(track, AliESDtrack::kTRDout); @@ -598,6 +651,9 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) { cerr<<"Number of seeds: "<Clear(); fNseeds=0; delete [] index; delete [] quality; @@ -641,22 +697,22 @@ Int_t AliTRDtracker::RefitInward(AliESD* event) continue; } nseed++; - 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.); - } +// if (1/seed2.Get1Pt()>1.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); -// } + for (Int_t i=0;iSetPIDsignals(seed2.GetPIDsignals(i),i); + pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i); + } UInt_t * indexes3 = pt->GetBackupIndexes(); for (Int_t i=0;i<200;i++) { @@ -669,17 +725,18 @@ Int_t AliTRDtracker::RefitInward(AliESD* event) if (t.GetNumberOfClusters() >= foundMin) { // UseClusters(&t); //CookLabel(pt, 1-fgkLabelFraction); - // t.CookdEdx(); + t.CookdEdx(); + CookdEdxTimBin(t); } found++; // cout<UpdateTrackParams(pt, AliESDtrack::kTRDrefit); - // for (Int_t i=0;iSetTRDsignals(pt->GetPIDsignals(i),i); -// seed->SetTRDTimBin(pt->GetPIDTimBin(i),i); -// } + for (Int_t i=0;iSetTRDsignals(pt->GetPIDsignals(i),i); + seed->SetTRDTimBin(pt->GetPIDTimBin(i),i); + } }else{ //if not prolongation to TPC - propagate without update AliTRDtrack* seed2 = new AliTRDtrack(*seed); @@ -690,10 +747,10 @@ Int_t AliTRDtracker::RefitInward(AliESD* event) pt2->CookdEdx(0.,1.); CookdEdxTimBin(*pt2); seed->UpdateTrackParams(pt2, AliESDtrack::kTRDrefit); - // for (Int_t i=0;iSetTRDsignals(pt2->GetPIDsignals(i),i); -// seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i); -// } + for (Int_t i=0;iSetTRDsignals(pt2->GetPIDsignals(i),i); + seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i); + } } delete pt2; } @@ -854,7 +911,7 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf) Int_t plane = fGeom->GetPlane(cl0->GetDetector()); if (plane>lastplane) continue; Int_t timebin = cl0->GetLocalTimeBin(); - AliTRDcluster * cl2= GetCluster(&t,plane, timebin); + AliTRDcluster * cl2= GetCluster(&t,plane, timebin,index); if (cl2) { cl =cl2; Double_t h01 = GetTiltFactor(cl); @@ -867,44 +924,6 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf) continue; } - /* - 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); - Double_t 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) > 12 * sz2) continue; - - Double_t h01 = GetTiltFactor(c); - Double_t chi2=t.GetPredictedChi2(c, h01); - - if (chi2 > maxChi2) continue; - maxChi2=chi2; - cl=c; - index=timeBin.GetIndex(i); - } - } - */ if (cl) { wYclosest = cl->GetY(); @@ -912,9 +931,7 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf) Double_t h01 = GetTiltFactor(cl); 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()); + t.SetSampledEdx(TMath::Abs(cl->GetQ()/dx)); Int_t det = cl->GetDetector(); Int_t plane = fGeom->GetPlane(det); @@ -945,13 +962,6 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) // layers confirms prolongation if a close cluster is found. // Returns the number of clusters expected to be found in sensitive layers - - Float_t wIndex, wTB, wChi2; - Float_t wYrt, wYclosest, wYcorrect, wYwindow; - Float_t wZrt, wZclosest, wZcorrect, wZwindow; - Float_t wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2; - - Int_t trackIndex = t.GetLabel(); Int_t tryAgain=fMaxGap; Double_t alpha=t.GetAlpha(); @@ -963,7 +973,8 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) 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; + //Double_t radLength, rho, x, dx, y, ymax = 0, z; + Double_t radLength, rho, x, dx, y, z; Bool_t lookForCluster; Int_t expectedNumberOfClusters = 0; @@ -971,18 +982,15 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning - Int_t nRefPlane = fgkFirstPlane; - Bool_t isNewLayer = kFALSE; - - Double_t chi2; - Double_t minDY; - Int_t zone =-10; + // Int_t zone =0; Int_t nr; + Float_t ratio0=0; + AliTRDtracklet tracklet; + // for (nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr0.000001) cross = kTRUE; //better to stop track - Int_t currentzone = fTrSec[s]->GetLayer(nr)->GetZone(z); - 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; - 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()>4) break; - } + //Int_t nrotate = t.GetNRotate(); + if (!AdjustSector(&t)) break; // // - s = t.GetSector(); - if (!t.PropagateTo(x,radLength,rho)) break; - y = t.GetY(); z = t.GetZ(); - - // Barrel Tracks [SR, 04.04.2003] - s = t.GetSector(); - if (fTrSec[s]->GetLayer(nr)->IsSensitive() != - fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) { - -// if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack); - } - - if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() && - ! fTrSec[s]->GetLayer(nr)->IsSensitive()) { - isNewLayer = kTRUE; - } else {isNewLayer = kFALSE;} - - y = t.GetY(); - z = t.GetZ(); // 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 - } +// if (nrotate!=t.GetNRotate()){ +// 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 (!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(); - // if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax); - // printf("label %d, pl %d, lookForCluster %d \n", - // trackIndex, nr+1, lookForCluster); - if(lookForCluster) { -// if (clusters[nr]==-1) { -// FindClusters(s,nr,nr+30,&t,clusters); -// } + if (clusters[nr]==-1) { + Float_t ncl = FindClusters(s,nr,nr+30,&t,clusters,tracklet); + ratio0 = ncl/Float_t(fTimeBinsPerPlane); + Float_t ratio1 = Float_t(t.fN+1)/Float_t(t.fNExpected+1.); + if (tracklet.GetChi2()<18.&&ratio0>0.8 && ratio1>0.6 && ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85&&t.fN>20){ + t.MakeBackupTrack(); // make backup of the track until is gold + } +// if (ncl>4){ +// t.PropagateTo(tracklet.GetX()); +// t.UpdateMI(tracklet); +// nr = fTrSec[0]->GetLayerNumber(t.GetX())+1; +// continue; +// } + } expectedNumberOfClusters++; t.fNExpected++; if (t.fX>345) t.fNExpectedLast++; - wIndex = (Float_t) t.GetLabel(); - wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex(); AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr+1)); Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt()); - Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl()); - if((t.GetSigmaY2() + sy2) < 0) break; + if((t.GetSigmaY2() + sy2) < 0) { + printf("problem\n"); + break; + } Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2); - Double_t y=t.GetY(), z=t.GetZ(); - - wYrt = (Float_t) y; - wZrt = (Float_t) z; - wYwindow = (Float_t) road; - wSigmaC2 = (Float_t) t.GetSigmaC2(); - wSigmaTgl2 = (Float_t) t.GetSigmaTgl2(); - wSigmaY2 = (Float_t) t.GetSigmaY2(); - wSigmaZ2 = (Float_t) t.GetSigmaZ2(); - wChi2 = -1; if (road>fgkWideRoad) { - if (t.GetNumberOfClusters()>4) - cerr< 0) maxChi2 = 3 * maxChi2; - - wYclosest = 12345678; - wYcorrect = 12345678; - wZclosest = 12345678; - wZcorrect = 12345678; - wZwindow = TMath::Sqrt(2.25 * 12 * sz2); - - // Find the closest correct cluster for debugging purposes - if (timeBin&&fVocal) { - minDY = 1000000; - for (Int_t i=0; iGetLabel(0) != trackIndex) && - (c->GetLabel(1) != trackIndex) && - (c->GetLabel(2) != trackIndex)) continue; - if(TMath::Abs(c->GetY() - y) > minDY) continue; - //minDY = TMath::Abs(c->GetY() - y); - minDY = c->GetY() - y; - wYcorrect = c->GetY(); - wZcorrect = c->GetZ(); - - Double_t h01 = GetTiltFactor(c); - wChi2 = t.GetPredictedChi2(c, h01); - } - } - // Now go for the real cluster search if (timeBin) { - /* - if (clusters[nr+1]>0) { + + if (clusters[nr+1]>0) { index = clusters[nr+1]; cl = (AliTRDcluster*)GetCluster(index); Double_t h01 = GetTiltFactor(cl); maxChi2=t.GetPredictedChi2(cl,h01); - } - */ - - 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(!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) > 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); - - if (chi2 > maxChi2) continue; - maxChi2=chi2; - cl=c; - index=timeBin.GetIndex(i); - } - } - + if (cl) { - wYclosest = cl->GetY(); - wZclosest = cl->GetZ(); if (cl->GetNPads()<5) - t.SetSampledEdx(cl->GetQ()/dx); + t.SetSampledEdx(TMath::Abs(cl->GetQ()/dx)); Double_t h01 = GetTiltFactor(cl); Int_t det = cl->GetDetector(); Int_t plane = fGeom->GetPlane(det); @@ -1205,18 +1086,26 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) } if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) { if(!t.Update(cl,maxChi2,index,h01)) { - if(!tryAgain--) return 0; + //if(!tryAgain--) return 0; } } else tryAgain=fMaxGap; - } + // + + if (cl->GetLocalTimeBin()==1&&t.fN>20 && float(t.fChi2)/float(t.fN)<5){ + Float_t ratio1 = Float_t(t.fN)/Float_t(t.fNExpected); + if (tracklet.GetChi2()<18&&ratio0>0.8&&ratio1>0.6 &&ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85){ + t.MakeBackupTrack(); // make backup of the track until is gold + } + } + + } else { - if (tryAgain==0) break; - tryAgain--; + // if (tryAgain==0) break; + //tryAgain--; } - - isNewLayer = kFALSE; - + + } } } @@ -1225,8 +1114,8 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t) else t.SetStop(kFALSE); return expectedNumberOfClusters; - - + + } //--------------------------------------------------------------------------- @@ -1318,7 +1207,7 @@ 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); - if (cl->GetNPads()<5) t.SetSampledEdx(cl->GetQ()/dx); + if (cl->GetNPads()<5) t.SetSampledEdx(TMath::Abs(cl->GetQ()/dx)); //t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); t.Update(cl,chi2,iCluster[nr-1],h01); @@ -1699,6 +1588,171 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn) } } } +//__________________________________________________________________________ +void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/) +{ + // + // Creates seeds using clusters between position inner plane and outer plane + // + + const Double_t maxtheta = 2; + const Double_t maxphi = 1.5; + Int_t maxSec=AliTRDgeometry::kNsect; + + // + // find the maximal and minimal layer for the planes + // fucking "object oriented" geometry - find the time bin range for different planes + // + Int_t layers[6][2]; + for (Int_t i=0;i<6;i++){layers[i][0]=10000; layers[i][1]=0;} + + for (Int_t ns=0;nsGetNumberOfLayers();ilayer++){ + AliTRDpropagationLayer& layer=*(fTrSec[ns]->GetLayer(ilayer)); + if (layer==0) continue; + Int_t det = layer[0]->GetDetector(); + Int_t plane = fGeom->GetPlane(det); + if (ilayerlayers[plane][1]) layers[plane][1] = ilayer; + } + } + // + + Int_t ilayer1 = layers[5][1]; // time bin in mplification region + Int_t ilayer2 = layers[3][1]; // + Int_t ilayerM = layers[4][1]; // + // + Double_t x1 = fTrSec[0]->GetX(ilayer1); + Double_t x2 = fTrSec[0]->GetX(ilayer2); + Double_t xm = fTrSec[0]->GetX(ilayerM); + Double_t dist = x2-x1; + // Int_t indexes1[20]; + //Int_t indexes2[20]; + AliTRDcluster *clusters1[15],*clusters2[15],*clustersM[15]; + // + // + for (Int_t ns=0; nsGetLayer(ilayer1)); //select propagation layers + AliTRDpropagationLayer& layer2=*(fTrSec[ns]->GetLayer(ilayer2)); + // + for (Int_t icl1=0;icl1GetY(); + Double_t z1 = cl1->GetZ(); + // + for (Int_t icl2=0;icl2GetY(); + Double_t z2 = cl2->GetZ(); + Double_t tanphi = (y2-y1)/dist; + Double_t tantheta = (z2-z1)/dist; + if (TMath::Abs(tanphi)>maxphi) continue; + if (TMath::Abs(tantheta)>maxtheta) continue; + // + clusters1[0] = cl1; + clusters2[0] = cl2; + Double_t road = 0.5+TMath::Abs(tanphi)*1; + Int_t ncl=0; + Double_t sum1=0, sumx1=0,sum2x1=0,sumxy1=0, sumy1=0; + Double_t sum2=0, sumx2=0,sum2x2=0,sumxy2=0, sumy2=0; + // + for (Int_t dlayer=1;dlayer<15;dlayer++){ + clusters1[dlayer]=0; + clusters2[dlayer]=0; + AliTRDpropagationLayer& layer1C=*(fTrSec[ns]->GetLayer(ilayer1-dlayer)); //select propagation layers + AliTRDpropagationLayer& layer2C=*(fTrSec[ns]->GetLayer(ilayer2-dlayer)); // + Double_t yy1 = y1+(tanphi) *(layer1C.GetX()-x1); + Double_t zz1 = z1+(tantheta)*(layer1C.GetX()-x1); + Double_t yy2 = y1+(tanphi) *(layer2C.GetX()-x1); + Double_t zz2 = z1+(tantheta)*(layer2C.GetX()-x1); + Int_t index1 = layer1C.FindNearestCluster(yy1,zz1,road); + Int_t index2 = layer2C.FindNearestCluster(yy2,zz2,road); + if (index1>=0) { + clusters1[dlayer]= (AliTRDcluster*)GetCluster(index1); + ncl++; + sum1++; + Double_t dx = layer1C.GetX()-x1; + sumx1 +=dx; + sum2x1+=dx*dx; + sumxy1+=dx*clusters1[dlayer]->GetY(); + sumy1 +=clusters1[dlayer]->GetY(); + } + if (index2>=0) { + clusters2[dlayer]= (AliTRDcluster*)GetCluster(index2); + ncl++; + sum2++; + Double_t dx = layer2C.GetX()-x2; + sumx2 +=dx; + sum2x2+=dx*dx; + sumxy2+=dx*clusters2[dlayer]->GetY(); + sumy2 +=clusters2[dlayer]->GetY(); + } + } + if (sum1<10) continue; + if (sum2<10) continue; + // + Double_t det1 = sum1*sum2x1-sumx1*sumx1; + Double_t angle1 = (sum1*sumxy1-sumx1*sumy1)/det1; + Double_t pos1 = (sum2x1*sumy1-sumx1*sumxy1)/det1; // at x1 + // + Double_t det2 = sum2*sum2x2-sumx2*sumx2; + Double_t angle2 = (sum2*sumxy2-sumx2*sumy2)/det2; + Double_t pos2 = (sum2x2*sumy2-sumx2*sumxy2)/det2; // at x2 + // + // + + Double_t sumM=0, sumxM=0,sum2xM=0,sumxyM=0, sumyM=0; + // + for (Int_t dlayer=1;dlayer<15;dlayer++){ + clustersM[dlayer]=0; + AliTRDpropagationLayer& layerM=*(fTrSec[ns]->GetLayer(ilayerM-dlayer)); //select propagation layers + Double_t yyM = y1+(tanphi) *(layerM.GetX()-x1); + Double_t zzM = z1+(tantheta)*(layerM.GetX()-x1); + Int_t indexM = layerM.FindNearestCluster(yyM,zzM,road); + if (indexM>=0) { + clustersM[dlayer]= (AliTRDcluster*)GetCluster(indexM); + ncl++; + sumM++; + Double_t dx = layerM.GetX()-xm; + sumxM +=dx; + sum2xM+=dx*dx; + sumxyM+=dx*clustersM[dlayer]->GetY(); + sumyM +=clustersM[dlayer]->GetY(); + } + } + Double_t detM = sumM*sum2xM-sumxM*sumxM; + Double_t posM=0, angleM=0; + if (TMath::Abs(detM)>0.0000001){ + angleM = (sumM*sumxyM-sumxM*sumyM)/detM; + posM = (sum2xM*sumyM-sumxM*sumxyM)/detM; // at xm + } + // + + if (ncl>15){ + TTreeSRedirector& cstream = *fDebugStreamer; + cstream<<"Seeds"<< + "Ncl="<UncheckedAt(iCluster); - if (c->GetNPads()>3&&(iCluster%3>0)) { - delete clusterArray->RemoveAt(iCluster); - continue; - } +// 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); @@ -1939,7 +1993,6 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I // // AliTRDtrackingSector Constructor // - AliTRDpadPlane *padPlane = 0; fGeom = geo; @@ -2097,11 +2150,12 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I } ymax = fGeom->GetChamberWidth(plane)/2.; - // // Modidified for new pad plane class, 22.04.05 (C.B.) // ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.; padPlane = fPar->GetPadPlane(plane,0); ymaxsensitive = (padPlane->GetColSize(1)*padPlane->GetNcols()-4)/2.; + + // ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.; for(Int_t ch = 0; ch < kNchambers; ch++) { zmax[ch] = fGeom->GetChamberLength(plane,ch)/2; @@ -2109,16 +2163,18 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I // Modidified for new pad plane class, 22.04.05 (C.B.) //Float_t pad = fPar->GetRowPadSize(plane,ch,0); Float_t pad = padPlane->GetRowSize(1); + //Float_t pad = fPar->GetRowPadSize(plane,ch,0); Float_t row0 = fPar->GetRow0(plane,ch,0); Int_t nPads = fPar->GetRowMax(plane,ch,0); zmaxsensitive[ch] = Float_t(nPads)*pad/2.; // zc[ch] = (pad * nPads)/2 + row0 - pad/2; - zc[ch] = (pad * nPads)/2 + row0; + // zc[ch] = (pad * nPads)/2 + row0; + zc[ch] = -(pad * nPads)/2 + row0; //zc[ch] = row0+zmax[ch]-AliTRDgeometry::RpadW(); } - dx = fPar->GetDriftVelocity() + dx = fgkDriftCorrection*fPar->GetDriftVelocity() / fPar->GetSamplingFrequency(); rho = 0.00295 * 0.85; radLength = 11.0; @@ -2130,7 +2186,7 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I steps = (Int_t) (dxAmp/dx); for(tb = 0; tb < steps; tb++) { - x = x0 + tb * dx + dx/2; + x = x0 + tb * dx + dx/2+ fgkOffsetX; tbIndex = CookTimeBinIndex(plane, -tb-1); ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex); ppl->SetYmax(ymax,ymaxsensitive); @@ -2148,12 +2204,13 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I InsertLayer(ppl); // Drift region - dx = fPar->GetDriftVelocity() + + dx = fgkDriftCorrection*fPar->GetDriftVelocity() / fPar->GetSamplingFrequency(); steps = (Int_t) (dxDrift/dx); for(tb = 0; tb < steps; tb++) { - x = x0 - tb * dx - dx/2; + x = x0 - tb * dx - dx/2 + fgkOffsetX; //temporary fix - fix it the parameters tbIndex = CookTimeBinIndex(plane, tb); ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex); @@ -2229,7 +2286,8 @@ Int_t AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region - Double_t dx = (Double_t) fPar->GetDriftVelocity() + + Double_t dx = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity() / fPar->GetSamplingFrequency(); Int_t tbAmp = fPar->GetTimeBefore(); @@ -2397,6 +2455,10 @@ Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const return m; } + + + + //______________________________________________________ void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive ) { @@ -2424,7 +2486,7 @@ void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes) -void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters( +Bool_t AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters( Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &radLength, Bool_t &lookForCluster) const { @@ -2467,15 +2529,15 @@ void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters( radLength = 36.66; } }else{ - rho = 2.7; radLength = 24.01; //aluminium in between + cross = kTRUE; rho = 2.7; radLength = 24.01; //aluminium in between } } // - if (fTimeBinIndex>=0) return; + if (fTimeBinIndex>=0) return cross; // // // check hole - if (fHole==kFALSE) return; + if (fHole==kFALSE) return cross; // for(Int_t ch = 0; ch < (Int_t) kZones; ch++) { if (TMath::Abs(z - fZc[ch]) < fZmax[ch]){ @@ -2486,7 +2548,7 @@ void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters( } } } - return; + return cross; } Int_t AliTRDtracker::AliTRDpropagationLayer::GetZone( Double_t z) const @@ -2543,15 +2605,42 @@ Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const { return m; } +Int_t AliTRDtracker::AliTRDpropagationLayer::FindNearestCluster(Double_t y, Double_t z, Double_t maxroad) const +{ + // + // Returns index of the cluster nearest to the given y,z + // + Int_t index = -1; + Int_t maxn = fN; + Double_t mindist = maxroad; + Float_t padlength =-1; + // + for (Int_t i=Find(y-maxroad); iGetSigmaZ2()*12); + } + // + if (c->GetY() > y+maxroad) break; + if((c->GetZ()-z)*(c->GetZ()-z) > padlength*0.75) continue; + if (TMath::Abs(c->GetY()-y)GetY()-y); + index = GetIndex(i); + } + } + return index; +} + + //--------------------------------------------------------- Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) { // // Returns correction factor for tilted pads geometry // - AliTRDpadPlane *padPlane = fPar->GetPadPlane(0,0); Double_t h01 = sin(TMath::Pi() / 180.0 * padPlane->GetTiltingAngle()); + //Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle()); Int_t det = c->GetDetector(); Int_t plane = fGeom->GetPlane(det); @@ -2630,169 +2719,448 @@ 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) +Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters,AliTRDtracklet&tracklet) { // // // 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); + Double_t sigmaz = TMath::Sqrt(track->GetSigmaZ2()); Int_t nall=0; Int_t nfound=0; + Double_t h01 =0; + Int_t plane =-1; + Float_t padlength=0; + AliTRDtrack track2(*track); + Float_t snpy = track->GetSnp(); + Float_t tany = TMath::Sqrt(snpy*snpy/(1.-snpy*snpy)); + if (snpy<0) tany*=-1; + // + Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt()); + Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl()); + Double_t road = 15.*sqrt(track->GetSigmaY2() + sy2); + if (road>6.) road=6.; - for (Int_t it=t0;itGetLayer(it)); + // + for (Int_t it=0;itGetLayer(it+t0)); 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; + track2.PropagateTo(x[it]); + yt[it] = track2.GetY(); + zt[it] = track2.GetZ(); + + Double_t y=yt[it],z=zt[it]; Double_t chi2 =1000000; nall++; // - // find nearest cluster at given pad + // find 2 nearest cluster at given time bin + // + // 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); iGetDetector(); + plane = fGeom->GetPlane(det); + padlength = TMath::Sqrt(c->GetSigmaZ2()*12.); + } + // if (c->GetLocalTimeBin()==0) continue; if (c->GetY() > 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((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue; + + Double_t dist = TMath::Abs(c->GetZ()-z); + if (dist> (0.5*padlength+6.*sigmaz)) continue; // 6 sigma boundary cut + Double_t cost = 0; + // + if (dist> (0.5*padlength-sigmaz)){ // sigma boundary cost function + cost = (dist-0.5*padlength)/(2.*sigmaz); + if (cost>-1) cost= (cost+1.)*(cost+1.); + else cost=0; + } + // Int_t label = TMath::Abs(track->GetLabel()); + // if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue; + chi2=track2.GetPredictedChi2(c,h01)+cost; + // + clfound++; + if (chi2 > maxChi2[1]) continue; + + for (Int_t ih=2;ih<9; ih++){ //store the clusters in the road + if (cl[ih][it]==0){ + cl[ih][it] = c; + indexes[ih][it] =timeBin.GetIndex(i); // index - 9 - reserved for outliers + break; + } } - } - if (cl[0][it]||cl[1][it]) nfound++; + // + if (chi2 t1) continue; + if (!cl[0][it+dt]) continue; + zmean[it]+=cl[0][it+dt]->GetZ(); + nmean[it]+=1.; + } + zmean[it]/=nmean[it]; + } + // + for (Int_t it=0; itGetZ(); - if (TMath::Abs(lastz-cl[ih][it]->GetZ())>1) { - lastz = cl[ih][it]->GetZ(); - changes[ih]++; + Float_t poscor = fgkCoef*(cl[ih][it]->GetLocalTimeBin() - fgkMean)+fgkOffset; + dz[ih][it] = cl[ih][it]->GetZ()- zt[it]; // calculate distance from track in z + dy[ih][it] = cl[ih][it]->GetY()+ dz[ih][it]*h01 - poscor -yt[it]; // in y + } + // minimize changes + if (!cl[0][it]) continue; + if (TMath::Abs(cl[0][it]->GetZ()-zmean[it])> padlength*0.8 &&cl[1][it]) + if (TMath::Abs(cl[1][it]->GetZ()-zmean[it])< padlength*0.5){ + best[0][it]=1; } - 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]; + } + // + // iterative choosing of "best path" + // + // + Int_t label = TMath::Abs(track->GetLabel()); + Int_t bestiter=0; + // + for (Int_t iter=0;iter<9;iter++){ + // + changes[iter]= 0; + sumz = 0; sum=0; sumdy=0;sumdy2=0;sumx=0;sumx2=0;sumxy=0;mpads=0; ngood[iter]=0; nbad[iter]=0; + // linear fit + for (Int_t it=0;itGetZ(); + Double_t zafter = cl[best[iter][it]][it]->GetZ(); + for (Int_t itd = it-1; itd>=0;itd--) { + if (cl[best[iter][itd]][itd]) { + zbefore= cl[best[iter][itd]][itd]->GetZ(); + break; + } + } + for (Int_t itd = it+1; itdGetZ(); + break; + } + } + if (TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore)>0.1&&TMath::Abs(cl[best[iter][it]][it]->GetZ()-zafter)>0.1) changes[iter]++; + // + Double_t dx = x[it]-xmean; // distance to reference x + sumz += cl[best[iter][it]][it]->GetZ(); 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]++; + sumdy += dy[best[iter][it]][it]; + sumdy2+= dy[best[iter][it]][it]*dy[best[iter][it]][it]; + sumx += dx; + sumx2 += dx*dx; + sumxy += dx*dy[best[iter][it]][it]; + mpads += cl[best[iter][it]][it]->GetNPads(); + if (cl[best[iter][it]][it]->GetLabel(0)==label || cl[best[iter][it]][it]->GetLabel(1)==label||cl[best[iter][it]][it]->GetLabel(2)==label){ + ngood[iter]++; } else{ - nbad[ih]++; + nbad[iter]++; } } - if (sumall[ih]>4){ - meanz[ih] = sumz/sum; - meany[ih] = sumdy/sum; - sigmay[ih] = TMath::Sqrt(sumdy2/sum-meany[ih]*meany[ih]); + // + // calculates line parameters + // + Double_t det = sum*sumx2-sumx*sumx; + angle[iter] = (sum*sumxy-sumx*sumdy)/det; + mean[iter] = (sumx2*sumdy-sumx*sumxy)/det; + meanz[iter] = sumz/sum; + moffset[iter] = sumdy/sum; + mpads /= sum; // mean number of pads + // + // + Double_t sigma2 = 0; // normalized residuals - for line fit + Double_t sigma1 = 0; // normalized residuals - constant fit + // + for (Int_t it=0;itfCyy; + Double_t weighty = (moffset[iter]/sigmatr2)/sweight; // weighted mean + Double_t sigmacl = TMath::Sqrt(sigma1*sigma1+track->fCyy); // + Double_t mindist=100000; + Int_t ihbest=0; + for (Int_t ih=0;ih<10;ih++){ + if (!cl[ih][it]) break; + Double_t dist2 = (dy[ih][it]-weighty)/sigmacl; + dist2*=dist2; //chi2 distance + if (dist2quality1){ - 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]; - } + // update best hypothesy if better chi2 according tracklet position and angle + // + Double_t sy2 = smean[iter] + track->fCyy; + Double_t sa2 = sangle[iter] + track->fCee; + Double_t say = track->fCey; + // Double_t chi20 = mean[bestiter]*mean[bestiter]/sy2+angle[bestiter]*angle[bestiter]/sa2; + // Double_t chi21 = mean[iter]*mean[iter]/sy2+angle[iter]*angle[iter]/sa2; + + Double_t detchi = sy2*sa2-say*say; + Double_t invers[3] = {sa2/detchi, sy2/detchi, -say/detchi}; //inverse value of covariance matrix - if (nbad[0]>4){ - nbad[0] = nfound; + Double_t chi20 = mean[bestiter]*mean[bestiter]*invers[0]+angle[bestiter]*angle[bestiter]*invers[1]+ + 2.*mean[bestiter]*angle[bestiter]*invers[2]; + Double_t chi21 = mean[iter]*mean[iter]*invers[0]+angle[iter]*angle[iter]*invers[1]+ + 2*mean[iter]*angle[iter]*invers[2]; + tchi2s[iter] =chi21; + // + if (changes[iter]<=changes[bestiter] && chi2125.) sigma2*=tchi2s[bestiter]/25.; + //if (tchi2s[bestiter]>25.) sigma2=1000.; // dont'accept + + Double_t expectederr = sigma2*sigma2+0.01*0.01; + if (mpads>3.5) expectederr += (mpads-3.5)*0.04; + if (changes[bestiter]>1) expectederr+= changes[bestiter]*0.01; + expectederr+=(0.03*(tany-fgkExB)*(tany-fgkExB))*15; + // if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.; + //expectederr+=10000; + for (Int_t it=0;itGetLocalTimeBin() - fgkMean)+fgkOffset; + cl[best[bestiter][it]][it]->SetSigmaY2(expectederr); // set cluster error + if (!cl[best[bestiter][it]][it]->IsUsed()){ + cl[best[bestiter][it]][it]->SetY( cl[best[bestiter][it]][it]->GetY()-poscor); // ExB corrction correction + cl[best[bestiter][it]][it]->Use(); + } + clusters[it+t0] = indexes[best[bestiter][it]][it]; + } + // + // set tracklet parameters + // + Double_t trackleterr2 = smoffset[bestiter]+0.01*0.01; + if (mpads>3.5) trackleterr2 += (mpads-3.5)*0.04; + trackleterr2+= changes[bestiter]*0.01; + trackleterr2*= TMath::Max(14.-nfound,1.); + trackleterr2+= 0.2*(tany-fgkExB)*(tany-fgkExB); + // + tracklet.Set(xmean, track2.GetY()+moffset[bestiter], meanz[bestiter], track2.GetAlpha(), trackleterr2); //set tracklet parameters + tracklet.SetTilt(h01); + tracklet.SetP0(mean[bestiter]); + tracklet.SetP1(angle[bestiter]); + tracklet.SetN(nfound); + tracklet.SetNCross(changes[bestiter]); + tracklet.SetPlane(plane); + tracklet.SetSigma2(expectederr); + tracklet.SetChi2(tchi2s[bestiter]); + track->fTracklets[plane] = tracklet; + track->fNWrong+=nbad[0]; + // + // Debuging part + // + TTreeSRedirector& cstream = *fDebugStreamer; + AliTRDcluster dummy; + Double_t dy0[100]; + Double_t dyb[100]; + for (Int_t it=0;it= fNclusters) return NULL; @@ -50,9 +56,9 @@ class AliTRDtracker : public AliTracker { Int_t ReadClusters(TObjArray *array, TTree *in) const; Int_t CookSectorIndex(Int_t gs) const { return kTrackingSectors - 1 - gs; } - AliTRDcluster * GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin); + AliTRDcluster * GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin, UInt_t &index); 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); + Int_t FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters, AliTRDtracklet& tracklet); Float_t GetSeedGap() const {return fgkSeedGap;} Int_t GetMaxGap() const {return fMaxGap;} @@ -101,11 +107,12 @@ class AliTRDtracker : public AliTracker { Double_t GetRho() const { return fRho; } Double_t GetX0() const { return fX0; } Int_t GetTimeBinIndex() const { return fTimeBinIndex; } - void GetPropagationParameters(Double_t y, Double_t z, + Bool_t GetPropagationParameters(Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &x0, - Bool_t &lookForCluster) const; + Bool_t &lookForCluster) const; Int_t GetZone( Double_t z) const; Int_t Find(Double_t y) const; + Int_t FindNearestCluster(Double_t y, Double_t z, Double_t maxroad) const; void SetZmax(Int_t cham, Double_t center, Double_t w) { fZc[cham] = center; fZmax[cham] = w; } void SetZ(Double_t* center, Double_t *w, Double_t *wsensitive); @@ -181,6 +188,8 @@ class AliTRDtracker : public AliTracker { protected: + friend class AliTRDtracker::AliTRDtrackingSector; + AliTRDgeometry *fGeom; // Pointer to TRD geometry AliTRDparameter *fPar; // Pointer to TRD parameter @@ -231,6 +240,13 @@ class AliTRDtracker : public AliTracker { static const Double_t fgkSeedErrorSZ; // sz parameter in MakeSeeds static const Float_t fgkLabelFraction; // min fraction of same label static const Float_t fgkWideRoad; // max road width in FindProlongation + // + static const Double_t fgkDriftCorrection; // correction coeficients for drift velocity + static const Double_t fgkOffset; // correction coeficients + static const Double_t fgkOffsetX; // correction coeficients offset in X + static const Double_t fgkCoef; // correction coeficients + static const Double_t fgkMean; // correction coeficients + static const Double_t fgkExB; // correction coeficients Bool_t fVocal; // Whatever... Bool_t fAddTRDseeds; // Something else @@ -249,6 +265,7 @@ class AliTRDtracker : public AliTracker { private: virtual void MakeSeeds(Int_t inner, Int_t outer, Int_t turn); + void MakeSeedsMI(Int_t inner, Int_t outer); Int_t FollowProlongation(AliTRDtrack& t, Int_t rf); Int_t FollowBackProlongation(AliTRDtrack& t); @@ -262,7 +279,7 @@ class AliTRDtracker : public AliTracker { void SetSZ2corr(Float_t w) {fSZ2corr = w;} Double_t ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const; Double_t ExpectedSigmaZ2(Double_t r, Double_t tgl) const; - + TTreeSRedirector *fDebugStreamer; //!debug streamer ClassDef(AliTRDtracker,2) // manager base class }; diff --git a/TRD/TRDrecLinkDef.h b/TRD/TRDrecLinkDef.h index 6deeeb4efb2..35dac45a2b8 100644 --- a/TRD/TRDrecLinkDef.h +++ b/TRD/TRDrecLinkDef.h @@ -21,6 +21,7 @@ #pragma link C++ class AliTRDclusterCorrection+; #pragma link C++ class AliTRDclusterizerMI+; +#pragma link C++ class AliTRDtracklet+; #pragma link C++ class AliTRDtrack+; #pragma link C++ class AliTRDtracker+; -- 2.39.3