#include "AliLog.h"
#include "AliMathBase.h"
-#include "AliTRDseedV1.h"
#include "AliTRDcluster.h"
-#include "AliTRDtrack.h"
+#include "AliTRDseedV1.h"
+#include "AliTRDtrackV1.h"
#include "AliTRDcalibDB.h"
#include "AliTRDchamberTimeBin.h"
#include "AliTRDtrackingChamber.h"
AliTRDseedV1::AliTRDseedV1(Int_t plane)
:AliTRDseed()
,fPlane(plane)
- ,fOwner(kFALSE)
,fMom(0.)
,fSnp(0.)
,fTgl(0.)
AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
:AliTRDseed((AliTRDseed&)ref)
,fPlane(ref.fPlane)
- ,fOwner(kFALSE)
,fMom(ref.fMom)
,fSnp(ref.fSnp)
,fTgl(ref.fTgl)
//
//AliInfo("");
- if(ref.fOwner) SetOwner();
for(int islice=0; islice < knSlices; islice++) fdEdx[islice] = ref.fdEdx[islice];
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = ref.fProb[ispec];
}
//AliInfo(Form("fOwner[%s]", fOwner?"YES":"NO"));
- if(fOwner)
+ if(IsOwner())
for(int itb=0; itb<knTimebins; itb++){
if(!fClusters[itb]) continue;
//AliInfo(Form("deleting c %p @ %d", fClusters[itb], itb));
//____________________________________________________________
-void AliTRDseedV1::Init(AliTRDtrack *track)
+Bool_t AliTRDseedV1::Init(AliTRDtrackV1 *track)
{
// Initialize this tracklet using the track information
//
//
Double_t y, z;
- track->GetProlongation(fX0, y, z);
+ if(!track->GetProlongation(fX0, y, z)) return kFALSE;
fYref[0] = y;
fYref[1] = track->GetSnp()/(1. - track->GetSnp()*track->GetSnp());
fZref[0] = z;
fZref[1] = track->GetTgl();
//printf("Tracklet ref x[%7.3f] y[%7.3f] z[%7.3f], snp[%f] tgl[%f]\n", fX0, fYref[0], fZref[0], track->GetSnp(), track->GetTgl());
+ return kTRUE;
}
Float_t clength = (/*.5 * */AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
AliTRDcluster *cluster = 0x0;
- for(int ic=0; ic<fgTimeBins; ic++){
+ for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
if(!(cluster = fClusters[ic])) continue;
Float_t x = cluster->GetX();
nclusters[slice]++;
} // End of loop over clusters
- // calculate mean charge per slice
- for(int is=0; is<nslices; is++) if(nclusters[is]) fdEdx[is] /= nclusters[is];
+ if(AliTRDReconstructor::RecoParam()->GetPIDMethod() == AliTRDrecoParam::kLQPID){
+ // calculate mean charge per slice (only LQ PID)
+ for(int is=0; is<nslices; is++){
+ if(nclusters[is]) fdEdx[is] /= nclusters[is];
+ }
+ }
}
+
//____________________________________________________________________
Float_t AliTRDseedV1::GetdQdl(Int_t ic) const
{
- return fClusters[ic] ? TMath::Abs(fClusters[ic]->GetQ()) /fdX / TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZfit[1]*fZfit[1]) : 0.;
+ return fClusters[ic] ? TMath::Abs(fClusters[ic]->GetQ()) /fdX / TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]) : 0.;
}
//____________________________________________________________________
return 0x0;
}
+ AliTRDrecoParam *rec = AliTRDReconstructor::RecoParam();
+ if (!rec) {
+ AliError("No TRD reco param.");
+ return 0x0;
+ }
+
// Retrieve the CDB container class with the parametric detector response
- const AliTRDCalPID *pd = calibration->GetPIDObject(AliTRDReconstructor::RecoParam()->GetPIDMethod());
+ const AliTRDCalPID *pd = calibration->GetPIDObject(rec->GetPIDMethod());
if (!pd) {
AliError("No access to AliTRDCalPID object");
return 0x0;
}
-
+ //AliInfo(Form("Method[%d] : %s", AliTRDReconstructor::RecoParam()->GetPIDMethod(), pd->IsA()->GetName()));
+
// calculate tracklet length TO DO
Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
/// TMath::Sqrt((1.0 - fSnp[iPlane]*fSnp[iPlane]) / (1.0 + fTgl[iPlane]*fTgl[iPlane]));
//calculate dE/dx
- CookdEdx(AliTRDReconstructor::RecoParam()->GetNdEdxSlices());
+ CookdEdx(rec->GetNdEdxSlices());
// Sets the a priori probabilities
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
if(!fClusters[ic]) continue;
fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
}
- fOwner = kTRUE;
+ SetBit(1);
} else {
- if(fOwner){
+ if(IsOwner()){
for(int ic=0; ic<knTimebins; ic++){
if(!fClusters[ic]) continue;
delete fClusters[ic];
//fClusters[ic] = tracker->GetClusters(index) TODO
}
}
- fOwner = kFALSE;
+ SetBit(1, kFALSE);
}
}
}
AliTRDchamberTimeBin *layer = 0x0;
- if(AliTRDReconstructor::StreamLevel()>=7 && c){
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel()>=7 && c){
TClonesArray clusters("AliTRDcluster", 24);
clusters.SetOwner(kTRUE);
AliTRDcluster *cc = 0x0;
Int_t det=-1, ncl, ncls = 0;
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+ for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
if(!(layer = chamber->GetTB(iTime))) continue;
if(!(ncl = Int_t(*layer))) continue;
for(int ic=0; ic<ncl; ic++){
// start seed update
for (Int_t iter = 0; iter < niter; iter++) {
ncl = 0;
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+ for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
if(!(layer = chamber->GetTB(iTime))) continue;
if(!Int_t(*layer)) continue;
fZ[iTime] = cl->GetZ();
ncl++;
}
- if(AliTRDReconstructor::StreamLevel()>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fPlane, ncl));
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel()>=7) AliInfo(Form("iter = %d ncl [%d] = %d", iter, fPlane, ncl));
if(ncl>1){
// calculate length of the time bin (calibration aware)
Int_t irp = 0; Float_t x[2]; Int_t tb[2];
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+ for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
if(!fClusters[iTime]) continue;
x[irp] = fClusters[iTime]->GetX();
tb[irp] = iTime;
fdX = (x[1] - x[0]) / (tb[0] - tb[1]);
// update X0 from the clusters (calibration/alignment aware)
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+ for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
if(!(layer = chamber->GetTB(iTime))) continue;
if(!layer->IsT0()) continue;
if(fClusters[iTime]){
fX0 = fClusters[iTime]->GetX();
break;
} else { // we have to infere the position of the anode wire from the other clusters
- for (Int_t jTime = iTime+1; jTime < fgTimeBins; jTime++) {
+ for (Int_t jTime = iTime+1; jTime < AliTRDtrackerV1::GetNTimeBins(); jTime++) {
if(!fClusters[jTime]) continue;
fX0 = fClusters[jTime]->GetX() + fdX * (jTime - iTime);
}
// TODO
// update x reference positions (calibration/alignment aware)
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+ for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
if(!fClusters[iTime]) continue;
fX[iTime] = fClusters[iTime]->GetX() - fX0;
}
AliTRDseed::Update();
}
- if(AliTRDReconstructor::StreamLevel()>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fPlane, fN2));
+ if(AliTRDReconstructor::RecoParam()->GetStreamLevel()>=7) AliInfo(Form("iter = %d nclFit [%d] = %d", iter, fPlane, fN2));
if(IsOK()){
tquality = GetQuality(kZcorr);
// Do cluster projection
AliTRDchamberTimeBin *layer = 0x0;
Int_t nYclusters = 0; Bool_t kEXIT = kFALSE;
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+ for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
if(!(layer = chamber->GetTB(iTime))) continue;
if(!Int_t(*layer)) continue;
// Select only one cluster/TimeBin
Int_t lastCluster = 0;
fN2 = 0;
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+ for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
ncl = tboundary[iTime] - lastCluster;
if(!ncl) continue;
Int_t iptr = lastCluster;
}
// number of minimum numbers of clusters expected for the tracklet
- Int_t kClmin = Int_t(AliTRDReconstructor::RecoParam()->GetFindableClusters()*fgTimeBins);
+ Int_t kClmin = Int_t(AliTRDReconstructor::RecoParam()->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
if (fN2 < kClmin){
AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
fN2 = 0;
// update used clusters
fNUsed = 0;
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
+ for (Int_t iTime = 0; iTime < AliTRDtrackerV1::GetNTimeBins(); iTime++) {
if(!fClusters[iTime]) continue;
if((fClusters[iTime]->IsUsed())) fNUsed++;
}
// 3. Do a Least Square Fit to the data
//
- //Float_t sigmaexp = 0.05 + TMath::Abs(fYref[1] * 0.25); // Expected r.m.s in y direction
- Float_t ycrosscor = fPadLength * fTilt * 0.5; // Y correction for crossing
- Float_t anglecor = fTilt * fZref[1]; // Correction to the angle
+ const Int_t kClmin = 8;
+ const Int_t kNtb = AliTRDtrackerV1::GetNTimeBins();
+ AliTRDtrackerV1::AliTRDLeastSquare fitterY, fitterZ;
+
+ // convertion factor from square to gauss distribution for sigma
+ Double_t convert = 1./TMath::Sqrt(12.);
+
+ // book cluster information
+ Double_t xc[knTimebins+1], yc[knTimebins], zc[knTimebins+1], sy[knTimebins], sz[knTimebins+1];
+ Int_t zRow[knTimebins];
+ AliTRDcluster *c = 0x0;
+ Int_t nc = 0;
+ for (Int_t ic=0; ic<kNtb; ic++) {
+ zRow[ic] = -1;
+ xc[ic] = -1.;
+ yc[ic] = 999.;
+ zc[ic] = 999.;
+ sy[ic] = 0.;
+ sz[ic] = 0.;
+ if(!(c = fClusters[ic])) continue;
+ if(!c->IsInChamber()) continue;
+ Float_t w = 1.;
+ if(c->GetNPads()>4) w = .5;
+ if(c->GetNPads()>5) w = .2;
+ zRow[nc] = c->GetPadRow();
+ xc[nc] = fX0 - c->GetX();
+ yc[nc] = c->GetY();
+ zc[nc] = c->GetZ();
+ sy[ic] = w; // all clusters have the same sigma
+ sz[ic] = fPadLength*convert;
+ fitterZ.AddPoint(&xc[ic], zc[ic], sz[ic]);
+ nc++;
+ }
+ // to few clusters
+ if (nc < kClmin) return kFALSE;
+
- // calculate residuals
- Float_t yres[knTimebins]; // y (r-phi) residuals
- Int_t zint[knTimebins], // Histograming of the z coordinate
- zout[2*knTimebins];//
+ Int_t zN[2*35];
+ Int_t nz = AliTRDtrackerV1::Freq(nc, zRow, zN, kFALSE);
+ // more than one pad row crossing
+ if(nz>2) return kFALSE;
- fN = 0;
- for (Int_t iTime = 0; iTime < fgTimeBins; iTime++) {
- if (!fClusters[iTime]) continue;
- if (!fClusters[iTime]->IsInChamber()) continue;
- yres[iTime] = fY[iTime] - fYref[0] - (fYref[1] + anglecor) * fX[iTime] + fTilt * (fZ[iTime] - fZref[0]);
- zint[fN] = Int_t(fZ[iTime]);
- fN++;
+ // estimate reference parameter at average x
+ Double_t y0 = fYref[0];
+ Double_t dydx = fYref[1];
+ Double_t dzdx = fZref[1];
+ zc[nc] = fZref[0];
+
+ // determine z offset of the fit
+ Int_t nchanges = 0, nCross = 0;
+ if(nz==2){ // tracklet is crossing pad row
+ // Find the break time allowing one chage on pad-rows
+ // with maximal number of accepted clusters
+ Int_t padRef = zRow[0];
+ for (Int_t ic=1; ic<nc; ic++) {
+ if(zRow[ic] == padRef) continue;
+
+ // debug
+ if(zRow[ic-1] == zRow[ic]){
+ printf("ERROR in pad row change!!!\n");
+ }
+
+ // evaluate parameters of the crossing point
+ Float_t sx = (xc[ic-1] - xc[ic])*convert;
+ xc[nc] = .5 * (xc[ic-1] + xc[ic]);
+ zc[nc] = .5 * (zc[ic-1] + zc[ic]);
+ sz[nc] = TMath::Max(dzdx * sx, .01);
+ dzdx = zc[ic-1] > zc[ic] ? 1. : -1.;
+ padRef = zRow[ic];
+ nCross = ic;
+ nchanges++;
+ }
}
- // calculate pad row boundary crosses
- Int_t kClmin = Int_t(AliTRDReconstructor::RecoParam()->GetFindableClusters()*fgTimeBins);
- Int_t nz = AliMathBase::Freq(fN, zint, zout, kFALSE);
- fZProb = zout[0];
- if(nz <= 1) zout[3] = 0;
- if(zout[1] + zout[3] < kClmin) {
- AliWarning(Form("Not enough clusters to fit the cross boundary tracklet %d [%d].", zout[1]+zout[3], kClmin));
+ // condition on nCross and reset nchanges TODO
+
+ if(nchanges==1){
+ if(dzdx * fZref[1] < 0.){
+ AliInfo("tracklet direction does not correspond to the track direction. TODO.");
+ }
+ SetBit(2, kTRUE); // mark pad row crossing
+ fCross[0] = xc[nc]; fCross[2] = zc[nc]; fCross[3] = sz[nc];
+ fitterZ.AddPoint(&xc[nc], zc[nc], sz[nc]);
+ fitterZ.Eval();
+ dzdx = fZref[1]; // we don't trust Parameter[1] ??;
+ zc[nc] = fitterZ.GetFunctionParameter(0);
+ } else if(nchanges > 1){ // debug
+ AliInfo("ERROR in n changes!!!");
return kFALSE;
}
- // Z distance bigger than pad - length
- if (TMath::Abs(zout[0]-zout[2]) > fPadLength) zout[3]=0;
-
-
- Double_t sumw = 0.,
- sumwx = 0.,
- sumwx2 = 0.,
- sumwy = 0.,
- sumwxy = 0.,
- sumwz = 0.,
- sumwxz = 0.;
- Int_t npads;
- fMPads = 0;
- fMeanz = 0.;
- // we will use only the clusters which are in the detector range
- for(int iTime=0; iTime<fgTimeBins; iTime++){
- fUsable[iTime] = kFALSE;
- if (!fClusters[iTime]) continue;
- npads = fClusters[iTime]->GetNPads();
-
- fUsable[iTime] = kTRUE;
- fN2++;
- fMPads += npads;
- Float_t weight = 1.0;
- if(npads > 5) weight = 0.2;
- else if(npads > 4) weight = 0.5;
- sumw += weight;
- sumwx += fX[iTime] * weight;
- sumwx2 += fX[iTime] * fX[iTime] * weight;
- sumwy += weight * yres[iTime];
- sumwxy += weight * yres[iTime] * fX[iTime];
- sumwz += weight * fZ[iTime];
- sumwxz += weight * fZ[iTime] * fX[iTime];
+
+
+ // estimate deviation from reference direction
+ dzdx *= fTilt;
+ for (Int_t ic=0; ic<nc; ic++) {
+ yc[ic] -= y0 + xc[ic]*(dydx + dzdx) + fTilt * (zc[ic] - zc[nc]);
+ fitterY.AddPoint(&xc[ic], yc[ic], sy[ic]);
}
- if (fN2 < kClmin){
- AliWarning(Form("Not enough clusters to fit the tracklet %d [%d].", fN2, kClmin));
- fN2 = 0;
- return kFALSE;
- }
- fMeanz = sumwz / sumw;
- fNChange = 0;
-
- // Tracklet on boundary
- Float_t correction = 0;
- if (fNChange > 0) {
- if (fMeanz < fZProb) correction = ycrosscor;
- if (fMeanz > fZProb) correction = -ycrosscor;
- }
+ fitterY.Eval();
+ fYfit[0] = y0+fitterY.GetFunctionParameter(0);
+ fYfit[1] = dydx+fitterY.GetFunctionParameter(1);
+ if(nchanges) fCross[1] = fYfit[0] + fCross[0] * fYfit[1];
- Double_t det = sumw * sumwx2 - sumwx * sumwx;
- fYfitR[0] = (sumwx2 * sumwy - sumwx * sumwxy) / det;
- fYfitR[1] = (sumw * sumwxy - sumwx * sumwy) / det;
-
- fSigmaY2 = 0;
- for (Int_t i = 0; i < fgTimeBins+1; i++) {
- if (!fUsable[i]) continue;
- Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
- fSigmaY2 += delta*delta;
- }
- fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
-
- fZfitR[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
- fZfitR[1] = (sumw * sumwxz - sumwx * sumwz) / det;
- fZfit[0] = (sumwx2 * sumwz - sumwx * sumwxz) / det;
- fZfit[1] = (sumw * sumwxz - sumwx * sumwz) / det;
- fYfitR[0] += fYref[0] + correction;
- fYfitR[1] += fYref[1];
- fYfit[0] = fYfitR[0];
- fYfit[1] = fYfitR[1];
+// printf("\nnz = %d\n", nz);
+// for(int ic=0; ic<35; ic++) printf("%d row[%d]\n", ic, zRow[ic]);
+//
+// for(int ic=0; ic<nz; ic++) printf("%d n[%d]\n", ic, zN[ic]);
return kTRUE;
}
+//___________________________________________________________________
+void AliTRDseedV1::Draw(Option_t*)
+{
+}
//___________________________________________________________________
-void AliTRDseedV1::Print()
+void AliTRDseedV1::Print(Option_t*) const
{
//
// Printing the seedstatus
//
- AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
- Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
printf("Seed status :\n");
printf(" fTilt = %f\n", fTilt);
printf(" fPadLength = %f\n", fPadLength);
printf(" fX0 = %f\n", fX0);
- for(int ic=0; ic<nTimeBins; ic++) {
+ for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++) {
const Char_t *isUsable = fUsable[ic]?"Yes":"No";
printf(" %d X[%f] Y[%f] Z[%f] Indexes[%d] clusters[%p] usable[%s]\n"
, ic
printf(" fCC =%f\n",fCC);
printf(" fChi2 =%f\n", fChi2);
printf(" fChi2Z =%f\n", fChi2Z);
-
}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff