Double_t fFreq = cp->GetSamplingFrequency();
//drift time corresponding to the center of the time bin
Double_t td = (fPadTime + .5)/fFreq; // [us]
- if(td < t0+2.5) return 0.; // do not calculate radial posion of clusters in the amplification region
-
// correction for t0
td -= t0;
-
- Double_t x = vd*td, xold=0.;
- Float_t tc0 = 0.244, // TRF rising time 0.2us
- dtcdx = 0.009, // diffusion contribution to the rising time of the signal
- kTC = 0.; // tail cancellation residual
- while(TMath::Abs(x-xold)>1.e-3){ // convergence on 10um level
- xold = x;
- Float_t tc = tc0 - dtcdx*x;
- Float_t tq = 0.;
- if(q && xq){
- for(Int_t iq=0; iq<3; iq++) tq += q[iq]*TMath::Exp(-kTC*(x - xq[iq]));
- }
- Float_t vdcorr = x/cp->TimeStruct(vd, x+.5*AliTRDgeometry::CamHght(), z);
- x = (td - tc - tq) * vdcorr;
+ // calculate radial posion of clusters in the drift region
+ if(td < .2 || td > 2.4) return 0.;
+ // correction for TRF rising time 0.2us
+ td -= 0.189;
+
+ // invert drift time function
+ Double_t xM= AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght(),
+ x = vd*td + .5*AliTRDgeometry::CamHght(),
+ t = cp->TimeStruct(vd, x, z), dx1=0.,dx2;
+ while(TMath::Abs(td-t)>1.e-4){ // convergence on 100ps
+ dx2 = vd*(td-t);
+ if(TMath::Abs(TMath::Abs(dx2)-TMath::Abs(dx1))<1.e-6){
+ x+=.5*dx2;
+ break;
+ } else x+=dx2;
+
+ if(x<0. || x>xM) return 0.;
+ t = cp->TimeStruct(vd, x, z);
+ dx1 = dx2;
}
- return x;
+
+ return x-.5*AliTRDgeometry::CamHght();
}
//_____________________________________________________________________________
// drift velocity correction TODO to be moved to the clusterizer
const Float_t cx[] = {
- -9.6280e-02, 1.3091e-01,-1.7415e-02,-9.9221e-02,-1.2040e-01,-9.5493e-02,
- -5.0041e-02,-1.6726e-02, 3.5756e-03, 1.8611e-02, 2.6378e-02, 3.3823e-02,
- 3.4811e-02, 3.5282e-02, 3.5386e-02, 3.6047e-02, 3.5201e-02, 3.4384e-02,
- 3.2864e-02, 3.1932e-02, 3.2051e-02, 2.2539e-02,-2.5154e-02,-1.2050e-01,
- -1.2050e-01
- };
+ 0.0000e+00, 6.0869e-02, -7.0366e-02, -1.4700e-01, -1.6228e-01, -1.3282e-01,
+ -8.7548e-02, -5.3547e-02, -3.2318e-02, -1.7403e-02, -9.6158e-03, -2.7985e-03,
+ -1.1035e-03, -5.1325e-04, 3.9906e-04, 7.6908e-04, 2.5395e-04, -1.7090e-04,
+ -1.8659e-03, -9.8477e-04, -2.2940e-03, -1.3164e-02, -6.6807e-02, -1.5843e-01,
+ 0.0000e+00, };
// PRF correction TODO to be replaced by the gaussian
// approximation with full error parametrization and // moved to the clusterizer
// The input matrix is stored in the vector c and the result in the vector d.
// Both arrays have to be initialized by the user with at least 3 elements. Return negative in case of failure.
//
-// For calculating the square root of the inverse
-// covariance matrix we use the following relation:
+// For calculating the square root of the symmetric matrix c
+// the following relation is used:
// BEGIN_LATEX
-// A^{1/2} = VD^{1/2}V^{-1}
+// C^{1/2} = VD^{1/2}V^{-1}
// END_LATEX
// with V being the matrix with the n eigenvectors as columns.
-// Thus the following are true:
-// - matrix D is diagonal with the diagonal given by the eigenvalues of A (A being symmetric)
+// In case C is symmetric the followings are true:
+// - matrix D is diagonal with the diagonal given by the eigenvalues of C
// - V = V^{-1}
//
// Author A.Bercuci <A.Bercuci@gsi.de>
// uncalibrated cluster correction
// TODO remove
- Double_t x, y; GetClusterXY(c, x, y);
+ Double_t x, y;
+ GetClusterXY(c, x, y);
+ //x = c->GetX(); y = c->GetY();
xc[n] = fX0 - x;
yc[n] = y;
zc[n] = c->GetZ();
Int_t row = i[0];
Int_t col = i[1];
- // Parameters to adjust the X position
- const Double_t kX0shift = 2.52 + 0.04273; //[cm]
- const Double_t kT0shift = 3.19e-3; //[us]
+ // Parameters to adjust the X position of clusters
+ const Double_t kX0shift = AliTRDgeometry::AnodePos(); //[cm]
+ // TRF rising time (fitted)
+ // It should be absorbed by the t0. For the moment t0 is 0 for simulations.
+ // A.Bercuci (Mar 26 2009)
+ const Double_t kT0shift = 0.189; //[us]
if (!fMatrix) {
// T0 correction
Double_t timeT0Cal = time - t0;
// Calculate the X-position,
- Double_t xLocal = ((timeT0Cal + 0.5) / fSamplingFrequency + kT0shift) * vdrift;
+ Double_t xLocal = ((timeT0Cal + 0.5) / fSamplingFrequency - kT0shift) * vdrift;
// Length of the amplification region
Double_t ampLength = (Double_t) AliTRDgeometry::CamHght();
// Invert the X-position,
// apply ExB correction to the Y-position
// and move to the Z-position relative to the middle of the chamber
- posLocal[0] = -xLocal;
+ posLocal[0] = kX0shift-xLocal;
posLocal[1] = (fPadPlane->GetColPos(col) + (0.5 + x[0]) * colSize) - driftLength * exbCorr;
posLocal[2] = (fPadPlane->GetRowPos(row) - 0.5 * rowSize) - fZShiftIdeal;
}
// Output values
- x[0] = posTracking[0] + kX0shift;
+ x[0] = posTracking[0];
x[1] = posTracking[1];
x[2] = posTracking[2];
x[3] = clusterCharge;
ClassImp(AliTRDclusterResolution)
-
+const Float_t AliTRDclusterResolution::fgkTimeBinLength = 1./ AliTRDCommonParam::Instance()->GetSamplingFrequency();
//_______________________________________________________
AliTRDclusterResolution::AliTRDclusterResolution(const char *name)
: AliTRDrecoTask(name, "Cluster Error Parametrization")
,fExB(0.)
,fVdrift(0.)
{
- // ideal equidistant binning
- //fAt = new TAxis(kNTB, -0.075, (kNTB-.5)*.15);
-
- // equidistant binning for scaled for X0-MC (track ref)
- fAt = new TAxis(kNTB, -0.075+.088, (kNTB-.5)*.15+.088);
-
- // non-equidistant binning after vdrift correction
-// const Double_t x[kNTB+1] = {
-// 0.000000, 0.185084, 0.400490, 0.519701, 0.554653, 0.653150, 0.805063, 0.990261,
-// 1.179610, 1.356406, 1.524094, 1.685499, 1.843083, 1.997338, 2.148077, 2.298274,
-// 2.448656, 2.598639, 2.747809, 2.896596, 3.045380, 3.195000, 3.340303, 3.461688,
-// 3.540094, 3.702677};
-// fAt = new TAxis(kNTB, x);
-
+ // time drift axis
+ fAt = new TAxis(kNTB, 0., kNTB*fgkTimeBinLength);
+ // z axis spans the drift cell 2.5 mm
fAd = new TAxis(kND, 0., .25);
// By default register all analysis
// The user can switch them off in his steering macro
- SetProcessCharge();
- SetProcessCenterPad();
- SetProcessMean();
- SetProcessSigma();
+ SetProcess(kQRes);
+ SetProcess(kCenter);
+ SetProcess(kMean);
+ SetProcess(kSigm);
}
//_______________________________________________________
Int_t ih = 0;
for(Int_t id=1; id<=fAd->GetNbins(); id++){
for(Int_t it=1; it<=fAt->GetNbins(); it++){
- arr->AddAt(h2 = new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*fAd->GetBinCenter(id)- 5.*fAd->GetBinWidth(id), 10.*fAd->GetBinCenter(id)+ 5.*fAd->GetBinWidth(id), 10.*fAt->GetBinCenter(it)- 5.*fAt->GetBinWidth(it), 10.*fAt->GetBinCenter(it)+ 5.*fAt->GetBinWidth(it)), 35, -.35, .35, 100, -.5, .5), ih++);
+ arr->AddAt(h2 = new TH2I(Form("hr_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] t_{drift}(%3.1f-%3.1f)[#mus]", 10.*fAd->GetBinLowEdge(id), 10.*fAd->GetBinUpEdge(id), fAt->GetBinLowEdge(it), fAt->GetBinUpEdge(it)), 35, -.35, .35, 100, -.5, .5), ih++);
h2->SetXTitle("tg#phi");
h2->SetYTitle("#Delta y[cm]");
h2->SetZTitle("entries");
ih = 0;
for(Int_t id=1; id<=fAd->GetNbins(); id++){
for(Int_t it=1; it<=fAt->GetNbins(); it++){
- arr->AddAt(h2 = new TH2I(Form("h_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] x_{drift}(%5.2f-%5.2f)[mm]", 10.*fAd->GetBinCenter(id)- 5.*fAd->GetBinWidth(id), 10.*fAd->GetBinCenter(id)+ 5.*fAd->GetBinWidth(id), 10.*fAt->GetBinCenter(it)- 5.*fAt->GetBinWidth(it), 10.*fAt->GetBinCenter(it)+ 5.*fAt->GetBinWidth(it)), 35, -.35, .35, 100, -.5, .5), ih++);
+ arr->AddAt(h2 = new TH2I(Form("hs_d%02dt%02d", id, it), Form("d_{wire}(%3.1f-%3.1f)[mm] t_{drift}(%3.1f-%3.1f)[#mus]", 10.*fAd->GetBinLowEdge(id), 10.*fAd->GetBinUpEdge(id), fAt->GetBinLowEdge(it), fAt->GetBinUpEdge(it)), 35, -.35, .35, 100, -.5, .5), ih++);
h2->SetXTitle("tg#phi-h*tg(#theta)");
h2->SetYTitle("#Delta y[cm]");
h2->SetZTitle("entries");
h2->Fill(cli->GetYDisplacement(), dy);
}
- Int_t it = fAt->FindBin(cli->GetDriftLength());
+ Int_t it = fAt->FindBin((t+.5)*fgkTimeBinLength);
if(it==0 || it == fAt->GetNbins()+1){
- AliWarning(Form("Drift length %f outside allowed range", cli->GetDriftLength()));
+ AliWarning(Form("Drift time %f outside allowed range", t));
continue;
}
Int_t id = fAd->FindBin(cli->GetAnisochronity());
}
arr->AddAt(h2, 0);
h2->SetXTitle("d [mm]");
- h2->SetYTitle("x_{drift} [cm]");
+ h2->SetYTitle("t_{drift} [#mus]");
h2->SetZTitle("#sigma_{x} [cm]");
arr->AddAt(h2 = (TH2F*)h2->Clone("hSY"), 1);
h2->SetZTitle("#sigma_{y} [cm]");
printf("ExB[%e] ExB2[%e]\n", fExB, exb2);
// process resolution dependency on charge
- if(HasProcessCharge()) ProcessCharge();
+ if(HasProcess(kQRes)) ProcessCharge();
// process resolution dependency on y displacement
- if(HasProcessCenterPad()) ProcessCenterPad();
+ if(HasProcess(kCenter)) ProcessCenterPad();
// process resolution dependency on drift legth and drift cell width
- if(HasProcessSigma()) ProcessSigma();
+ if(HasProcess(kSigm)) ProcessSigma();
// process systematic shift on drift legth and drift cell width
- if(HasProcessMean()) ProcessMean();
+ if(HasProcess(kMean)) ProcessMean();
return kTRUE;
}
public:
enum EAxisBinning { // bins in z and x direction
kNTB = 25
- ,kND = 5
- ,kN = kND*kNTB
+ ,kND = 5
+ ,kN = kND*kNTB
};
enum EResultContainer { // results container type
kQRes = 0
- ,kCenter= 1
- ,kSigm = 2
- ,kMean = 3
+ ,kCenter= 1
+ ,kSigm = 2
+ ,kMean = 3
};
enum ECheckBits { // force setting the ExB
- kExB = BIT(23)
- };
- enum ESteeringBits { // steering bits for task
- kSaveAs = 0
- ,kProcCharge = 1
- ,kProcCenterPad = 2
- ,kProcSigma = 3
- ,kProcMean = 4
+ kSaveAs = BIT(22)
+ ,kExB = BIT(23)
};
AliTRDclusterResolution(const char *name="ClErrParam");
virtual ~AliTRDclusterResolution();
- void ConnectInputData(Option_t *);
- void CreateOutputObjects();
- void Exec(Option_t *);
- Int_t GetDetector() const { return fDet; }
+ void ConnectInputData(Option_t *);
+ void CreateOutputObjects();
+ void Exec(Option_t *);
+ Int_t GetDetector() const { return fDet; }
inline Float_t GetExB() const;
inline Float_t GetVdrift() const;
- Bool_t GetRefFigure(Int_t ifig);
- Bool_t HasProcessCharge() const {return TESTBIT(fStatus, kProcCharge);}
- Bool_t HasProcessCenterPad() const {return TESTBIT(fStatus, kProcCenterPad);}
- Bool_t HasExB() const { return TestBit(kExB);}
- Bool_t HasProcessMean() const {return TESTBIT(fStatus, kProcMean);}
- Bool_t HasProcessSigma() const {return TESTBIT(fStatus, kProcSigma);}
-
- TObjArray* Histos();
-
- Bool_t IsVisual() const {return Bool_t(fCanvas);}
- Bool_t IsSaveAs() const {return TESTBIT(fStatus, kSaveAs);}
-
- Bool_t PostProcess();
- Bool_t SetExB(Int_t det=-1, Int_t c = 70, Int_t r = 7);
- void SetVisual();
- void SetProcessCharge(Bool_t v = kTRUE) {v ? SETBIT(fStatus, kProcCharge) : CLRBIT(fStatus, kProcCharge);}
- void SetProcessCenterPad(Bool_t v = kTRUE) {v ? SETBIT(fStatus, kProcCenterPad) : CLRBIT(fStatus, kProcCenterPad);}
- void SetProcessMean(Bool_t v = kTRUE) {v ? SETBIT(fStatus, kProcMean) : CLRBIT(fStatus, kProcMean);}
- void SetProcessSigma(Bool_t v = kTRUE) {v ? SETBIT(fStatus, kProcSigma) : CLRBIT(fStatus, kProcSigma);}
- void SetSaveAs(Bool_t v = kTRUE) {v ? SETBIT(fStatus, kSaveAs) : CLRBIT(fStatus, kSaveAs);}
+ Bool_t GetRefFigure(Int_t ifig);
+ Bool_t HasProcess(EResultContainer bit) const {return TESTBIT(fStatus, bit);}
+ Bool_t HasExB() const { return TestBit(kExB);}
+
+ TObjArray* Histos();
+
+ Bool_t IsVisual() const {return Bool_t(fCanvas);}
+ Bool_t IsSaveAs() const {return TestBit(kSaveAs);}
+
+ Bool_t PostProcess();
+ Bool_t SetExB(Int_t det=-1, Int_t c = 70, Int_t r = 7);
+ void SetVisual();
+ void SetProcess(EResultContainer bit, Bool_t v = kTRUE) {v ? SETBIT(fStatus, bit) : CLRBIT(fStatus, bit);}
+ void SetSaveAs(Bool_t v = kTRUE) {SetBit(kSaveAs, v);}
+ inline void ResetProcess();
void Terminate(Option_t *){};
Short_t fDet; // detector (-1 for all)
Float_t fExB; // tg of the Lorentz angle
Float_t fVdrift; // mean drift velocity
+ static const Float_t fgkTimeBinLength;// time bin length (invers of sampling frequency)
ClassDef(AliTRDclusterResolution, 1) // cluster resolution
};
}
return fVdrift;
}
+
+//___________________________________________________
+inline void AliTRDclusterResolution::ResetProcess()
+{
+ CLRBIT(fStatus, kQRes);
+ CLRBIT(fStatus, kCenter);
+ CLRBIT(fStatus, kSigm);
+ CLRBIT(fStatus, kMean);
+}
+
#endif
UChar_t s;
Int_t pdg = fMC->GetPDG(), det=-1;
Int_t label = fMC->GetLabel();
- Double_t x, y, z, pt, dydx, dzdx;
+ Double_t xAnode, x, y, z, pt, dydx, dzdx;
Float_t p, pt0, x0, y0, z0, dx, dy, dz, dydx0, dzdx0;
Double_t covR[3]/*, cov[3]*/;
//radial shift with respect to the MC reference (radial position of the pad plane)
x= fTracklet->GetX();
if(!fMC->GetDirections(x0, y0, z0, dydx0, dzdx0, pt0, s)) continue;
+ xAnode = fTracklet->GetX0();
+
// MC track position at reference radial position
dx = x0 - x;
if(fDebugLevel>=1){
Float_t zt = z0 - dx*dzdx0;
p = pt0*(1.+dzdx0*dzdx0); // pt -> p
- // add Kalman residuals for y, z and pt
- dx = fTracklet->GetX0() - x;
+ // Kalman position at reference radial position
+ dx = xAnode - x;
y = fTracklet->GetYref(0) - dx*fTracklet->GetYref(1);
dy = yt - y;
z = fTracklet->GetZref(0) - dx*fTracklet->GetZref(1);
// recalculate tracklet based on the MC info
AliTRDseedV1 tt(*fTracklet);
- tt.SetZref(0, z0 - (x0-tt.GetX0())*dzdx0);
+ tt.SetZref(0, z0 - (x0-xAnode)*dzdx0);
tt.SetZref(1, dzdx0);
tt.Fit(kTRUE);
x= tt.GetX();y= tt.GetY();z= tt.GetZ();
clInfo->SetGlobalPosition(yt, zt, dydx0, dzdx0);
clInfo->SetResolution(dy);
clInfo->SetAnisochronity(d);
- clInfo->SetDriftLength(dx-.5*AliTRDgeometry::CamHght());
+ clInfo->SetDriftLength(((c->GetPadTime()+0.5)*.1)*1.5);
+ //dx-.5*AliTRDgeometry::CamHght());
clInfo->SetTilt(tilt);
// Fill Debug Tree
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff