tracklets[ip].SetChi2(chi2);
}
}
-
//update track points array
if(np && points){
Float_t xyz[3];
for(int ip=0; ip<np; ip++){
points[ip].GetXYZ(xyz);
- xyz[1] = y0 - (y0>0.?1.:-1.)*TMath::Sqrt(R*R-(xyz[0]-x0)*(xyz[0]-x0));
+ xyz[1] = TMath::Abs(xyz[0] - x0) > R ? 100. : y0 - (y0>0.?1.:-1.)*TMath::Sqrt(R*R-(xyz[0]-x0)*(xyz[0]-x0));
xyz[2] = z0 + dzdx * (xyz[0] - xref);
points[ip].SetXYZ(xyz);
}
//
// Author : A.Bercuci@gsi.de
- //printf("Start track @ x[%f]\n", track->GetX());
+ // printf("Start track @ x[%f]\n", track->GetX());
//prepare marker points along the track
Int_t ip = np ? 0 : 1;
while(ip < np){
//don't do anything if next marker is after next update point.
if((up?-1:1) * (points[ip].GetX() - x) - fgkMaxStep < 0) break;
-
- //printf("Propagate to x[%d] = %f\n", ip, points[ip].GetX());
-
if(((up?-1:1) * (points[ip].GetX() - track->GetX()) < 0) && !PropagateToX(*track, points[ip].GetX(), fgkMaxStep)) return -1.;
Double_t xyz[3]; // should also get the covariance
- track->GetXYZ(xyz); points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]);
+ track->GetXYZ(xyz);
+ track->Global2LocalPosition(xyz, track->GetAlpha());
+ points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]);
ip++;
}
- //printf("plane[%d] tracklet[%p] x[%f]\n", iplane, ptrTracklet, x);
+ // printf("plane[%d] tracklet[%p] x[%f]\n", iplane, ptrTracklet, x);
- //Propagate closer to the next update point
+ // Propagate closer to the next update point
if(((up?-1:1) * (x - track->GetX()) + fgkMaxStep < 0) && !PropagateToX(*track, x + (up?-1:1)*fgkMaxStep, fgkMaxStep)) return -1.;
if(!AdjustSector(track)) return -1;
xyz1[0] = x * TMath::Cos(alpha) - y * TMath::Sin(alpha);
xyz1[1] = +x * TMath::Sin(alpha) + y * TMath::Cos(alpha);
xyz1[2] = z;
+ if((xyz0[0] - xyz1[9] < 1e-3) && (xyz0[0] - xyz1[9] < 1e-3)) continue; // check wheter we are at the same global x position
Double_t param[7];
- if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param)<=0.) break;
+ if(AliTracker::MeanMaterialBudget(xyz0, xyz1, param) <=0.) break;
Double_t xrho = param[0]*param[4]; // density*length
Double_t xx0 = param[1]; // radiation length
//Update track
Double_t chi2 = track->GetPredictedChi2(ptrTracklet);
if(chi2<1e+10) track->Update(ptrTracklet, chi2);
-
if(!up) continue;
//Reset material budget if 2 consecutive gold
if(((up?-1:1) * (points[ip].GetX() - track->GetX()) < 0) && !PropagateToX(*track, points[ip].GetX(), fgkMaxStep)) return -1.;
Double_t xyz[3]; // should also get the covariance
- track->GetXYZ(xyz); points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]);
+ track->GetXYZ(xyz);
+ track->Global2LocalPosition(xyz, track->GetAlpha());
+ points[ip].SetXYZ(xyz[0], xyz[1], xyz[2]);
ip++;
}
Double_t alpha = AliTRDgeometry::GetAlpha();
Double_t y = track->GetY();
Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
-
+
if (y > ymax) {
if (!track->Rotate( alpha)) {
return kFALSE;
#include "AliTRDcluster.h"
#include "AliESDHeader.h"
#include "AliESDRun.h"
+#include "AliESDtrack.h"
#include "AliTRDgeometry.h"
+#include "AliTRDpadPlane.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
+#include "AliTRDtrackerV1.h"
+#include "AliTRDReconstructor.h"
#include "AliTrackReference.h"
+#include "AliTrackPointArray.h"
+#include "AliTracker.h"
#include "TTreeStream.h"
#include "AliTRDtrackInfo/AliTRDtrackInfo.h"
AliTRDrecoTask("DetChecker", "Basic Detector Checker")
,fEventInfo(0x0)
,fTriggerNames(0x0)
+ ,fReconstructor(0x0)
{
//
// Default constructor
//
DefineInput(1,AliTRDeventInfo::Class());
+ fReconstructor = new AliTRDReconstructor;
+ fReconstructor->SetRecoParam(AliTRDrecoParam::GetLowFluxParam());
InitFunctorList();
}
//
if(fEventInfo) delete fEventInfo;
if(fTriggerNames) delete fTriggerNames;
+ delete fReconstructor;
}
//_______________________________________________________
fContainer->AddAt(new TH1F("hEventsTriggerTracks", "Trigger Class (Tracks)", 100, 0, 100), kNEventsTriggerTracks);
fContainer->AddAt(new TH1F("hNcls", "Nr. of clusters per track", 181, -0.5, 180.5), kNclustersHist);
fContainer->AddAt(new TH1F("hNtls", "Nr. tracklets per track", 7, -0.5, 6.5), kNtrackletsHist);
+ fContainer->AddAt(new TH1F("hNtlsFindable", "Ratio of found/findable Tracklets" , 11, -0.05, 1.05), kNTrackletsVsFindable);
fContainer->AddAt(new TH1F("hNclTls","Mean Number of clusters per tracklet", 31, -0.5, 30.5), kNclusterTrackletHist);
fContainer->AddAt(new TH1F("hChi2", "Chi2", 200, 0, 20), kChi2);
fContainer->AddAt(new TH1F("hChi2n", "Norm. Chi2 (tracklets)", 50, 0, 5), kChi2Normalized);
if(fDebugLevel > 2){
Int_t crossing = tracklet->GetNChange();
AliTRDcluster *c = 0x0;
- for(Int_t itime = 0; itime < kNTimeBins; itime++){
+ for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
if(!(c = tracklet->GetClusters(itime))) continue;
break;
}
AliWarning("No Histogram defined.");
return 0x0;
}
- Int_t nTracklets = 0;
+ Int_t nTracklets = GetNTracklets(track);
+ h->Fill(nTracklets);
+ if(fDebugLevel > 3){
+ if(nTracklets == 1){
+ // If we have one Tracklet, check in which layer this happens
+ Int_t layer = -1;
+ AliTRDseedV1 *tracklet = 0x0;
+ for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+ if((tracklet = track->GetTracklet(il)) && tracklet->IsOK()){layer = il; break;}
+ }
+ (*fDebugStream) << "PlotNTracklets"
+ << "Layer=" << layer
+ << "\n";
+ }
+ }
+ return h;
+}
+
+//_______________________________________________________
+TH1 *AliTRDcheckDetector::PlotTrackletsVsFindable(const AliTRDtrackV1 *track){
+ //
+ // Plots the ratio of number of tracklets vs.
+ // number of findable tracklets
+ //
+ // Findable tracklets are defined as track prolongation
+ // to layer i does not hit the dead area +- epsilon
+ //
+ const Float_t epsilon = 0.01; // dead area tolerance
+ const Float_t epsilon_R = 1;
+ Double_t x_anode[AliTRDgeometry::kNlayer] = {300.2, 312.8, 325.4, 338.0, 350.6, 363.2}; // Take the default X0
+
+ if(track) fTrack = track;
+ if(!fTrack){
+ AliWarning("No Track defined.");
+ return 0x0;
+ }
+ TH1 *h = 0x0;
+ if(!(h = dynamic_cast<TH1F *>(fContainer->At(kNTrackletsVsFindable)))){
+ AliWarning("No Histogram defined.");
+ return 0x0;
+ }
+ Int_t nFound = 0, nFindable = 0;
+ Int_t stack = -1;
+ Double_t ymin = 0., ymax = 0., zmin = 0., zmax = 0.;
+ Double_t y = 0., z = 0.;
AliTRDseedV1 *tracklet = 0x0;
- for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
- if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
- nTracklets++;
+ AliTRDpadPlane *pp;
+ AliTRDgeometry *fGeo = new AliTRDgeometry;
+ for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+ if((tracklet = fTrack->GetTracklet(il)) && tracklet->IsOK()){
+ tracklet->SetReconstructor(fReconstructor);
+ nFound++;
+ // printout tracklet
+ }
}
- h->Fill(nTracklets);
+ // 2 Different cases:
+ // 1st stand alone: here we cannot propagate, but be can do a Tilted Rieman Fit
+ // 2nd barrel track: here we propagate the track to the layers
+ AliTrackPoint points[6];
+ Float_t xyz[3];
+ memset(xyz, 0, sizeof(Float_t) * 3);
+ if(((fESD->GetStatus() & AliESDtrack::kTRDout) > 0) && !((fESD->GetStatus() & AliESDtrack::kTRDin) > 0)){
+ // stand alone track
+ for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+ xyz[0] = x_anode[il];
+ points[il].SetXYZ(xyz);
+ }
+ AliTRDtrackerV1::FitRiemanTilt(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kTRUE, 6, points);
+ } else {
+ // barrel track
+ for(Int_t il = AliTRDgeometry::kNlayer; il--;){
+ // The track points have to be in reverse order for the Kalman Filter (since we move back)
+ xyz[0] = x_anode[AliTRDgeometry::kNlayer - il - 1] - epsilon_R;
+ points[il].SetXYZ(xyz);
+ }
+ for(Int_t ipt = 0; ipt < AliTRDgeometry::kNlayer; ipt++)
+ printf("%d. X = %f\n", ipt, points[ipt].GetX());
+ AliTRDtrackerV1::FitKalman(const_cast<AliTRDtrackV1 *>(fTrack), 0x0, kFALSE, 6, points);
+ AliTrackPoint tempPointCont[AliTRDgeometry::kNlayer];
+ memcpy(tempPointCont, points, sizeof(AliTrackPoint) * AliTRDgeometry::kNlayer);
+ for(Int_t il = AliTRDgeometry::kNlayer; il--; ){
+ tempPointCont[il].GetXYZ(xyz);
+ points[AliTRDgeometry::kNlayer - il -1].SetXYZ(xyz);
+ }
+ }
+ for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+ y = points[il].GetY();
+ z = points[il].GetZ();
+ if((stack = fGeo->GetStack(z, il)) < 0) continue; // Not findable
+ pp = fGeo->GetPadPlane(il, stack);
+ ymin = pp->GetCol0() + epsilon;
+ ymax = pp->GetColEnd() - epsilon;
+ zmin = pp->GetRowEnd() + epsilon;
+ zmax = pp->GetRow0() - epsilon;
+ // ignore y-crossing (material)
+ if((z > zmin && z < zmax) && (y > ymin && y < ymax)) nFindable++;
+/*
+ if(fDebugLevel > 3){
+ Double_t pos_tracklet[2] = {tracklet ? tracklet->GetYfit(0) : 0, tracklet ? tracklet->GetMeanz() : 0};
+ Int_t hasTracklet = tracklet ? 1 : 0;
+ (*fDebugStream) << "GetFindableTracklets"
+ << "layer=" << il
+ << "ytracklet=" << pos_tracklet[0]
+ << "ytrack=" << y
+ << "ztracklet=" << pos_tracklet[1]
+ << "ztrack=" << z
+ << "tracklet=" << hasTracklet
+ << "\n";
+ }
+*/
+ }
+
+ delete fGeo;
+ h->Fill(nFindable > 0 ? TMath::Min(nFound/static_cast<Double_t>(nFindable), 1.) : 1);
+ if(fDebugLevel > 2) AliInfo(Form("Findable[Found]: %d[%d|%f]", nFindable, nFound, nFound/static_cast<Float_t>(nFindable > 0 ? nFindable : 1)));
return h;
}
AliTRDcluster *c = 0x0;
for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
- for(Int_t itime = 0; itime < kNTimeBins; itime++){
+ for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
if(!(c = tracklet->GetClusters(itime))) continue;
Int_t localtime = c->GetLocalTimeBin();
Double_t absolute_charge = TMath::Abs(c->GetQ());
AliTRDcluster *c = 0x0;
for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK())continue;
- for(Int_t itime = 0; itime < kNTimeBins; itime++){
+ for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
if(!(c = tracklet->GetClusters(itime))) continue;
h->Fill(c->GetQ());
}
AliTRDseedV1 *tracklet = 0x0;
AliTRDcluster *c = 0x0, *c1 = 0x0; // c1 for the Debug Stream
Double_t Qtot = 0;
+ Int_t nTracklets = 0;
+ if(fDebugLevel > 3)
+ nTracklets = GetNTracklets(fTrack); // fill NTracklet to the Debug Stream
for(Int_t itl = 0x0; itl < AliTRDgeometry::kNlayer; itl++){
if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
Qtot = 0;
c1 = 0x0;
- for(Int_t itime = 0; itime < kNTimeBins; itime++){
+ for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
if(!(c = tracklet->GetClusters(itime))) continue;
if(!c1) c1 = c;
Qtot += TMath::Abs(c->GetQ());
<< "Sector=" << sector
<< "crossing=" << crossing
<< "momentum=" << momentum
+ << "nTracklets="<< nTracklets
<< "pdg=" << pdg
<< "theta=" << theta
<< "phi=" << phi
Int_t sector = -1;
for(Int_t itl = 0; itl < AliTRDgeometry::kNlayer; itl++){
if(!(tracklet = fTrack->GetTracklet(itl)) || !tracklet->IsOK()) continue;
- for(Int_t itime = 0; itime < kNTimeBins; itime++){
+ for(Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins(); itime++){
if(!(c = tracklet->GetClusters(itime))) continue;
sector = static_cast<Int_t>(c->GetDetector()/AliTRDgeometry::kNdets);
}
return h;
}
+//_______________________________________________________
+Int_t AliTRDcheckDetector::GetNTracklets(const AliTRDtrackV1 *track){
+ //
+ // Count the number of tracklets per track
+ //
+ if(!track) return 0;
+ Int_t nTracklets = 0;
+ AliTRDseedV1 *tracklet = 0x0;
+ for(Int_t il = 0; il < AliTRDgeometry::kNlayer; il++){
+ if((tracklet = track->GetTracklet(il)) && tracklet->IsOK()) nTracklets++;
+ }
+ return nTracklets;
+}
+
+//________________________________________________________
+void AliTRDcheckDetector::SetRecoParam(AliTRDrecoParam *r)
+{
+
+ fReconstructor->SetRecoParam(r);
+}
class TH1;
class TMap;
class AliESDHeader;
+class AliTRDReconstructor;
+class AliTRDrecoParam;
+
class AliTRDcheckDetector : public AliTRDrecoTask{
-// common constants
-enum{
- kNTimeBins = 30
-};
-// The Histogram number
-enum{
- kNTracksEventHist=0,
- kNEventsTriggerTracks=1,
- kNclustersHist=2,
- kNtrackletsHist=3,
- kNclusterTrackletHist=4,
- kChi2=5,
- kChi2Normalized=6,
- kNTracksSectorHist=7,
- kPulseHeight=8,
- kClusterCharge=9,
- kChargeDeposit=10,
- kNEventsTrigger=11,
- kPurity = 12
-};
-public:
- AliTRDcheckDetector();
- virtual ~AliTRDcheckDetector();
+ // The Histogram number
+ typedef enum{
+ kNTracksEventHist=0,
+ kNEventsTriggerTracks=1,
+ kNclustersHist=2,
+ kNtrackletsHist=3,
+ kNTrackletsVsFindable = 4,
+ kNclusterTrackletHist=5,
+ kChi2=6,
+ kChi2Normalized=7,
+ kNTracksSectorHist=8,
+ kPulseHeight=9,
+ kClusterCharge=10,
+ kChargeDeposit=11,
+ kNEventsTrigger=12,
+ kPurity = 12
+ }HistType_t;
+ public:
+ AliTRDcheckDetector();
+ virtual ~AliTRDcheckDetector();
- virtual void ConnectInputData(const Option_t *);
- virtual void CreateOutputObjects();
- virtual void Exec(Option_t *);
- virtual void Terminate(Option_t *);
+ virtual void ConnectInputData(const Option_t *);
+ virtual void CreateOutputObjects();
+ virtual void Exec(Option_t *);
+ virtual void Terminate(Option_t *);
- virtual TObjArray *Histos();
+ virtual TObjArray *Histos();
- // Plotting Functions:
- TH1 *PlotMeanNClusters(const AliTRDtrackV1 *t = 0x0);
- TH1 *PlotNClusters(const AliTRDtrackV1 *t = 0x0);
- TH1 *PlotNTracklets(const AliTRDtrackV1 *t = 0x0);
- TH1 *PlotTracksSector(const AliTRDtrackV1 *t = 0x0);
- TH1 *PlotPulseHeight(const AliTRDtrackV1 *t = 0x0);
- TH1 *PlotChi2(const AliTRDtrackV1 *t = 0x0);
- TH1 *PlotNormalizedChi2(const AliTRDtrackV1 *t = 0x0);
- TH1 *PlotClusterCharge(const AliTRDtrackV1 *t = 0x0);
- TH1 *PlotChargeDeposit(const AliTRDtrackV1 *t = 0x0);
-
- virtual Bool_t PostProcess();
- virtual void GetRefFigure(Int_t ifig);
+ // Plotting Functions:
+ TH1 *PlotMeanNClusters(const AliTRDtrackV1 *t = 0x0);
+ TH1 *PlotNClusters(const AliTRDtrackV1 *t = 0x0);
+ TH1 *PlotNTracklets(const AliTRDtrackV1 *t = 0x0);
+ TH1 *PlotTrackletsVsFindable(const AliTRDtrackV1 *track = 0x0);
+ TH1 *PlotTracksSector(const AliTRDtrackV1 *t = 0x0);
+ TH1 *PlotPulseHeight(const AliTRDtrackV1 *t = 0x0);
+ TH1 *PlotChi2(const AliTRDtrackV1 *t = 0x0);
+ TH1 *PlotNormalizedChi2(const AliTRDtrackV1 *t = 0x0);
+ TH1 *PlotClusterCharge(const AliTRDtrackV1 *t = 0x0);
+ TH1 *PlotChargeDeposit(const AliTRDtrackV1 *t = 0x0);
+
+ virtual Bool_t PostProcess();
+ virtual void GetRefFigure(Int_t ifig);
+
+ void SetRecoParam(AliTRDrecoParam *r);
-private:
- AliTRDcheckDetector(const AliTRDcheckDetector &);
- AliTRDcheckDetector& operator=(const AliTRDcheckDetector &);
- AliTRDeventInfo *fEventInfo; //! ESD Header
- TMap *fTriggerNames; //! Containing trigger class names
+ private:
+ AliTRDcheckDetector(const AliTRDcheckDetector &);
+ AliTRDcheckDetector& operator=(const AliTRDcheckDetector &);
+ Int_t GetNTracklets(const AliTRDtrackV1 *track);
+ AliTRDeventInfo *fEventInfo; //! ESD Header
+ TMap *fTriggerNames; //! Containing trigger class names
+ AliTRDReconstructor *fReconstructor; //
+
ClassDef(AliTRDcheckDetector, 1)
};
#endif