///////////////////////////////////////////////////////////////////////////////
// //
-// The standard TRD tracker //
+// The standard TRD tracker //
+// Based on Kalman filltering approach //
// //
///////////////////////////////////////////////////////////////////////////////
#include <TObjArray.h>
#include "AliTRDgeometry.h"
-#include "AliTRDparameter.h"
#include "AliTRDpadPlane.h"
-#include "AliTRDgeometryFull.h"
+#include "AliTRDgeometry.h"
#include "AliTRDcluster.h"
#include "AliTRDtrack.h"
+#include "AliTRDseed.h"
#include "AliESD.h"
#include "AliTRDcalibDB.h"
#include "AliRieman.h"
#include "AliTrackPointArray.h"
#include "AliAlignObj.h"
-
-//
+#include "AliTRDReconstructor.h"
ClassImp(AliTRDtracker)
-ClassImp(AliTRDseed)
-
-
const Float_t AliTRDtracker::fgkMinClustersInTrack = 0.5;
const Float_t AliTRDtracker::fgkLabelFraction = 0.8;
const Double_t AliTRDtracker::fgkMaxChi2 = 12.;
- const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // correspond to tan = 3
- const Double_t AliTRDtracker::fgkMaxStep = 2.; // maximal step size in propagation
-
-
-//
-
+ const Double_t AliTRDtracker::fgkMaxSnp = 0.95; // correspond to tan = 3
+ const Double_t AliTRDtracker::fgkMaxStep = 2.; // maximal step size in propagation
-
-
-//____________________________________________________________________
+//_____________________________________________________________________________
AliTRDtracker::AliTRDtracker():AliTracker(),
fGeom(0),
fNclusters(0),
fAddTRDseeds(kFALSE),
fNoTilt(kFALSE)
{
+ //
// Default constructor
+ //
for(Int_t i=0;i<kTrackingSectors;i++) fTrSec[i]=0;
for(Int_t j=0;j<5;j++)
- for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
+ for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
fDebugStreamer = 0;
+
}
-//____________________________________________________________________
+
+//_____________________________________________________________________________
AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
{
//
}
else {
printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
- fGeom = new AliTRDgeometryFull();
- fGeom->SetPHOShole();
- fGeom->SetRICHhole();
+ fGeom = new AliTRDgeometry();
}
+ fGeom->ReadGeoMatrices();
savedir->cd();
-
fNclusters = 0;
fClusters = new TObjArray(2000);
fNseeds = 0;
fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
savedir->cd();
+
}
-//___________________________________________________________________
+//_____________________________________________________________________________
AliTRDtracker::~AliTRDtracker()
{
//
//fDebugStreamer->Close();
delete fDebugStreamer;
}
-}
-
-//_____________________________________________________________________
+}
-Int_t AliTRDtracker::LocalToGlobalID(Int_t lid){
+//_____________________________________________________________________________
+Int_t AliTRDtracker::LocalToGlobalID(Int_t lid)
+{
//
- // transform internal TRD ID to global detector ID
+ // Transform internal TRD ID to global detector ID
//
+
Int_t isector = fGeom->GetSector(lid);
Int_t ichamber= fGeom->GetChamber(lid);
Int_t iplan = fGeom->GetPlane(lid);
};
Int_t modId = isector*fGeom->Ncham()+ichamber;
UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
+
return volid;
+
}
-Int_t AliTRDtracker::GlobalToLocalID(Int_t gid){
+//_____________________________________________________________________________
+Int_t AliTRDtracker::GlobalToLocalID(Int_t gid)
+{
//
- // transform global detector ID to local detector ID
+ // Transform global detector ID to local detector ID
//
+
Int_t modId=0;
AliAlignObj::ELayerID layerId = AliAlignObj::VolUIDToLayer(gid, modId);
Int_t isector = modId/fGeom->Ncham();
}
if (iLayer<0) return -1;
Int_t lid = fGeom->GetDetector(iLayer,ichamber,isector);
+
return lid;
-}
+}
-Bool_t AliTRDtracker::Transform(AliTRDcluster * cluster){
+//_____________________________________________________________________________
+Bool_t AliTRDtracker::Transform(AliTRDcluster * cluster)
+{
+ //
+ // Transform something ... whatever ...
//
+
+ const Double_t kX0shift = 2.52; // magic constants for geo manager transformation
+ const Double_t kX0shift5 = 3.05; //
//
- const Double_t kDriftCorrection = 1.01; // drift coeficient correction
- const Double_t kExBcor = 0.001; // ExB coef correction
- const Double_t kTime0Cor = 0.32; // time0 correction
//
// apply alignment and calibration to transform cluster
//
//
+ Int_t detector = cluster->GetDetector();
+ Int_t plane = fGeom->GetPlane(cluster->GetDetector());
+ Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
+ Int_t sector = fGeom->GetSector(cluster->GetDetector());
+
Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.); // drift distance
//
- Int_t plane = fGeom->GetPlane(cluster->GetDetector());
- Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
- cluster->SetX(xplane- kDriftCorrection*(cluster->GetX()-kTime0Cor));
- //
// ExB correction
//
Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
//
- cluster->SetY(cluster->GetY() - driftX*(exB+ kExBcor));
+ AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
+ AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
+ Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
+ Double_t localPos[3], localPosTracker[3];
+ localPos[0] = -cluster->GetX();
+ localPos[1] = cluster->GetY() - driftX*exB;
+ localPos[2] = cluster->GetZ() -zshiftIdeal;
+ //
+ cluster->SetY(cluster->GetY() - driftX*exB);
+ Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
+ cluster->SetX(xplane- cluster->GetX());
+ //
+ TGeoHMatrix * matrix = fGeom->GetCorrectionMatrix(cluster->GetDetector());
+ if (!matrix){
+ // no matrix found - if somebody used geometry with holes
+ AliError("Invalid Geometry - Default Geometry used\n");
+ return kTRUE;
+ }
+ matrix->LocalToMaster(localPos, localPosTracker);
+ //
+ //
+ //
+ if (AliTRDReconstructor::StreamLevel()>1){
+ (*fDebugStreamer)<<"Transform"<<
+ "Cl.="<<cluster<<
+ "matrix.="<<matrix<<
+ "Detector="<<detector<<
+ "Sector="<<sector<<
+ "Plane="<<plane<<
+ "Chamber="<<chamber<<
+ "lx0="<<localPosTracker[0]<<
+ "ly0="<<localPosTracker[1]<<
+ "lz0="<<localPosTracker[2]<<
+ "\n";
+ }
+ //
+ if (plane==5)
+ cluster->SetX(localPosTracker[0]+kX0shift5);
+ else
+ cluster->SetX(localPosTracker[0]+kX0shift);
+
+ cluster->SetY(localPosTracker[1]);
+ cluster->SetZ(localPosTracker[2]);
+
return kTRUE;
+
}
-Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
+//_____________________________________________________________________________
+// Bool_t AliTRDtracker::Transform(AliTRDcluster * cluster)
+//{
+// //
+// //
+// const Double_t kDriftCorrection = 1.01; // drift coeficient correction
+// const Double_t kTime0Cor = 0.32; // time0 correction
+// //
+// const Double_t kX0shift = 2.52;
+// const Double_t kX0shift5 = 3.05;
+
+// //
+// // apply alignment and calibration to transform cluster
+// //
+// //
+// Int_t detector = cluster->GetDetector();
+// Int_t plane = fGeom->GetPlane(cluster->GetDetector());
+// Int_t chamber = fGeom->GetChamber(cluster->GetDetector());
+// Int_t sector = fGeom->GetSector(cluster->GetDetector());
+
+// Double_t dxAmp = (Double_t) fGeom->CamHght(); // Amplification region
+// Double_t driftX = TMath::Max(cluster->GetX()-dxAmp*0.5,0.); // drift distance
+// //
+// // ExB correction
+// //
+// Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
+// Double_t exB = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
+// //
+
+// AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
+// AliTRDpadPlane * padPlane = commonParam->GetPadPlane(plane,chamber);
+// Double_t zshiftIdeal = 0.5*(padPlane->GetRow0()+padPlane->GetRowEnd());
+// Double_t localPos[3], globalPos[3], localPosTracker[3], localPosTracker2[3];
+// localPos[2] = -cluster->GetX();
+// localPos[0] = cluster->GetY() - driftX*exB;
+// localPos[1] = cluster->GetZ() -zshiftIdeal;
+// TGeoHMatrix * matrix = fGeom->GetGeoMatrix(cluster->GetDetector());
+// matrix->LocalToMaster(localPos, globalPos);
+
+// Double_t sectorAngle = 20.*(sector%18)+10;
+// TGeoHMatrix rotSector;
+// rotSector.RotateZ(sectorAngle);
+// rotSector.LocalToMaster(globalPos, localPosTracker);
+// //
+// //
+// TGeoHMatrix matrix2(*matrix);
+// matrix2.MultiplyLeft(&rotSector);
+// matrix2.LocalToMaster(localPos,localPosTracker2);
+// //
+// //
+// //
+// cluster->SetY(cluster->GetY() - driftX*exB);
+// Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane);
+// cluster->SetX(xplane- kDriftCorrection*(cluster->GetX()-kTime0Cor));
+// (*fDebugStreamer)<<"Transform"<<
+// "Cl.="<<cluster<<
+// "matrix.="<<matrix<<
+// "matrix2.="<<&matrix2<<
+// "Detector="<<detector<<
+// "Sector="<<sector<<
+// "Plane="<<plane<<
+// "Chamber="<<chamber<<
+// "lx0="<<localPosTracker[0]<<
+// "ly0="<<localPosTracker[1]<<
+// "lz0="<<localPosTracker[2]<<
+// "lx2="<<localPosTracker2[0]<<
+// "ly2="<<localPosTracker2[1]<<
+// "lz2="<<localPosTracker2[2]<<
+// "\n";
+// //
+// if (plane==5)
+// cluster->SetX(localPosTracker[0]+kX0shift5);
+// else
+// cluster->SetX(localPosTracker[0]+kX0shift);
+
+// cluster->SetY(localPosTracker[1]);
+// cluster->SetZ(localPosTracker[2]);
+// return kTRUE;
+// }
+
+//_____________________________________________________________________________
+Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track)
+{
//
// Rotates the track when necessary
//
}
return kTRUE;
-}
+}
-AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin, UInt_t &index){
+//_____________________________________________________________________________
+AliTRDcluster *AliTRDtracker::GetCluster(AliTRDtrack *track, Int_t plane
+ , Int_t timebin, UInt_t &index)
+{
//
- //try to find cluster in the backup list
+ // Try to find cluster in the backup list
//
+
AliTRDcluster * cl =0;
- UInt_t *indexes = track->GetBackupIndexes();
+ Int_t *indexes = track->GetBackupIndexes();
for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
if (indexes[i]==0) break;
AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
break;
}
}
+
return cl;
-}
+}
-Int_t AliTRDtracker::GetLastPlane(AliTRDtrack * track){
+//_____________________________________________________________________________
+Int_t AliTRDtracker::GetLastPlane(AliTRDtrack * track)
+{
+ //
+ // Return last updated plane
//
- //return last updated plane
+
Int_t lastplane=0;
- UInt_t *indexes = track->GetBackupIndexes();
+ Int_t *indexes = track->GetBackupIndexes();
for (UInt_t i=0;i<kMaxTimeBinIndex;i++){
AliTRDcluster * cli = (AliTRDcluster*)fClusters->UncheckedAt(indexes[i]);
if (!cli) break;
lastplane = iplane;
}
}
+
return lastplane;
+
}
-//___________________________________________________________________
+
+//_____________________________________________________________________________
Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
{
//
Float_t foundMin = fgkMinClustersInTrack * timeBins;
Int_t nseed = 0;
Int_t found = 0;
- Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
+ // Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
Int_t n = event->GetNumberOfTracks();
for (Int_t i=0; i<n; i++) {
//seed2->ResetCovariance();
AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
AliTRDtrack &t=*pt;
- FollowProlongation(t, innerTB);
+ FollowProlongation(t);
if (t.GetNumberOfClusters() >= foundMin) {
UseClusters(&t);
CookLabel(pt, 1-fgkLabelFraction);
cout<<"Total number of found tracks: "<<found<<endl;
return 0;
+
}
-
-
//_____________________________________________________________________________
-Int_t AliTRDtracker::PropagateBack(AliESD* event) {
+Int_t AliTRDtracker::PropagateBack(AliESD* event)
+{
//
// Gets seeds from ESD event. The seeds are AliTPCtrack's found and
// backpropagated by the TPC tracker. Each seed is first propagated
Float_t p4 = track->GetC();
//
Int_t expectedClr = FollowBackProlongation(*track);
- /*
- // only debug purpose
- if (track->GetNumberOfClusters()<expectedClr/3){
- AliTRDtrack *track1 = new AliTRDtrack(*seed);
- track1->SetSeedLabel(lbl);
- FollowBackProlongation(*track1);
- AliTRDtrack *track2= new AliTRDtrack(*seed);
- track->SetSeedLabel(lbl);
- FollowBackProlongation(*track2);
- delete track1;
- delete track2;
- }
- */
if (TMath::Abs(track->GetC()-p4)/TMath::Abs(p4)<0.2 || TMath::Abs(track->GetPt())>0.8 ) {
//
//make backup for back propagation
// Debug part of tracking
TTreeSRedirector& cstream = *fDebugStreamer;
Int_t eventNr = event->GetEventNumber();
- if (track->GetBackupTrack()){
- cstream<<"Tracks"<<
- "EventNr="<<eventNr<<
- "ESD.="<<seed<<
- "trd.="<<track<<
- "trdback.="<<track->GetBackupTrack()<<
- "\n";
- }else{
- cstream<<"Tracks"<<
- "EventNr="<<eventNr<<
- "ESD.="<<seed<<
- "trd.="<<track<<
- "trdback.="<<track<<
- "\n";
+ if (AliTRDReconstructor::StreamLevel()>0){
+ if (track->GetBackupTrack()){
+ cstream<<"Tracks"<<
+ "EventNr="<<eventNr<<
+ "ESD.="<<seed<<
+ "trd.="<<track<<
+ "trdback.="<<track->GetBackupTrack()<<
+ "\n";
+ }else{
+ cstream<<"Tracks"<<
+ "EventNr="<<eventNr<<
+ "ESD.="<<seed<<
+ "trd.="<<track<<
+ "trdback.="<<track<<
+ "\n";
+ }
}
//
//Propagation to the TOF (I.Belikov)
if (track->PropagateTo(xtof)) {
seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
- for (Int_t i=0;i<kNPlane;i++) {
- seed->SetTRDsignals(track->GetPIDsignals(i),i);
- seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
+ for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
+ for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
+ seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
+ }
+ seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
}
// seed->SetTRDtrack(new AliTRDtrack(*track));
if (track->GetNumberOfClusters()>foundMin) found++;
if (track->GetNumberOfClusters()>15&&track->GetNumberOfClusters()>0.5*expectedClr){
seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
//seed->SetStatus(AliESDtrack::kTRDStop);
- for (Int_t i=0;i<kNPlane;i++) {
- seed->SetTRDsignals(track->GetPIDsignals(i),i);
- seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
+ for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
+ for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
+ seed->SetTRDsignals(track->GetPIDsignals(i,j),i,j);
+ }
+ seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
}
//seed->SetTRDtrack(new AliTRDtrack(*track));
found++;
cerr<<"Number of seeds: "<<fNseeds<<endl;
cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
- // MakeSeedsMI(3,5,event); //new seeding
-
+ if (AliTRDReconstructor::SeedingOn()) MakeSeedsMI(3,5,event); //new seeding
fSeeds->Clear(); fNseeds=0;
delete [] index;
Float_t foundMin = fgkMinClustersInTrack * timeBins;
Int_t nseed = 0;
Int_t found = 0;
- Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
+ // Int_t innerTB = fTrSec[0]->GetInnerTimeBin();
AliTRDtrack seed2;
Int_t n = event->GetNumberOfTracks();
// }
AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
- UInt_t * indexes2 = seed2.GetIndexes();
- for (Int_t i=0;i<kNPlane;i++) {
- pt->SetPIDsignals(seed2.GetPIDsignals(i),i);
+ Int_t * indexes2 = seed2.GetIndexes();
+ for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
+ for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
+ pt->SetPIDsignals(seed2.GetPIDsignals(i,j),i,j);
+ }
pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
}
- UInt_t * indexes3 = pt->GetBackupIndexes();
+ Int_t * indexes3 = pt->GetBackupIndexes();
for (Int_t i=0;i<200;i++) {
if (indexes2[i]==0) break;
indexes3[i] = indexes2[i];
}
//AliTRDtrack *pt = seed2;
AliTRDtrack &t=*pt;
- FollowProlongation(t, innerTB);
+ FollowProlongation(t);
if (t.GetNumberOfClusters() >= foundMin) {
// UseClusters(&t);
//CookLabel(pt, 1-fgkLabelFraction);
Double_t xTPC = 250;
if(PropagateToX(t,xTPC,fgkMaxStep)) {
seed->UpdateTrackParams(pt, AliESDtrack::kTRDrefit);
- for (Int_t i=0;i<kNPlane;i++) {
- seed->SetTRDsignals(pt->GetPIDsignals(i),i);
+ for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
+ for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
+ seed->SetTRDsignals(pt->GetPIDsignals(i,j),i,j);
+ }
seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
}
}else{
pt2->CookdEdx( ); // Modification by PS
CookdEdxTimBin(*pt2);
seed->UpdateTrackParams(pt2, AliESDtrack::kTRDrefit);
- for (Int_t i=0;i<kNPlane;i++) {
- seed->SetTRDsignals(pt2->GetPIDsignals(i),i);
+ for (Int_t i=0;i<AliESDtrack::kNPlane;i++) {
+ for (Int_t j=0;j<AliESDtrack::kNSlice;j++) {
+ seed->SetTRDsignals(pt2->GetPIDsignals(i,j),i,j);
+ }
seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
}
}
}
-
-
-
-//---------------------------------------------------------------------------
-Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
+//_____________________________________________________________________________
+Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t)
{
+ //
// Starting from current position on track=t this function tries
// to extrapolate the track up to timeBin=0 and to confirm prolongation
// if a close cluster is found. Returns the number of clusters
// expected to be found in sensitive layers
// GeoManager used to estimate mean density
+ //
+
Int_t sector;
Int_t lastplane = GetLastPlane(&t);
Double_t radLength = 0.0;
Int_t expectedNumberOfClusters = 0;
//
//
- Double_t alpha=AliTRDgeometry::GetAlpha();
- Double_t tanmax = TMath::Tan(0.5*alpha);
-
- for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr>rf; nr--) {
+ //
+ for (Int_t iplane = lastplane; iplane>=0; iplane--){
//
+ Int_t row0 = GetGlobalTimeBin(0, iplane,GetTimeBinsPerPlane()-1);
+ Int_t rowlast = GetGlobalTimeBin(0, iplane,0);
//
- Int_t currentplane = fTrSec[0]->GetLayer(nr)->GetPlane();
- Double_t currentx = fTrSec[0]->GetLayer(nr)->GetX();
+ // propagate track close to the plane if neccessary
+ //
+ Double_t currentx = fTrSec[0]->GetLayer(rowlast)->GetX();
if (currentx < -fgkMaxStep +t.GetX()){
//propagate closer to chamber - safety space fgkMaxStep
if (!PropagateToX(t, currentx+fgkMaxStep, fgkMaxStep)) break;
}
if (!AdjustSector(&t)) break;
+ //
+ // get material budget
+ //
Double_t xyz0[3],xyz1[3],param[7],x,y,z;
t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //starting global position
+ // end global position
+ x = fTrSec[0]->GetLayer(row0)->GetX();
+ if (!t.GetProlongation(x,y,z)) break;
+ xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
+ xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
+ xyz1[2] = z;
+ AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+ rho = param[0];
+ radLength = param[1]; // get mean propagation parameters
//
+ // propagate nad update
//
- // propagate and update track in active layers
- //
- Int_t nr0 = nr; //first active layer
- if (nr >rf && (fTrSec[0]->GetLayer(nr)->IsSensitive())){
- //
- // get all time bins at given plane
- //
- while (nr >rf && ((fTrSec[0]->GetLayer(nr)->IsSensitive())) && fTrSec[0]->GetLayer(nr)->GetPlane() == currentplane){
- x = fTrSec[0]->GetLayer(nr)->GetX();
- nr--;
- if (!t.GetProlongation(x,y,z)) break;
- if (TMath::Abs(y)>x*tanmax){
- nr--;
- break;
- }
- }
- nr++;
- x = fTrSec[0]->GetLayer(nr)->GetX();
- if (!t.GetProlongation(x,y,z)) break;
- xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
- xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
- xyz1[2] = z;
- // end global position
- AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
- rho = param[0];
- radLength = param[1]; // get mean propagation parameters
- }
- //
- // propagate and update
- if (nr0-nr< fTimeBinsPerPlane/2 ){
- // short tracklet - do not update - edge effect
- continue;
- }
sector = t.GetSector();
- //
- //
- for (Int_t ilayer=nr0;ilayer>=nr;ilayer--) {
+ // for (Int_t itime=GetTimeBinsPerPlane()-1;itime>=0;itime--) {
+ for (Int_t itime=0 ;itime<GetTimeBinsPerPlane();itime++) {
+ Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
expectedNumberOfClusters++;
t.fNExpected++;
if (t.fX>345) t.fNExpectedLast++;
AliTRDcluster *cl=0;
UInt_t index=0;
Double_t maxChi2=fgkMaxChi2;
- //dx = (fTrSec[sector]->GetLayer(ilayer+1))->GetX()-timeBin.GetX();
x = timeBin.GetX();
- // t.PropagateTo(x,radLength,rho);
if (timeBin) {
AliTRDcluster * cl0 = timeBin[0];
if (!cl0) continue; // no clusters in given time bin
cl =cl2;
Double_t h01 = GetTiltFactor(cl);
maxChi2=t.GetPredictedChi2(cl,h01);
- }
-
+ }
if (cl) {
// if (cl->GetNPads()<5)
Double_t dxsample = timeBin.GetdX();
Double_t xcluster = cl->GetX();
t.PropagateTo(xcluster,radLength,rho);
if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
- if(!t.Update(cl,maxChi2,index,h01)) {
- }
- }
+ }
}
}
- }
+ }
}
- return expectedNumberOfClusters;
-
-
-}
-
+ return expectedNumberOfClusters;
+}
-//___________________________________________________________________
+//_____________________________________________________________________________
Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
{
-
+ //
// Starting from current radial position of track <t> this function
// extrapolates the track up to outer timebin and in the sensitive
// layers confirms prolongation if a close cluster is found.
// Returns the number of clusters expected to be found in sensitive layers
// Use GEO manager for material Description
+ //
Int_t sector;
Int_t clusters[1000];
for (Int_t i=0;i<1000;i++) clusters[i]=-1;
- Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
Double_t radLength = 0.0;
Double_t rho = 0.0;
- Double_t x;
Int_t expectedNumberOfClusters = 0;
- x = t.GetX();
-
- Double_t alpha=AliTRDgeometry::GetAlpha(); // note: change in meaning
- Double_t tanmax = TMath::Tan(0.5*alpha);
- Int_t nr;
Float_t ratio0=0;
AliTRDtracklet tracklet;
//
//
-
- for (nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB; nr++) {
+ for (Int_t iplane = 0; iplane<AliESDtrack::kNPlane; iplane++){
+ Int_t row0 = GetGlobalTimeBin(0, iplane,GetTimeBinsPerPlane()-1);
+ Int_t rowlast = GetGlobalTimeBin(0, iplane,0);
+ //
+ Double_t currentx = fTrSec[0]->GetLayer(row0)->GetX();
+ if (currentx<t.GetX()) continue;
//
- // propagate to current X
+ // propagate closer to chamber if neccessary
//
- Int_t currentplane = fTrSec[0]->GetLayer(nr)->GetPlane();
- Double_t currentx = fTrSec[0]->GetLayer(nr)->GetX();
if (currentx > fgkMaxStep +t.GetX()){
- //propagate closter to chamber
if (!PropagateToX(t, currentx-fgkMaxStep, fgkMaxStep)) break;
}
if (!AdjustSector(&t)) break;
if (TMath::Abs(t.GetSnp())>fgkMaxSnp) break;
- Double_t xyz0[3],xyz1[3],param[7],x,y,z;
- t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //starting global position
- //
- //
- //
- Int_t nr0 = nr;
- if (nr <outerTB && (fTrSec[0]->GetLayer(nr)->IsSensitive())){
- //
- // get all time bins at given plane
- //
- while (nr <outerTB && fTrSec[0]->GetLayer(nr)->GetPlane() == currentplane){
- x = fTrSec[0]->GetLayer(nr)->GetX();
- nr++;
- if (!t.GetProlongation(x,y,z)) break;
- if (TMath::Abs(y)>(x*tanmax)){
- nr++;
- break;
- }
- }
- nr--;
- //
- //
- //
- x = fTrSec[0]->GetLayer(nr)->GetX();
- if (!t.GetProlongation(x,y,z)) break;
- // minimal mean and maximal budget scan
- Float_t minbudget =10000;
- Float_t meanbudget =0;
- Float_t maxbudget =-1;
- // Float_t normbudget =0;
- // for (Int_t idy=-1;idy<=1;idy++)
- // for (Int_t idz=-1;idz<=1;idz++){
- for (Int_t idy=0;idy<1;idy++)
- for (Int_t idz=0;idz<1;idz++){
- Double_t y2 = y+idy*TMath::Min(TMath::Sqrt(t.GetSigmaY2()),1.);
- Double_t z2 = z+idz*TMath::Min(TMath::Sqrt(t.GetSigmaZ2()),1.);
-
- xyz1[0] = x*TMath::Cos(t.GetAlpha())-y2*TMath::Sin(t.GetAlpha());
- xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y2*TMath::Cos(t.GetAlpha());
- xyz1[2] = z2;
- AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
- Float_t budget = param[0]*param[4];
- meanbudget+=budget;
- if (budget<minbudget) minbudget=budget;
- if (budget>maxbudget) maxbudget=budget;
- }
- t.fBudget[0]+=minbudget;
- t.fBudget[1]+=meanbudget/9.;
- t.fBudget[2]+=minbudget;
- //
- xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
- xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
- xyz1[2] = z;
- // end global position
- AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
- rho = param[0];
- radLength = param[1]; // get mean propagation parameters
- }
- //
//
+ // get material budget inside of chamber
//
- if (nr-nr0< fTimeBinsPerPlane/2){
- // short tracklet - do not update - edge effect
- continue;
- }
+ Double_t xyz0[3],xyz1[3],param[7],x,y,z;
+ t.GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]); //starting global position
+ // end global position
+ x = fTrSec[0]->GetLayer(rowlast)->GetX();
+ if (!t.GetProlongation(x,y,z)) break;
+ xyz1[0] = x*TMath::Cos(t.GetAlpha())-y*TMath::Sin(t.GetAlpha());
+ xyz1[1] = +x*TMath::Sin(t.GetAlpha())+y*TMath::Cos(t.GetAlpha());
+ xyz1[2] = z;
+ AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+ rho = param[0];
+ radLength = param[1]; // get mean propagation parameters
//
+ // Find clusters
//
sector = t.GetSector();
- Float_t ncl = FindClusters(sector,nr0,nr,&t,clusters,tracklet);
- if (tracklet.GetN()-2*tracklet.GetNCross()<10) continue;
+ Float_t ncl = FindClusters(sector,row0,rowlast,&t,clusters,tracklet);
+ if (tracklet.GetN()<GetTimeBinsPerPlane()/3) continue;
//
+ // Propagate and update track
//
- for (Int_t ilayer=nr0;ilayer<=nr;ilayer++) {
+ for (Int_t itime= GetTimeBinsPerPlane()-1;itime>=0;itime--) {
+ Int_t ilayer = GetGlobalTimeBin(0, iplane,itime);
expectedNumberOfClusters++;
t.fNExpected++;
if (t.fX>345) t.fNExpectedLast++;
if(!t.Update(cl,maxChi2,index,h01)) {
}
}
- //
-
-// 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
-// }
-// }
+ //
// reset material budget if 2 consecutive gold
if (plane>0)
if (t.fTracklets[plane].GetN()+t.fTracklets[plane-1].GetN()>20){
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 (nr<outerTB)
- t.SetStop(kTRUE);
- else
- t.SetStop(kFALSE);
- return expectedNumberOfClusters;
-}
-
-
-
+ return expectedNumberOfClusters;
+}
+//_____________________________________________________________________________
Int_t AliTRDtracker::PropagateToX(AliTRDtrack& t, Double_t xToGo, Double_t maxStep)
{
+ //
// Starting from current radial position of track <t> this function
// extrapolates the track up to radial position <xToGo>.
// Returns 1 if track reaches the plane, and 0 otherwise
+ //
+
const Double_t kEpsilon = 0.00001;
// Double_t tanmax = TMath::Tan(0.5*AliTRDgeometry::GetAlpha());
Double_t xpos = t.GetX();
AdjustSector(&t);
xpos = t.GetX();
}
+
return 1;
}
-
-
//_____________________________________________________________________________
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"<<endl;
if (ReadClusters(fClusters,cTree)) {
Error("LoadClusters","Problem with reading the clusters !");
Transform(c);
fTrSec[trackingSector]->GetLayer(layer)->InsertCluster(c,index);
}
+
return 0;
+
}
//_____________________________________________________________________________
}
-//__________________________________________________________________________
+//_____________________________________________________________________________
void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD * esd)
{
//
// Creates seeds using clusters between position inner plane and outer plane
//
- const Double_t maxtheta = 1;
- const Double_t maxphi = 2.0;
+
+ const Double_t kMaxTheta = 1;
+ const Double_t kMaxPhi = 2.0;
//
const Double_t kRoad0y = 6; // road for middle cluster
const Double_t kRoad0z = 8.5; // road for middle cluster
//
const Double_t kRoad2y = 3; // road in y for extrapolated cluster
const Double_t kRoad2z = 20; // road in z for extrapolated cluster
- const Int_t maxseed = 3000;
+ const Int_t kMaxSeed = 3000;
Int_t maxSec=AliTRDgeometry::kNsect;
//
//
//
// registered seed
- AliTRDseed *pseed = new AliTRDseed[maxseed*6];
- AliTRDseed *seed[maxseed];
- for (Int_t iseed=0;iseed<maxseed;iseed++) seed[iseed]= &pseed[iseed*6];
+ AliTRDseed *pseed = new AliTRDseed[kMaxSeed*6];
+ AliTRDseed *seed[kMaxSeed];
+ for (Int_t iseed=0;iseed<kMaxSeed;iseed++) seed[iseed]= &pseed[iseed*6];
AliTRDseed *cseed = seed[0];
//
- Double_t seedquality[maxseed];
- Double_t seedquality2[maxseed];
- Double_t seedparams[maxseed][7];
- Int_t seedlayer[maxseed];
+ Double_t seedquality[kMaxSeed];
+ Double_t seedquality2[kMaxSeed];
+ Double_t seedparams[kMaxSeed][7];
+ Int_t seedlayer[kMaxSeed];
Int_t registered =0;
- Int_t sort[maxseed];
+ Int_t sort[kMaxSeed];
//
// seeding part
//
padlength[sLayer+3] = TMath::Sqrt(cl3->GetSigmaZ2()*12.);
ycl[sLayer+3] = cl3->GetY();
zcl[sLayer+3] = cl3->GetZ();
- Float_t yymin0 = ycl[sLayer+3] - 1- maxphi *(xcl[sLayer+3]-xcl[sLayer+0]);
- Float_t yymax0 = ycl[sLayer+3] + 1+ maxphi *(xcl[sLayer+3]-xcl[sLayer+0]);
+ Float_t yymin0 = ycl[sLayer+3] - 1- kMaxPhi *(xcl[sLayer+3]-xcl[sLayer+0]);
+ Float_t yymax0 = ycl[sLayer+3] + 1+ kMaxPhi *(xcl[sLayer+3]-xcl[sLayer+0]);
Int_t maxn0 = layer0; //
for (Int_t icl0=layer0.Find(yymin0);icl0<maxn0;icl0++){
AliTRDcluster *cl0 = layer0[icl0];
zcl[sLayer+0] = cl0->GetZ();
if ( ycl[sLayer+0]>yymax0) break;
Double_t tanphi = (ycl[sLayer+3]-ycl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
- if (TMath::Abs(tanphi)>maxphi) continue;
+ if (TMath::Abs(tanphi)>kMaxPhi) continue;
Double_t tantheta = (zcl[sLayer+3]-zcl[sLayer+0])/(xcl[sLayer+3]-xcl[sLayer+0]);
- if (TMath::Abs(tantheta)>maxtheta) continue;
+ if (TMath::Abs(tantheta)>kMaxTheta) continue;
padlength[sLayer+0] = TMath::Sqrt(cl0->GetSigmaZ2()*12.);
//
// expected position in 1 layer
if (cl0->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
if (cl1->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
if (cl2->GetLabel(0)!=cl3->GetLabel(0)) isFake = kTRUE;
- if ((!isFake) || (icl3%10)==0 ){ //debugging print
- TTreeSRedirector& cstream = *fDebugStreamer;
- cstream<<"Seeds0"<<
- "isFake="<<isFake<<
- "Cl0.="<<cl0<<
- "Cl1.="<<cl1<<
- "Cl2.="<<cl2<<
- "Cl3.="<<cl3<<
- "Xref="<<xref<<
- "X0="<<xcl[sLayer+0]<<
- "X1="<<xcl[sLayer+1]<<
- "X2="<<xcl[sLayer+2]<<
- "X3="<<xcl[sLayer+3]<<
- "Y2exp="<<y2exp<<
- "Z2exp="<<z2exp<<
- "Chi2R="<<chi2R<<
- "Chi2Z="<<chi2Z<<
- "Seed0.="<<&cseed[sLayer+0]<<
- "Seed1.="<<&cseed[sLayer+1]<<
- "Seed2.="<<&cseed[sLayer+2]<<
- "Seed3.="<<&cseed[sLayer+3]<<
- "Zmin="<<minmax[0]<<
- "Zmax="<<minmax[1]<<
- "\n";
+ if (AliTRDReconstructor::StreamLevel()>0){
+ if ((!isFake) || (icl3%10)==0 ){ //debugging print
+ TTreeSRedirector& cstream = *fDebugStreamer;
+ cstream<<"Seeds0"<<
+ "isFake="<<isFake<<
+ "Cl0.="<<cl0<<
+ "Cl1.="<<cl1<<
+ "Cl2.="<<cl2<<
+ "Cl3.="<<cl3<<
+ "Xref="<<xref<<
+ "X0="<<xcl[sLayer+0]<<
+ "X1="<<xcl[sLayer+1]<<
+ "X2="<<xcl[sLayer+2]<<
+ "X3="<<xcl[sLayer+3]<<
+ "Y2exp="<<y2exp<<
+ "Z2exp="<<z2exp<<
+ "Chi2R="<<chi2R<<
+ "Chi2Z="<<chi2Z<<
+ "Seed0.="<<&cseed[sLayer+0]<<
+ "Seed1.="<<&cseed[sLayer+1]<<
+ "Seed2.="<<&cseed[sLayer+2]<<
+ "Seed3.="<<&cseed[sLayer+3]<<
+ "Zmin="<<minmax[0]<<
+ "Zmax="<<minmax[1]<<
+ "\n";
+ }
}
//<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
2.*TMath::Abs(tseed.fMeanz-tseed.fZref[0])/padlength[jLayer];
- if (iter==0 && tseed.isOK()) {
+ if (iter==0 && tseed.IsOK()) {
cseed[sLayer+jLayer] = tseed;
quality = tquality;
if (tquality<5) break;
}
- if (tseed.isOK() && tquality<quality)
+ if (tseed.IsOK() && tquality<quality)
cseed[sLayer+jLayer] = tseed;
}
- if (!cseed[sLayer+jLayer].isOK()){
+ if (!cseed[sLayer+jLayer].IsOK()){
isOK = kFALSE;
break;
}
if (!isOK) continue;
nclusters=0;
for (Int_t iLayer=0;iLayer<4;iLayer++){
- if (cseed[sLayer+iLayer].isOK()){
+ if (cseed[sLayer+iLayer].IsOK()){
nclusters+=cseed[sLayer+iLayer].fN2;
}
}
//
for (Int_t iLayer=0;iLayer<2;iLayer++){
Int_t jLayer = tLayer[iLayer]; // set tracking layer
- if ( (jLayer==0) && !(cseed[1].isOK())) continue; // break not allowed
- if ( (jLayer==5) && !(cseed[4].isOK())) continue; // break not allowed
+ if ( (jLayer==0) && !(cseed[1].IsOK())) continue; // break not allowed
+ if ( (jLayer==5) && !(cseed[4].IsOK())) continue; // break not allowed
Float_t zexp = cseed[jLayer].fZref[0];
Double_t zroad = padlength[jLayer]*0.5+1.;
//
tseed.fZ[iTime] = cl->GetZ(); // register cluster
}
tseed.Update();
- if (tseed.isOK()){
+ if (tseed.IsOK()){
Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
Float_t tquality = (18-tseed.fN2)/2. + TMath::Abs(dangle)/0.1+
TMath::Abs(tseed.fYfit[0]-tseed.fYref[0])/0.2+
}
zroad*=2.;
}
- if ( cseed[jLayer].isOK()){
+ if ( cseed[jLayer].IsOK()){
cseed[jLayer].CookLabels();
cseed[jLayer].UpdateUsed();
nusedf+= cseed[jLayer].fNUsed;
Float_t squality[6];
Int_t sortindexes[6];
for (Int_t jLayer=0;jLayer<6;jLayer++){
- if (bseed[jLayer].isOK()){
+ if (bseed[jLayer].IsOK()){
AliTRDseed &tseed = bseed[jLayer];
Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
tseed.fZ[iTime] = cl->GetZ(); // register cluster
}
tseed.Update();
- if (tseed.isOK()) {
+ if (tseed.IsOK()) {
Float_t dangle = tseed.fYfit[1]-tseed.fYref[1];
Double_t zcor = tseed.fTilt*(tseed.fZProb-tseed.fZref[0]);
//
for (Int_t iLayer=0;iLayer<6;iLayer++) {
if (TMath::Abs(cseed[iLayer].fYref[0]/cseed[iLayer].fX0)<0.15)
findable++;
- if (cseed[iLayer].isOK()){
+ if (cseed[iLayer].IsOK()){
nclusters+=cseed[iLayer].fN2;
nlayers++;
}
if (nlayers<3) continue;
rieman.Reset();
for (Int_t iLayer=0;iLayer<6;iLayer++){
- if (cseed[iLayer].isOK()) rieman.AddPoint(xcl[iLayer],cseed[iLayer].fYfitR[0],
+ if (cseed[iLayer].IsOK()) rieman.AddPoint(xcl[iLayer],cseed[iLayer].fYfitR[0],
cseed[iLayer].fZProb,1,10);
}
rieman.Update();
chi2RF =0;
chi2ZF =0;
for (Int_t iLayer=0;iLayer<6;iLayer++){
- if (cseed[iLayer].isOK()){
+ if (cseed[iLayer].IsOK()){
cseed[iLayer].fYref[0] = rieman.GetYat(xcl[iLayer]);
chi2RF += (cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0])*
(cseed[iLayer].fYref[0]-cseed[iLayer].fYfitR[0]);
fitterT2.ClearPoints();
rieman2.Reset();
for (Int_t iLayer=0; iLayer<6;iLayer++){
- if (!cseed[iLayer].isOK()) continue;
+ if (!cseed[iLayer].IsOK()) continue;
for (Int_t itime=0;itime<25;itime++){
if (!cseed[iLayer].fUsable[itime]) continue;
Double_t x = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2; // x relative to the midle chamber
//
Bool_t acceptablez =kTRUE;
for (Int_t iLayer=0; iLayer<6;iLayer++){
- if (cseed[iLayer].isOK()){
+ if (cseed[iLayer].IsOK()){
Double_t zT2 = rpolz0+rpolz1*(xcl[iLayer] - xref2);
if (TMath::Abs(cseed[iLayer].fZProb-zT2)>padlength[iLayer]*0.5+1)
acceptablez = kFALSE;
//
Double_t aC = fitterTC.GetParameter(0);
Double_t bC = fitterTC.GetParameter(1);
- Double_t CC = aC/TMath::Sqrt(bC*bC+1.); // curvature
+ Double_t cC = aC/TMath::Sqrt(bC*bC+1.); // curvature
//
Double_t aR = fitterT2.GetParameter(0);
Double_t bR = fitterT2.GetParameter(1);
Double_t dR = fitterT2.GetParameter(2);
- Double_t CR = 1+bR*bR-dR*aR;
+ Double_t cR = 1+bR*bR-dR*aR;
Double_t dca = 0.;
- if (CR>0){
+ if (cR>0){
dca = -dR/(TMath::Sqrt(1+bR*bR-dR*aR)+TMath::Sqrt(1+bR*bR));
- CR = aR/TMath::Sqrt(CR);
+ cR = aR/TMath::Sqrt(cR);
}
//
Double_t chi2ZT2=0, chi2ZTC=0;
for (Int_t iLayer=0; iLayer<6;iLayer++){
- if (cseed[iLayer].isOK()){
+ if (cseed[iLayer].IsOK()){
Double_t zT2 = rpolz0+rpolz1*(xcl[iLayer] - xref2);
Double_t zTC = polz0c+polz1c*(xcl[iLayer] - xref2);
chi2ZT2 += TMath::Abs(cseed[iLayer].fMeanz-zT2);
AliTRDseed::FitRiemanTilt(cseed, kTRUE);
Float_t sumdaf = 0;
for (Int_t iLayer=0;iLayer<6;iLayer++){
- if (cseed[iLayer].isOK())
+ if (cseed[iLayer].IsOK())
sumdaf += TMath::Abs((cseed[iLayer].fYfit[1]-cseed[iLayer].fYref[1])/cseed[iLayer].fSigmaY2);
}
sumdaf /= Float_t (nlayers-2.);
// if (nlayers == findable && TMath::Log(0.000000001+seedquality2[index])<-4.) isGold =kTRUE; // gold
// if (isGold &&nusedf<10){
// for (Int_t jLayer=0;jLayer<6;jLayer++){
-// if ( seed[index][jLayer].isOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
+// if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.1)
// seed[index][jLayer].UseClusters(); //sign gold
// }
// }
//
//
Int_t index0=0;
- if (!cseed[0].isOK()){
+ if (!cseed[0].IsOK()){
index0 = 1;
- if (!cseed[1].isOK()) index0 = 2;
+ if (!cseed[1].IsOK()) index0 = 2;
}
seedparams[registered][0] = cseed[index0].fX0;
seedparams[registered][1] = cseed[index0].fYref[0];
seedparams[registered][2] = cseed[index0].fZref[0];
- seedparams[registered][5] = CR;
- seedparams[registered][3] = cseed[index0].fX0*CR - TMath::Sin(TMath::ATan(cseed[0].fYref[1]));
+ seedparams[registered][5] = cR;
+ seedparams[registered][3] = cseed[index0].fX0*cR - TMath::Sin(TMath::ATan(cseed[0].fYref[1]));
seedparams[registered][4] = cseed[index0].fZref[1]/
TMath::Sqrt(1+cseed[index0].fYref[1]*cseed[index0].fYref[1]);
seedparams[registered][6] = ns;
Int_t labels[12], outlab[24];
Int_t nlab=0;
for (Int_t iLayer=0;iLayer<6;iLayer++){
- if (!cseed[iLayer].isOK()) continue;
+ if (!cseed[iLayer].IsOK()) continue;
if (cseed[iLayer].fLabels[0]>=0) {
labels[nlab] = cseed[iLayer].fLabels[0];
nlab++;
Int_t frequency = outlab[1];
for (Int_t iLayer=0;iLayer<6;iLayer++){
cseed[iLayer].fFreq = frequency;
- cseed[iLayer].fC = CR;
- cseed[iLayer].fCC = CC;
+ cseed[iLayer].fC = cR;
+ cseed[iLayer].fCC = cC;
cseed[iLayer].fChi2 = chi2TR;
cseed[iLayer].fChi2Z = chi2ZF;
}
if (1||(!isFake)){ //debugging print
Float_t zvertex = GetZ();
TTreeSRedirector& cstream = *fDebugStreamer;
+ if (AliTRDReconstructor::StreamLevel()>0)
cstream<<"Seeds1"<<
"isFake="<<isFake<<
"Vertex="<<zvertex<<
"C="<<curv<< // non constrained - no tilt correction
"DR="<<dR<< // DR parameter - tilt correction
"DCA="<<dca<< // DCA - tilt correction
- "CR="<<CR<< // non constrained curvature - tilt correction
- "CC="<<CC<< // constrained curvature
+ "CR="<<cR<< // non constrained curvature - tilt correction
+ "CC="<<cC<< // constrained curvature
"Polz0="<<polz0c<<
"Polz1="<<polz1c<<
"RPolz0="<<rpolz0<<
"sLayer="<<sLayer<<
"\n";
}
- if (registered<maxseed-1) {
+ if (registered<kMaxSeed-1) {
registered++;
cseed = seed[registered];
}
// choos best
//
TMath::Sort(registered,seedquality2,sort,kTRUE);
- Bool_t signedseed[maxseed];
+ Bool_t signedseed[kMaxSeed];
for (Int_t i=0;i<registered;i++){
signedseed[i]= kFALSE;
}
for (Int_t jLayer=0;jLayer<6;jLayer++){
if (TMath::Abs(seed[index][jLayer].fYref[0]/xcl[jLayer])<0.15)
findable++;
- if (seed[index][jLayer].isOK()){
+ if (seed[index][jLayer].IsOK()){
seed[index][jLayer].UpdateUsed();
ncl +=seed[index][jLayer].fN2;
nused +=seed[index][jLayer].fNUsed;
Int_t labels[1000], outlab[1000];
Int_t nlab=0;
for (Int_t iLayer=0;iLayer<6;iLayer++){
- if (seed[index][iLayer].isOK()){
+ if (seed[index][iLayer].IsOK()){
if (seed[index][iLayer].fLabels[0]>=0) {
labels[nlab] = seed[index][iLayer].fLabels[0];
nlab++;
Float_t ratio = Float_t(nused)/Float_t(ncl);
if (ratio<0.25){
for (Int_t jLayer=0;jLayer<6;jLayer++){
- if ( seed[index][jLayer].isOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.2 )
+ if ( seed[index][jLayer].IsOK()&&TMath::Abs(seed[index][jLayer].fYfit[1]-seed[index][jLayer].fYfit[1])<0.2 )
seed[index][jLayer].UseClusters(); //sign gold
}
}
esdtrack.SetLabel(label);
esd->AddTrack(&esdtrack);
TTreeSRedirector& cstream = *fDebugStreamer;
- cstream<<"Tracks"<<
- "EventNr="<<eventNr<<
- "ESD.="<<&esdtrack<<
- "trd.="<<track<<
- "trdback.="<<track<<
- "\n";
+ if (AliTRDReconstructor::StreamLevel()>0)
+ cstream<<"Tracks"<<
+ "EventNr="<<eventNr<<
+ "ESD.="<<&esdtrack<<
+ "trd.="<<track<<
+ "trdback.="<<track<<
+ "\n";
}
-
- cstream<<"Seeds2"<<
+ if (AliTRDReconstructor::StreamLevel()>0)
+ cstream<<"Seeds2"<<
"Iter="<<iter<<
"Track.="<<track<<
"Like="<<seedquality[index]<<
}
}
} // end of loop over sectors
+
delete [] pseed;
+
}
//_____________________________________________________________________________
// from the file. The names of the cluster tree and branches
// should match the ones used in AliTRDclusterizer::WriteClusters()
//
+
Int_t nsize = Int_t(ClusterTree->GetTotBytes()/(sizeof(AliTRDcluster)));
TObjArray *clusterArray = new TObjArray(nsize+1000);
delete clusterArray;
return 0;
+
}
-//__________________________________________________________________
+//_____________________________________________________________________________
Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint& p) const
{
//
}
-//__________________________________________________________________
+//_____________________________________________________________________________
void AliTRDtracker::CookLabel(AliKalmanTrack* pt, Float_t wrong) const
{
//
}
-
-//__________________________________________________________________
+//_____________________________________________________________________________
void AliTRDtracker::UseClusters(const AliKalmanTrack* t, Int_t from) const
{
//
// Use clusters, but don't abuse them!
//
+
const Float_t kmaxchi2 =18;
const Float_t kmincl =10;
AliTRDtrack * track = (AliTRDtrack*)t;
if (track->fTracklets[iplane].GetN()<kmincl) continue;
if (!(c->IsUsed())) c->Use();
}
-}
+}
-//_____________________________________________________________________
+//_____________________________________________________________________________
Double_t AliTRDtracker::ExpectedSigmaY2(Double_t , Double_t , Double_t ) const
{
+ //
// Parametrised "expected" error of the cluster reconstruction in Y
+ //
Double_t s = 0.08 * 0.08;
return s;
+
}
-//_____________________________________________________________________
+//_____________________________________________________________________________
Double_t AliTRDtracker::ExpectedSigmaZ2(Double_t , Double_t ) const
{
+ //
// Parametrised "expected" error of the cluster reconstruction in Z
+ //
Double_t s = 9 * 9 /12.;
return s;
+
}
-//_____________________________________________________________________
+//_____________________________________________________________________________
Double_t AliTRDtracker::GetX(Int_t sector, Int_t plane, Int_t localTB) const
{
//
}
-
-//_______________________________________________________
-AliTRDtracker::AliTRDpropagationLayer::AliTRDpropagationLayer(Double_t x,
- Double_t dx, Double_t rho, Double_t radLength, Int_t tbIndex, Int_t plane)
+//_____________________________________________________________________________
+AliTRDtracker::AliTRDpropagationLayer
+ ::AliTRDpropagationLayer(Double_t x, Double_t dx, Double_t rho
+ , Double_t radLength, Int_t tbIndex, Int_t plane)
{
//
// AliTRDpropagationLayer constructor
}
-//_______________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer::SetHole(
- Double_t Zmax, Double_t Ymax, Double_t rho,
- Double_t radLength, Double_t Yc, Double_t Zc)
+//_____________________________________________________________________________
+void AliTRDtracker::AliTRDpropagationLayer
+ ::SetHole(Double_t Zmax, Double_t Ymax, Double_t rho
+ , Double_t radLength, Double_t Yc, Double_t Zc)
{
//
// Sets hole in the layer
//
+
fHole = kTRUE;
fHoleZc = Zc;
fHoleZmax = Zmax;
fHoleYmax = Ymax;
fHoleRho = rho;
fHoleX0 = radLength;
+
}
-
-//_______________________________________________________
-AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs)
+//_____________________________________________________________________________
+AliTRDtracker::AliTRDtrackingSector
+ ::AliTRDtrackingSector(AliTRDgeometry* geo, Int_t gs)
{
//
// AliTRDtrackingSector Constructor
//
+
AliTRDpadPlane *padPlane = 0;
fGeom = geo;
}
-//______________________________________________________
-
-Int_t AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t localTB) const
+//_____________________________________________________________________________
+Int_t AliTRDtracker::AliTRDtrackingSector
+ ::CookTimeBinIndex(Int_t plane, Int_t localTB) const
{
//
// depending on the digitization parameters calculates "global"
// time bin index for timebin <localTB> in plane <plane>
//
//
+
Int_t tbPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
- Int_t gtb = (plane+1) * tbPerPlane - localTB;
+ Int_t gtb = (plane+1) * tbPerPlane - localTB -1;
if (localTB<0) return -1;
- if (gtb<0) return -1;
+ if (gtb<0) return -1;
+
return gtb;
-}
-//______________________________________________________
+}
-void AliTRDtracker::AliTRDtrackingSector::MapTimeBinLayers()
+//_____________________________________________________________________________
+void AliTRDtracker::AliTRDtrackingSector
+ ::MapTimeBinLayers()
{
//
// For all sensitive time bins sets corresponding layer index
fTimeBinIndex[index] = i;
}
- Double_t x1, dx1, x2, dx2, gap;
-
- for(Int_t i = 0; i < fN-1; i++) {
- x1 = fLayers[i]->GetX();
- dx1 = fLayers[i]->GetdX();
- x2 = fLayers[i+1]->GetX();
- dx2 = fLayers[i+1]->GetdX();
- gap = (x2 - dx2/2) - (x1 + dx1/2);
-// if(gap < -0.01) {
-// printf("*** warning: layers %d and %d are overlayed:\n",i,i+1);
-// printf(" %f + %f + %f > %f\n", x1, dx1/2, dx2/2, x2);
-// }
-// if(gap > 0.01) {
-// printf("*** warning: layers %d and %d have a large gap:\n",i,i+1);
-// printf(" (%f - %f) - (%f + %f) = %f\n",
-// x2, dx2/2, x1, dx1, gap);
-// }
- }
}
-
-
-//______________________________________________________
-
-Int_t AliTRDtracker::AliTRDtrackingSector::GetLayerNumber(Double_t x) const
+//_____________________________________________________________________________
+Int_t AliTRDtracker::AliTRDtrackingSector
+ ::GetLayerNumber(Double_t x) const
{
//
// Returns the number of time bin which in radial position is closest to <x>
}
if(TMath::Abs(x - fLayers[m]->GetX()) >
TMath::Abs(x - fLayers[m+1]->GetX())) return m+1;
+
else return m;
}
-//______________________________________________________
-
-Int_t AliTRDtracker::AliTRDtrackingSector::GetInnerTimeBin() const
+//_____________________________________________________________________________
+Int_t AliTRDtracker::AliTRDtrackingSector
+ ::GetInnerTimeBin() const
{
//
// Returns number of the innermost SENSITIVE propagation layer
//
return GetLayerNumber(0);
-}
-//______________________________________________________
+}
-Int_t AliTRDtracker::AliTRDtrackingSector::GetOuterTimeBin() const
+//_____________________________________________________________________________
+Int_t AliTRDtracker::AliTRDtrackingSector
+ ::GetOuterTimeBin() const
{
//
// Returns number of the outermost SENSITIVE time bin
//
return GetLayerNumber(GetNumberOfTimeBins() - 1);
-}
-//______________________________________________________
+}
-Int_t AliTRDtracker::AliTRDtrackingSector::GetNumberOfTimeBins() const
+//_____________________________________________________________________________
+Int_t AliTRDtracker::AliTRDtrackingSector
+ ::GetNumberOfTimeBins() const
{
//
// Returns number of SENSITIVE time bins
layer = GetLayerNumber(tb);
if(layer>=0) break;
}
+
return tb+1;
-}
-//______________________________________________________
+}
-void AliTRDtracker::AliTRDtrackingSector::InsertLayer(AliTRDpropagationLayer* pl)
+//_____________________________________________________________________________
+void AliTRDtracker::AliTRDtrackingSector
+ ::InsertLayer(AliTRDpropagationLayer* pl)
{
//
// Insert layer <pl> in fLayers array.
// Layers are sorted according to X coordinate.
+ //
if ( fN == ((Int_t) kMaxLayersPerSector)) {
printf("AliTRDtrackingSector::InsertLayer(): Too many layers !\n");
}
-//______________________________________________________
-
-Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
+//_____________________________________________________________________________
+Int_t AliTRDtracker::AliTRDtrackingSector
+ ::Find(Double_t x) const
{
//
// Returns index of the propagation layer nearest to X
if (x > fLayers[m]->GetX()) b=m+1;
else e=m;
}
- return m;
-}
-
-
+ return m;
+}
-//______________________________________________________
-void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
+//_____________________________________________________________________________
+void AliTRDtracker::AliTRDpropagationLayer
+ ::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
{
//
// set centers and the width of sectors
+ //
+
for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
fZc[icham] = center[icham];
fZmax[icham] = w[icham];
fZmaxSensitive[icham] = wsensitive[icham];
- // printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
+ // printf("chamber\t%d\tzc\t%f\tzmax\t%f\tzsens\t%f\n",icham,fZc[icham],fZmax[icham],fZmaxSensitive[icham]);
}
+
}
-//______________________________________________________
+//_____________________________________________________________________________
void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
{
//
// set centers and the width of sectors
+ //
+
fHole = kFALSE;
for (Int_t icham=0;icham< AliTRDgeometry::kNcham;icham++){
fIsHole[icham] = holes[icham];
if (holes[icham]) fHole = kTRUE;
}
-}
-
+}
-
-
-//______________________________________________________
-
-void AliTRDtracker::AliTRDpropagationLayer::InsertCluster(AliTRDcluster* c,
- UInt_t index) {
-
-// Insert cluster in cluster array.
-// Clusters are sorted according to Y coordinate.
+//_____________________________________________________________________________
+void AliTRDtracker::AliTRDpropagationLayer
+ ::InsertCluster(AliTRDcluster* c, UInt_t index)
+{
+ //
+ // Insert cluster in cluster array.
+ // Clusters are sorted according to Y coordinate.
+ //
if(fTimeBinIndex < 0) {
printf("*** attempt to insert cluster into non-sensitive time bin!\n");
memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTRDcluster*));
memmove(fIndex +i+1 ,fIndex +i,(fN-i)*sizeof(UInt_t));
fIndex[i]=index; fClusters[i]=c; fN++;
-}
-
-//______________________________________________________
-Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const {
+}
-// Returns index of the cluster nearest in Y
+//_____________________________________________________________________________
+Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Float_t y) const
+{
+ //
+ // Returns index of the cluster nearest in Y
+ //
if (fN<=0) return 0;
if (y <= fClusters[0]->GetY()) return 0;
if (y > fClusters[m]->GetY()) b=m+1;
else e=m;
}
+
return m;
+
}
-Int_t AliTRDtracker::AliTRDpropagationLayer::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad, Float_t maxroadz) const
+//_____________________________________________________________________________
+Int_t AliTRDtracker::AliTRDpropagationLayer
+ ::FindNearestCluster(Float_t y, Float_t z, Float_t maxroad
+ , Float_t maxroadz) const
{
//
// Returns index of the cluster nearest to the given y,z
//
+
Int_t index = -1;
Int_t maxn = fN;
Float_t mindist = maxroad;
index = fIndex[i];
}
}
+
return index;
-}
+}
-//---------------------------------------------------------
+//_____________________________________________________________________________
+Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c)
+{
+ //
+ // Returns correction factor for tilted pads geometry
+ //
-Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
-//
-// Returns correction factor for tilted pads geometry
-//
Int_t det = c->GetDetector();
Int_t plane = fGeom->GetPlane(det);
AliTRDpadPlane *padPlane = AliTRDCommonParam::Instance()->GetPadPlane(plane,0);
Double_t h01 = TMath::Tan(-TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
if(fNoTilt) h01 = 0;
+
return h01;
-}
+}
+//_____________________________________________________________________________
void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack)
{
+ //
// *** ADDED TO GET MORE INFORMATION FOR TRD PID ---- PS
// This is setting fdEdxPlane and fTimBinPlane
// Sums up the charge in each plane for track TRDtrack and also get the
// Time bin for Max. Cluster
// Prashant Shukla (shukla@physi.uni-heidelberg.de)
+ //
+
+ Double_t clscharge[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
+ Double_t maxclscharge[AliESDtrack::kNPlane];
+ Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
+ Int_t timebin[AliESDtrack::kNPlane];
- Double_t clscharge[kNPlane], maxclscharge[kNPlane];
- Int_t nCluster[kNPlane], timebin[kNPlane];
+ //Initialization of cluster charge per plane.
+ for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
+ for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
+ clscharge[iPlane][iSlice] = 0.0;
+ nCluster[iPlane][iSlice] = 0;
+ }
+ }
//Initialization of cluster charge per plane.
- for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
- clscharge[iPlane] = 0.0;
- nCluster[iPlane] = 0;
+ for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
timebin[iPlane] = -1;
maxclscharge[iPlane] = 0.0;
}
for (Int_t iClus = 0; iClus < nClus; iClus++) {
Double_t charge = TRDtrack.GetClusterdQdl(iClus);
Int_t index = TRDtrack.GetClusterIndex(iClus);
- AliTRDcluster *TRDcluster = (AliTRDcluster *) GetCluster(index);
- if (!TRDcluster) continue;
- Int_t tb = TRDcluster->GetLocalTimeBin();
+ AliTRDcluster *pTRDcluster = (AliTRDcluster *) GetCluster(index);
+ if (!pTRDcluster) continue;
+ Int_t tb = pTRDcluster->GetLocalTimeBin();
if (!tb) continue;
- Int_t detector = TRDcluster->GetDetector();
+ Int_t detector = pTRDcluster->GetDetector();
Int_t iPlane = fGeom->GetPlane(detector);
- clscharge[iPlane] = clscharge[iPlane]+charge;
+ Int_t iSlice = tb*AliESDtrack::kNSlice/AliTRDtrack::kNtimeBins;
+ clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice]+charge;
if(charge > maxclscharge[iPlane]) {
maxclscharge[iPlane] = charge;
timebin[iPlane] = tb;
}
- nCluster[iPlane]++;
+ nCluster[iPlane][iSlice]++;
} // end of loop over cluster
// Setting the fdEdxPlane and fTimBinPlane variabales
- Double_t Total_ch = 0;
- for (Int_t iPlane = 0; iPlane < kNPlane; iPlane++) {
- // Quality control of TRD track.
- if (nCluster[iPlane]<= 5) {
- clscharge[iPlane]=0.0;
- timebin[iPlane]=-1;
+ Double_t totalCharge = 0;
+
+ for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
+ for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
+ if (nCluster[iPlane][iSlice]) clscharge[iPlane][iSlice] /= nCluster[iPlane][iSlice];
+ TRDtrack.SetPIDsignals(clscharge[iPlane][iSlice], iPlane, iSlice);
+ totalCharge= totalCharge+clscharge[iPlane][iSlice];
}
- if (nCluster[iPlane]) clscharge[iPlane] /= nCluster[iPlane];
- TRDtrack.SetPIDsignals(clscharge[iPlane], iPlane);
- TRDtrack.SetPIDTimBin(timebin[iPlane], iPlane);
- Total_ch= Total_ch+clscharge[iPlane];
+ TRDtrack.SetPIDTimBin(timebin[iPlane], iPlane);
}
+
// Int_t i;
// Int_t nc=TRDtrack.GetNumberOfClusters();
// Float_t dedx=0;
// TRDtrack.SetPIDTimBin(timbin[iPlane], iPlane);
// }
-} // end of function
-
+}
-Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters,AliTRDtracklet&tracklet)
+//_____________________________________________________________________________
+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
+ // Try to find nearest clusters to the track in timebins from t0 to t1
//
//
//
//
//
- for (Int_t it=0;it<=t1-t0; it++){
+ for (Int_t it=0;it<100; it++){
x[it]=0;
yt[it]=0;
zt[it]=0;
- clusters[it+t0]=-2;
+ clusters[it]=-2;
zmean[it]=0;
nmean[it]=0;
//
Int_t nfound=0;
Double_t h01 =0;
Int_t plane =-1;
+ Int_t detector =-1;
Float_t padlength=0;
AliTRDtrack track2(*track);
Float_t snpy = track->GetSnp();
AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
h01 = GetTiltFactor(c);
if (plane<0){
- Int_t det = c->GetDetector();
+ Int_t det = c->GetDetector();
plane = fGeom->GetPlane(det);
padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
}
//
clfound++;
if (chi2 > maxChi2[1]) continue;
+ detector = c->GetDetector();
for (Int_t ih=2;ih<9; ih++){ //store the clusters in the road
if (cl[ih][it]==0){
//
Double_t sigmas[10];
Double_t tchi2s[10]; // chi2s for tracklet
+
+ for (Int_t it=0;it<10;it++) {
+
+ ngood[it] = 0;
+ nbad[it] = 0;
+ //
+ meanz[it] = 0;
+ moffset[it] = 0; // mean offset
+ mean[it] = 0; // mean value
+ angle[it] = 0; // angle
+ //
+ smoffset[it] = 1e10; // sigma of mean offset
+ smean[it] = 1e10; // sigma of mean value
+ sangle[it] = 1e10; // sigma of angle
+ smeanangle[it] = 0; // correlation
+ //
+ sigmas[it] = 1e10;
+ tchi2s[it] = 1e10; // chi2s for tracklet
+
+ }
+
//
// calculate zmean
//
TGraph graphz(t1-t0,x,zt);
//
//
+ if (AliTRDReconstructor::StreamLevel()>0)
cstream<<"tracklet"<<
"track.="<<track<< // track parameters
"tany="<<tany<< // tangent of the local track angle
"clfound="<<clfound<< // total number of found clusters in road
"mpads="<<mpads<< // mean number of pads per cluster
"plane="<<plane<< // plane number
+ "detector="<<detector<< // detector number
"road="<<road<< // the width of the used road
"graph0.="<<&graph0<< // x - y = dy for closest cluster
"graph1.="<<&graph1<< // x - y = dy for second closest cluster
//
//
return nfound;
-}
+}
-Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist, Int_t *outlist, Bool_t down)
+//_____________________________________________________________________________
+Int_t AliTRDtracker::Freq(Int_t n, const Int_t *inlist
+ , Int_t *outlist, Bool_t down)
{
//
// Sort eleements according occurancy
// The size of output array has is 2*n
//
+
Int_t * sindexS = new Int_t[n]; // temp array for sorting
Int_t * sindexF = new Int_t[2*n];
for (Int_t i=0;i<n;i++) sindexF[i]=0;
delete [] sindexF;
return countPos;
+
}
-AliTRDtrack * AliTRDtracker::RegisterSeed(AliTRDseed * seeds, Double_t * params)
+//_____________________________________________________________________________
+AliTRDtrack *AliTRDtracker::RegisterSeed(AliTRDseed * seeds, Double_t * params)
{
//
- //
+ // Register a seed
//
Double_t alpha=AliTRDgeometry::GetAlpha();
Double_t shift=AliTRDgeometry::GetAlpha()/2.;
Int_t index =0;
AliTRDcluster *cl =0;
for (Int_t ilayer=0;ilayer<6;ilayer++){
- if (seeds[ilayer].isOK()){
+ if (seeds[ilayer].IsOK()){
for (Int_t itime=22;itime>0;itime--){
if (seeds[ilayer].fIndexes[itime]>0){
index = seeds[ilayer].fIndexes[itime];
CookdEdxTimBin(*track);
CookLabel(track, 0.9);
}
- return track;
-}
-
-
-
-
-
-
-AliTRDseed::AliTRDseed()
-{
- //
- //
- fTilt =0; // tilting angle
- fPadLength = 0; // pad length
- fX0 = 0; // x0 position
- for (Int_t i=0;i<25;i++){
- fX[i]=0; // !x position
- fY[i]=0; // !y position
- fZ[i]=0; // !z position
- fIndexes[i]=0; // !indexes
- fClusters[i]=0; // !clusters
- }
- for (Int_t i=0;i<2;i++){
- fYref[i]=0; // reference y
- fZref[i]=0; // reference z
- fYfit[i]=0; // y fit position +derivation
- fYfitR[i]=0; // y fit position +derivation
- fZfit[i]=0; // z fit position
- fZfitR[i]=0; // z fit position
- fLabels[i]=0; // labels
- }
- fSigmaY = 0;
- fSigmaY2 = 0;
- fMeanz=0; // mean vaue of z
- fZProb=0; // max probbable z
- fMPads=0;
- //
- fN=0; // number of associated clusters
- fN2=0; // number of not crossed
- fNUsed=0; // number of used clusters
- fNChange=0; // change z counter
-}
-
-void AliTRDseed::Reset(){
- //
- // reset seed
- //
- for (Int_t i=0;i<25;i++){
- fX[i]=0; // !x position
- fY[i]=0; // !y position
- fZ[i]=0; // !z position
- fIndexes[i]=0; // !indexes
- fClusters[i]=0; // !clusters
- fUsable[i] = kFALSE;
- }
- for (Int_t i=0;i<2;i++){
- fYref[i]=0; // reference y
- fZref[i]=0; // reference z
- fYfit[i]=0; // y fit position +derivation
- fYfitR[i]=0; // y fit position +derivation
- fZfit[i]=0; // z fit position
- fZfitR[i]=0; // z fit position
- fLabels[i]=-1; // labels
- }
- fSigmaY =0; //"robust" sigma in y
- fSigmaY2=0; //"robust" sigma in y
- fMeanz =0; // mean vaue of z
- fZProb =0; // max probbable z
- fMPads =0;
- //
- fN=0; // number of associated clusters
- fN2=0; // number of not crossed
- fNUsed=0; // number of used clusters
- fNChange=0; // change z counter
-}
-
-void AliTRDseed::CookLabels(){
- //
- // cook 2 labels for seed
- //
- Int_t labels[200];
- Int_t out[200];
- Int_t nlab =0;
- for (Int_t i=0;i<25;i++){
- if (!fClusters[i]) continue;
- for (Int_t ilab=0;ilab<3;ilab++){
- if (fClusters[i]->GetLabel(ilab)>=0){
- labels[nlab] = fClusters[i]->GetLabel(ilab);
- nlab++;
- }
- }
- }
- Int_t nlab2 = AliTRDtracker::Freq(nlab,labels,out,kTRUE);
- fLabels[0] = out[0];
- if (nlab2>1 && out[3]>1) fLabels[1] =out[2];
-}
-
-void AliTRDseed::UseClusters()
-{
- //
- // use clusters
- //
- for (Int_t i=0;i<25;i++){
- if (!fClusters[i]) continue;
- if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
- }
-}
-
-
-void AliTRDseed::Update(){
- //
- //
- //
- const Float_t ratio = 0.8;
- const Int_t kClmin = 6;
- const Float_t kmaxtan = 2;
- if (TMath::Abs(fYref[1])>kmaxtan) return; // too much inclined track
- //
- 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
- fNChange =0;
- //
- Double_t sumw, sumwx,sumwx2;
- Double_t sumwy, sumwxy, sumwz,sumwxz;
- Int_t zints[25]; // histograming of the z coordinate - get 1 and second max probable coodinates in z
- Int_t zouts[50]; //
- Float_t allowedz[25]; // allowed z for given time bin
- Float_t yres[25]; // residuals from reference
- Float_t anglecor = fTilt*fZref[1]; //correction to the angle
- //
- //
- fN=0; fN2 =0;
- for (Int_t i=0;i<25;i++){
- yres[i] =10000;
- if (!fClusters[i]) continue;
- yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i]; // residual y
- zints[fN] = Int_t(fZ[i]);
- fN++;
- }
- if (fN<kClmin) return;
- Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
- fZProb = zouts[0];
- if (nz<=1) zouts[3]=0;
- if (zouts[1]+zouts[3]<kClmin) return;
- //
- if (TMath::Abs(zouts[0]-zouts[2])>12.) zouts[3]=0; // z distance bigger than pad - length
- //
- Int_t breaktime = -1;
- Bool_t mbefore = kFALSE;
- Int_t cumul[25][2];
- Int_t counts[2]={0,0};
- //
- if (zouts[3]>=3){
- //
- // find the break time allowing one chage on pad-rows with maximal numebr of accepted clusters
- //
- fNChange=1;
- for (Int_t i=0;i<25;i++){
- cumul[i][0] = counts[0];
- cumul[i][1] = counts[1];
- if (TMath::Abs(fZ[i]-zouts[0])<2) counts[0]++;
- if (TMath::Abs(fZ[i]-zouts[2])<2) counts[1]++;
- }
- Int_t maxcount = 0;
- for (Int_t i=0;i<24;i++) {
- Int_t after = cumul[24][0]-cumul[i][0];
- Int_t before = cumul[i][1];
- if (after+before>maxcount) {
- maxcount=after+before;
- breaktime=i;
- mbefore=kFALSE;
- }
- after = cumul[24][1]-cumul[i][1];
- before = cumul[i][0];
- if (after+before>maxcount) {
- maxcount=after+before;
- breaktime=i;
- mbefore=kTRUE;
- }
- }
- breaktime-=1;
- }
- for (Int_t i=0;i<25;i++){
- if (i>breaktime) allowedz[i] = mbefore ? zouts[2]:zouts[0];
- if (i<=breaktime) allowedz[i] = (!mbefore) ? zouts[2]:zouts[0];
- }
- if ( (allowedz[0]>allowedz[24] && fZref[1]<0) || (allowedz[0]<allowedz[24] && fZref[1]>0)){
- //
- // tracklet z-direction not in correspondance with track z direction
- //
- fNChange =0;
- for (Int_t i=0;i<25;i++){
- allowedz[i] = zouts[0]; //only longest taken
- }
- }
- //
- if (fNChange>0){
- //
- // cross pad -row tracklet - take the step change into account
- //
- for (Int_t i=0;i<25;i++){
- if (!fClusters[i]) continue;
- if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
- yres[i] = fY[i]-fYref[0]-(fYref[1]+anglecor)*fX[i]; // residual y
- if (TMath::Abs(fZ[i]-fZProb)>2){
- if (fZ[i]>fZProb) yres[i]+=fTilt*fPadLength;
- if (fZ[i]<fZProb) yres[i]-=fTilt*fPadLength;
- }
- }
- }
- //
- Double_t yres2[25];
- Double_t mean,sigma;
- for (Int_t i=0;i<25;i++){
- if (!fClusters[i]) continue;
- if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
- yres2[fN2] = yres[i];
- fN2++;
- }
- if (fN2<kClmin){
- fN2 = 0;
- return;
- }
- EvaluateUni(fN2,yres2,mean,sigma,Int_t(fN2*ratio-2));
- if (sigma<sigmaexp*0.8) sigma=sigmaexp;
- fSigmaY = sigma;
- //
- //
- // reset sums
- sumw=0; sumwx=0; sumwx2=0;
- sumwy=0; sumwxy=0; sumwz=0;sumwxz=0;
- fN2 =0;
- fMeanz =0;
- fMPads =0;
- //
- for (Int_t i=0;i<25;i++){
- fUsable[i]=kFALSE;
- if (!fClusters[i]) continue;
- if (TMath::Abs(fZ[i]-allowedz[i])>2) continue;
- if (TMath::Abs(yres[i]-mean)>4.*sigma) continue;
- fUsable[i] = kTRUE;
- fN2++;
- fMPads+=fClusters[i]->GetNPads();
- Float_t weight =1;
- if (fClusters[i]->GetNPads()>4) weight=0.5;
- if (fClusters[i]->GetNPads()>5) weight=0.2;
- //
- Double_t x = fX[i];
- sumw+=weight; sumwx+=x*weight; sumwx2+=x*x*weight;
- sumwy+=weight*yres[i]; sumwxy+=weight*(yres[i])*x;
- sumwz+=weight*fZ[i]; sumwxz+=weight*fZ[i]*x;
- }
- if (fN2<kClmin){
- fN2 = 0;
- return;
- }
- fMeanz = sumwz/sumw;
- Float_t correction =0;
- if (fNChange>0){
- // tracklet on boundary
- if (fMeanz<fZProb) correction = ycrosscor;
- if (fMeanz>fZProb) correction = -ycrosscor;
- }
- 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<25;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];
- //
- //
- UpdateUsed();
-}
-
-
-
-
-
-
-void AliTRDseed::UpdateUsed(){
- //
- fNUsed =0;
- for (Int_t i=0;i<25;i++){
- if (!fClusters[i]) continue;
- if ((fClusters[i]->IsUsed())) fNUsed++;
- }
-}
-
-
-void AliTRDseed::EvaluateUni(Int_t nvectors, Double_t *data, Double_t &mean, Double_t &sigma, Int_t hh)
-{
- //
- // robust estimator in 1D case MI version
- //
- //for the univariate case
- //estimates of location and scatter are returned in mean and sigma parameters
- //the algorithm works on the same principle as in multivariate case -
- //it finds a subset of size hh with smallest sigma, and then returns mean and
- //sigma of this subset
-
- if (hh==0)
- hh=(nvectors+2)/2;
- Double_t faclts[]={2.6477,2.5092,2.3826,2.2662,2.1587,2.0589,1.9660,1.879,1.7973,1.7203,1.6473};
- Int_t *index=new Int_t[nvectors];
- TMath::Sort(nvectors, data, index, kFALSE);
- //
- Int_t nquant = TMath::Min(Int_t(Double_t(((hh*1./nvectors)-0.5)*40))+1, 11);
- Double_t factor = faclts[nquant-1];
- //
- //
- Double_t sumx =0;
- Double_t sumx2 =0;
- Int_t bestindex = -1;
- Double_t bestmean = 0;
- Double_t bestsigma = data[index[nvectors-1]]-data[index[0]]; // maximal possible sigma
- for (Int_t i=0; i<hh; i++){
- sumx += data[index[i]];
- sumx2 += data[index[i]]*data[index[i]];
- }
- //
- Double_t norm = 1./Double_t(hh);
- Double_t norm2 = 1./Double_t(hh-1);
- for (Int_t i=hh; i<nvectors; i++){
- Double_t cmean = sumx*norm;
- Double_t csigma = (sumx2 - hh*cmean*cmean)*norm2;
- if (csigma<bestsigma){
- bestmean = cmean;
- bestsigma = csigma;
- bestindex = i-hh;
- }
- //
- //
- sumx += data[index[i]]-data[index[i-hh]];
- sumx2 += data[index[i]]*data[index[i]]-data[index[i-hh]]*data[index[i-hh]];
- }
-
- Double_t bstd=factor*TMath::Sqrt(TMath::Abs(bestsigma));
- mean = bestmean;
- sigma = bstd;
- delete [] index;
-}
+ return track;
-Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror){
- //
- //
- //
- TLinearFitter fitterT2(4,"hyp4"); // fitting with tilting pads - kz not fixed
- fitterT2.StoreData(kTRUE);
- Float_t xref2 = (cseed[2].fX0+cseed[3].fX0)*0.5; // reference x0 for z
- //
- Int_t npointsT =0;
- fitterT2.ClearPoints();
- for (Int_t iLayer=0; iLayer<6;iLayer++){
- if (!cseed[iLayer].isOK()) continue;
- Double_t tilt = cseed[iLayer].fTilt;
-
- for (Int_t itime=0;itime<25;itime++){
- if (!cseed[iLayer].fUsable[itime]) continue;
- Double_t x = cseed[iLayer].fX[itime]+cseed[iLayer].fX0-xref2; // x relative to the midle chamber
- Double_t y = cseed[iLayer].fY[itime];
- Double_t z = cseed[iLayer].fZ[itime];
- // tilted rieman
- //
- Double_t uvt[6];
- Double_t x2 = cseed[iLayer].fX[itime]+cseed[iLayer].fX0; // global x
- Double_t t = 1./(x2*x2+y*y);
- uvt[1] = t; // t
- uvt[0] = 2.*x2*uvt[1]; // u
- uvt[2] = 2.0*tilt*uvt[1];
- uvt[3] = 2.0*tilt*x*uvt[1];
- uvt[4] = 2.0*(y+tilt*z)*uvt[1];
- //
- Double_t error = 2*uvt[1];
- if (terror) error*=cseed[iLayer].fSigmaY;
- else {error *=0.2;} //default error
- fitterT2.AddPoint(uvt,uvt[4],error);
- npointsT++;
- }
- }
- fitterT2.Eval();
- Double_t rpolz0 = fitterT2.GetParameter(3);
- Double_t rpolz1 = fitterT2.GetParameter(4);
- //
- // linear fitter - not possible to make boundaries
- // non accept non possible z and dzdx combination
- //
- Bool_t acceptablez =kTRUE;
- for (Int_t iLayer=0; iLayer<6;iLayer++){
- if (cseed[iLayer].isOK()){
- Double_t zT2 = rpolz0+rpolz1*(cseed[iLayer].fX0 - xref2);
- if (TMath::Abs(cseed[iLayer].fZProb-zT2)>cseed[iLayer].fPadLength*0.5+1)
- acceptablez = kFALSE;
- }
- }
- if (!acceptablez){
- Double_t zmf = cseed[2].fZref[0]+cseed[2].fZref[1]*(xref2-cseed[2].fX0);
- Double_t dzmf = (cseed[2].fZref[1]+ cseed[3].fZref[1])*0.5;
- fitterT2.FixParameter(3,zmf);
- fitterT2.FixParameter(4,dzmf);
- fitterT2.Eval();
- fitterT2.ReleaseParameter(3);
- fitterT2.ReleaseParameter(4);
- rpolz0 = fitterT2.GetParameter(3);
- rpolz1 = fitterT2.GetParameter(4);
- }
- //
- Double_t chi2TR = fitterT2.GetChisquare()/Float_t(npointsT);
- Double_t params[3];
- params[0] = fitterT2.GetParameter(0);
- params[1] = fitterT2.GetParameter(1);
- params[2] = fitterT2.GetParameter(2);
- Double_t CR = 1+params[1]*params[1]-params[2]*params[0];
- for (Int_t iLayer = 0; iLayer<6;iLayer++){
- Double_t x = cseed[iLayer].fX0;
- Double_t y=0,dy=0, z=0, dz=0;
- // y
- Double_t res2 = (x*params[0]+params[1]);
- res2*=res2;
- res2 = 1.-params[2]*params[0]+params[1]*params[1]-res2;
- if (res2>=0){
- res2 = TMath::Sqrt(res2);
- y = (1-res2)/params[0];
- }
- //dy
- Double_t x0 = -params[1]/params[0];
- if (-params[2]*params[0]+params[1]*params[1]+1>0){
- Double_t Rm1 = params[0]/TMath::Sqrt(-params[2]*params[0]+params[1]*params[1]+1);
- if ( 1./(Rm1*Rm1)-(x-x0)*(x-x0)>0){
- Double_t res = (x-x0)/TMath::Sqrt(1./(Rm1*Rm1)-(x-x0)*(x-x0));
- if (params[0]<0) res*=-1.;
- dy = res;
- }
- }
- z = rpolz0+rpolz1*(x-xref2);
- dz = rpolz1;
- cseed[iLayer].fYref[0] = y;
- cseed[iLayer].fYref[1] = dy;
- cseed[iLayer].fZref[0] = z;
- cseed[iLayer].fZref[1] = dz;
- cseed[iLayer].fC = CR;
- //
- }
- return chi2TR;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff