#include "TMath.h"
#include "TGeoManager.h"
#include "TTreeStream.h"
#include "TGraphErrors.h"
#include "AliLog.h"
#include "AliMathBase.h"
#include "AliRieman.h"
#include "AliCDBManager.h"
#include "AliTRDReconstructor.h"
#include "AliTRDpadPlane.h"
#include "AliTRDtransform.h"
#include "AliTRDcluster.h"
#include "AliTRDseedV1.h"
#include "AliTRDtrackV1.h"
#include "AliTRDcalibDB.h"
#include "AliTRDchamberTimeBin.h"
#include "AliTRDtrackingChamber.h"
#include "AliTRDtrackerV1.h"
#include "AliTRDrecoParam.h"
#include "AliTRDCommonParam.h"
#include "AliTRDtrackletOflHelper.h"
#include "Cal/AliTRDCalTrkAttach.h"
#include "Cal/AliTRDCalPID.h"
#include "Cal/AliTRDCalROC.h"
#include "Cal/AliTRDCalDet.h"
class AliTracker;
ClassImp(AliTRDseedV1)
AliTRDseedV1::AliTRDseedV1(Int_t det)
:AliTRDtrackletBase()
,fkReconstructor(NULL)
,fClusterIter(NULL)
,fExB(0.)
,fVD(0.)
,fT0(0.)
,fS2PRF(0.)
,fDiffL(0.)
,fDiffT(0.)
,fClusterIdx(0)
,fErrorMsg(0)
,fN(0)
,fDet(det)
,fPt(0.)
,fdX(0.)
,fX0(0.)
,fX(0.)
,fY(0.)
,fZ(0.)
,fS2Y(0.)
,fS2Z(0.)
,fChi2(0.)
{
memset(fIndexes,0xFF,kNclusters*sizeof(fIndexes[0]));
memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
memset(fPad, 0, 4*sizeof(Float_t));
fYref[0] = 0.; fYref[1] = 0.;
fZref[0] = 0.; fZref[1] = 0.;
fYfit[0] = 0.; fYfit[1] = 0.;
fZfit[0] = 0.; fZfit[1] = 0.;
memset(fdEdx, 0, kNdEdxSlices*sizeof(Float_t));
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
fLabels[0]=-1; fLabels[1]=-1;
fLabels[2]=0;
memset(fRefCov, 0, 7*sizeof(Double_t));
fC[0] = 0.; fC[1] = 0.;
fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
SetStandAlone(kFALSE);
}
AliTRDseedV1::AliTRDseedV1(const AliTRDseedV1 &ref)
:AliTRDtrackletBase((AliTRDtrackletBase&)ref)
,fkReconstructor(NULL)
,fClusterIter(NULL)
,fExB(0.)
,fVD(0.)
,fT0(0.)
,fS2PRF(0.)
,fDiffL(0.)
,fDiffT(0.)
,fClusterIdx(0)
,fErrorMsg(0)
,fN(0)
,fDet(-1)
,fPt(0.)
,fdX(0.)
,fX0(0.)
,fX(0.)
,fY(0.)
,fZ(0.)
,fS2Y(0.)
,fS2Z(0.)
,fChi2(0.)
{
if(this != &ref){
ref.Copy(*this);
}
SetBit(kOwner, kFALSE);
SetStandAlone(ref.IsStandAlone());
}
AliTRDseedV1& AliTRDseedV1::operator=(const AliTRDseedV1 &ref)
{
if(this != &ref){
ref.Copy(*this);
}
SetBit(kOwner, kFALSE);
return *this;
}
AliTRDseedV1::~AliTRDseedV1()
{
if(IsOwner()) {
for(int itb=0; itb<kNclusters; itb++){
if(!fClusters[itb]) continue;
delete fClusters[itb];
fClusters[itb] = NULL;
}
}
}
void AliTRDseedV1::Copy(TObject &ref) const
{
AliTRDseedV1 &target = (AliTRDseedV1 &)ref;
target.fkReconstructor = fkReconstructor;
target.fClusterIter = NULL;
target.fExB = fExB;
target.fVD = fVD;
target.fT0 = fT0;
target.fS2PRF = fS2PRF;
target.fDiffL = fDiffL;
target.fDiffT = fDiffT;
target.fClusterIdx = 0;
target.fErrorMsg = fErrorMsg;
target.fN = fN;
target.fDet = fDet;
target.fPt = fPt;
target.fdX = fdX;
target.fX0 = fX0;
target.fX = fX;
target.fY = fY;
target.fZ = fZ;
target.fS2Y = fS2Y;
target.fS2Z = fS2Z;
target.fChi2 = fChi2;
memcpy(target.fIndexes, fIndexes, kNclusters*sizeof(Int_t));
memcpy(target.fClusters, fClusters, kNclusters*sizeof(AliTRDcluster*));
memcpy(target.fPad, fPad, 4*sizeof(Float_t));
target.fYref[0] = fYref[0]; target.fYref[1] = fYref[1];
target.fZref[0] = fZref[0]; target.fZref[1] = fZref[1];
target.fYfit[0] = fYfit[0]; target.fYfit[1] = fYfit[1];
target.fZfit[0] = fZfit[0]; target.fZfit[1] = fZfit[1];
memcpy(target.fdEdx, fdEdx, kNdEdxSlices*sizeof(Float_t));
memcpy(target.fProb, fProb, AliPID::kSPECIES*sizeof(Float_t));
memcpy(target.fLabels, fLabels, 3*sizeof(Int_t));
memcpy(target.fRefCov, fRefCov, 7*sizeof(Double_t));
target.fC[0] = fC[0]; target.fC[1] = fC[1];
memcpy(target.fCov, fCov, 3*sizeof(Double_t));
TObject::Copy(ref);
}
void AliTRDseedV1::Init(const AliRieman *rieman)
{
fZref[0] = rieman->GetZat(fX0);
fZref[1] = rieman->GetDZat(fX0);
fYref[0] = rieman->GetYat(fX0);
fYref[1] = rieman->GetDYat(fX0);
if(fkReconstructor && fkReconstructor->IsHLT()){
fRefCov[0] = 1;
fRefCov[2] = 10;
}else{
fRefCov[0] = rieman->GetErrY(fX0);
fRefCov[2] = rieman->GetErrZ(fX0);
}
fC[0] = rieman->GetC();
fChi2 = rieman->GetChi2();
}
Bool_t AliTRDseedV1::Init(const AliTRDtrackV1 *track)
{
Double_t y, z;
if(!track->GetProlongation(fX0, y, z)) return kFALSE;
Update(track);
return kTRUE;
}
void AliTRDseedV1::Reset(Option_t *opt)
{
for(Int_t ic=kNclusters; ic--;) fIndexes[ic] = -1;
memset(fClusters, 0, kNclusters*sizeof(AliTRDcluster*));
fN=0; SetBit(kRowCross, kFALSE);
if(strcmp(opt, "c")==0) return;
fExB=0.;fVD=0.;fT0=0.;fS2PRF=0.;
fDiffL=0.;fDiffT=0.;
fClusterIdx=0;
fErrorMsg = 0;
fDet=-1;
fPt=0.;
fdX=0.;fX0=0.; fX=0.; fY=0.; fZ=0.;
fS2Y=0.; fS2Z=0.;
fC[0]=0.; fC[1]=0.;
fChi2 = 0.;
memset(fPad, 0, 4*sizeof(Float_t));
fYref[0] = 0.; fYref[1] = 0.;
fZref[0] = 0.; fZref[1] = 0.;
fYfit[0] = 0.; fYfit[1] = 0.;
fZfit[0] = 0.; fZfit[1] = 0.;
memset(fdEdx, 0, kNdEdxSlices*sizeof(Float_t));
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) fProb[ispec] = -1.;
fLabels[0]=-1; fLabels[1]=-1;
fLabels[2]=0;
memset(fRefCov, 0, 7*sizeof(Double_t));
fCov[0] = 4.e-4; fCov[1] = 0.; fCov[2] = 5.3;
}
void AliTRDseedV1::Update(const AliTRDtrackV1 *trk)
{
Double_t fSnp = trk->GetSnp();
Double_t fTgl = trk->GetTgl();
fPt = trk->Pt();
Double_t norm =1./TMath::Sqrt((1.-fSnp)*(1.+fSnp));
fYref[1] = fSnp*norm;
fZref[1] = fTgl*norm;
SetCovRef(trk->GetCovariance());
Double_t dx = trk->GetX() - fX0;
fYref[0] = trk->GetY() - dx*fYref[1];
fZref[0] = trk->GetZ() - dx*fZref[1];
}
void AliTRDseedV1::UpdateUsed()
{
Int_t nused = 0, nshared = 0;
for (Int_t i = kNclusters; i--; ) {
if (!fClusters[i]) continue;
if(fClusters[i]->IsUsed()){
nused++;
} else if(fClusters[i]->IsShared()){
if(IsStandAlone()) nused++;
else nshared++;
}
}
SetNUsed(nused);
SetNShared(nshared);
}
void AliTRDseedV1::UseClusters()
{
AliTRDcluster **c = &fClusters[0];
for (Int_t ic=kNclusters; ic--; c++) {
if(!(*c)) continue;
if(IsStandAlone()){
if((*c)->IsShared() || (*c)->IsUsed()){
if((*c)->IsShared()) SetNShared(GetNShared()-1);
else SetNUsed(GetNUsed()-1);
(*c) = NULL;
fIndexes[ic] = -1;
SetN(GetN()-1);
continue;
}
} else {
if((*c)->IsUsed() || IsKink()){
(*c)->SetShared();
continue;
}
}
(*c)->Use();
}
}
void AliTRDseedV1::CookdEdx(Int_t nslices)
{
memset(fdEdx, 0, kNdEdxSlices*sizeof(Float_t));
const Double_t kDriftLength = (.5 * AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick());
AliTRDcluster *c(NULL);
for(int ic=0; ic<AliTRDtrackerV1::GetNTimeBins(); ic++){
if(!(c = fClusters[ic]) && !(c = fClusters[ic+kNtb])) continue;
Float_t dx = TMath::Abs(fX0 - c->GetX());
Int_t slice;
if(dx<kDriftLength){
slice = Int_t(dx * nslices / kDriftLength);
} else slice = c->GetX() < fX0 ? nslices-1 : 0;
Float_t w = 1.;
fdEdx[slice] += w * GetdQdl(ic);
}
}
void AliTRDseedV1::CookLabels()
{
Int_t labels[200];
Int_t out[200];
Int_t nlab = 0;
for (Int_t i = 0; i < kNclusters; 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++;
}
}
}
fLabels[2] = AliMathBase::Freq(nlab,labels,out,kTRUE);
fLabels[0] = out[0];
if ((fLabels[2] > 1) && (out[3] > 1)) fLabels[1] = out[2];
}
Float_t AliTRDseedV1::GetAnodeWireOffset(Float_t zt)
{
Float_t d = fPad[3] - zt;
if(d<0.){
AliError(Form("Fail AnodeWireOffset calculation z0[%+7.2f] zt[%+7.2f] d[%+7.2f].", fPad[3], zt, d));
return 0.125;
}
d -= ((Int_t)(2 * d)) / 2.0;
if(d > 0.25) d = 0.5 - d;
return d;
}
Float_t AliTRDseedV1::GetCharge(Bool_t useOutliers) const
{
AliTRDcluster *c(NULL); Float_t qt(0.);
for(int ic=0; ic<kNclusters; ic++){
if(!(c=fClusters[ic])) continue;
if(!c->IsInChamber() && !useOutliers) continue;
qt += TMath::Abs(c->GetQ());
}
return qt;
}
Int_t AliTRDseedV1::GetChargeGaps(Float_t sz[kNtb], Float_t pos[kNtb], Int_t isz[kNtb]) const
{
Bool_t gap(kFALSE);
Int_t n(0);
Int_t ipos[kNtb]; memset(isz, 0, kNtb*sizeof(Int_t));memset(ipos, 0, kNtb*sizeof(Int_t));
for(int ic(0); ic<kNtb; ic++){
if(fClusters[ic] || fClusters[ic+kNtb]){
if(gap) n++;
continue;
}
gap = kTRUE;
isz[n]++;
ipos[n] = ic;
}
if(!n) return 0;
AliTRDcluster fake;
for(Int_t igap(0); igap<n; igap++){
sz[igap] = isz[igap]*fVD/AliTRDCommonParam::Instance()->GetSamplingFrequency();
fake.SetPadTime(ipos[igap]);
pos[igap] = fake.GetXloc(fT0, fVD);
if(isz[igap]>1){
fake.SetPadTime(ipos[igap]-isz[igap]+1);
pos[igap] += fake.GetXloc(fT0, fVD);
pos[igap] /= 2.;
}
}
return n;
}
Double_t AliTRDseedV1::EstimatedCrossPoint(AliTRDpadPlane *pp, Float_t bz)
{
Int_t row[] = {-1, -1};
Double_t zoff(0.5 * (pp->GetRow0() + pp->GetRowEnd())), sx(0.), mean(0.5*pp->GetNrows()-0.5);
AliTRDcluster *c(NULL);
fS2Y = 0.;
if(!IsRowCross()){
for(int ic=0; ic<kNtb; ic++){
if(!(c=fClusters[ic])) continue;
if(!c->IsInChamber()) continue;
row[0] = c->GetPadRow();
fZfit[0] = Int_t(mean-row[0])*pp->GetLengthIPad() +
0.5*(mean-row[0]>0.?1.:-1.)*(row[0]>0&&row[0]<pp->GetNrows()-1?pp->GetLengthIPad():pp->GetLengthOPad());
break;
}
} else {
Float_t tbm[2] = {0.};
Int_t tb[kNtb]={0},
nc[2] = {0},
mc(0);
Bool_t w[2] = {kFALSE, kFALSE};
for(int ic(0); ic<kNtb; ic++){
tb[ic]= -1;
if(!(c=fClusters[ic]) || !c->IsInChamber()) continue;
if(row[0]<0) row[0] = c->GetPadRow();
tb[nc[0]++] = ic; tbm[0] += ic;
}
if(nc[0]>2){
tbm[0] /= nc[0];
w[0] = kTRUE;
}
for(int ic(kNtb), jc(0); ic<kNclusters; ic++, jc++){
if(!(c=fClusters[ic]) || !c->IsInChamber()) continue;
if(row[1]<0) row[1] = c->GetPadRow();
tbm[1] += jc; nc[1]++;
for(Int_t kc(0); kc<nc[0]; kc++)
if(tb[kc]==jc){
tb[kc] += 100;
mc++;
break;
}
}
if(nc[1]>2){
tbm[1] /= nc[1];
w[1] = kTRUE;
}
if(!w[0] && !w[1]){
AliError("Too few clusters to estimate tracklet.");
return -1;
}
if(!w[0] || !w[1]){
SetBit(kRowCross, kFALSE);
if(w[1]) row[0] = row[1];
fZfit[0] = Int_t(mean-row[0])*pp->GetLengthIPad() +
0.5*(mean-row[0]>0.?1.:-1.)*(row[0]>0&&row[0]<pp->GetNrows()-1?pp->GetLengthIPad():pp->GetLengthOPad());
}else{
fZfit[0] = Int_t(mean-0.5*(row[0]+row[1]))*pp->GetLengthIPad();
Int_t itb(0), dtb(0);
if(!mc) {
itb = Int_t(0.5*(tbm[0] + tbm[1]));
dtb = Int_t(0.5*TMath::Abs(tbm[0] - tbm[1]));
} else {
Double_t rmax(100.); Int_t itbStart(-1), itbStop(0);
for(Int_t jc(0); jc<nc[0]; jc++){
if(tb[jc] < 100) continue;
Int_t ltb(tb[jc]-100);
Double_t r = (1. - ltb/tbm[0])*(1. - ltb/tbm[1]);
if(TMath::Abs(r)<rmax){ rmax = TMath::Abs(r); itb = ltb; }
if(itbStart<0) itbStart = ltb;
itbStop = ltb;
}
dtb = itbStop-itbStart+1;
}
AliTRDCommonParam *cp = AliTRDCommonParam::Instance();
Double_t freq(cp?cp->GetSamplingFrequency():10.);
fS2Y = ((itb-0.5)/freq - fT0 - 0.189)*fVD;
sx = dtb*0.288675134594812921/freq; sx *= sx; sx += 1.56e-2; sx *= fVD*fVD;
}
}
Float_t dx(fX0-fS2Y);
fZfit[1] = (fZfit[0]+zoff)/dx;
UnbiasDZDX(IsRowCross(), bz);
if(IsRowCross()){
const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam();
if(recoParam){
fS2Y += recoParam->GetCorrDZDXxcross()*TMath::Abs(fZfit[1]);
sx += recoParam->GetCorrDZDXxcross()*recoParam->GetCorrDZDXxcross()*GetS2DZDX(fZfit[1]);
}
sx += GetS2XcrossDZDX(TMath::Abs(fZfit[1]));
fZfit[0] += fZfit[1]*fS2Y;
fS2Z = fZfit[1]*fZfit[1]*sx+fS2Y*fS2Y*GetS2DZDX(fZfit[1]);
}
return sx;
}
void AliTRDseedV1::UnbiasDZDX(Bool_t rc, Float_t bz)
{
const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam();
if(!recoParam) return;
fZfit[1] *= recoParam->GetCorrDZDX(rc)-(bz>0?0.01:0.);
if(rc) fZfit[1] += recoParam->GetCorrDZDXbiasRC(fZfit[1]<0);
}
Double_t AliTRDseedV1::UnbiasY(Bool_t rc, Float_t bz)
{
const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam();
if(!recoParam) return 0.;
Double_t par[3]={0.};
Int_t idx(2*(rc?1:0)+Int_t(bz>0));
recoParam->GetYcorrTailCancel(idx, par);
return par[0]*TMath::Sin(par[1]*fYref[1])+par[2];
}
Float_t AliTRDseedV1::GetQperTB(Int_t tb) const
{
Float_t q = 0;
if(fClusters[tb] )
q += TMath::Abs(fClusters[tb]->GetQ());
if(fClusters[tb+kNtb] )
q += TMath::Abs(fClusters[tb+kNtb]->GetQ());
return q/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
}
Float_t AliTRDseedV1::GetdQdl() const
{
Float_t Q = GetCharge(kTRUE);
return Q/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
}
Float_t AliTRDseedV1::GetdQdl(Int_t ic, Float_t *dl) const
{
// #frac{dq}{dl} = #frac{q_{c}}{dx * #sqrt{1 + #(){#frac{dy}{dx}}^{2}_{fit} + #(){#frac{dz}{dx}}^{2}_{ref}}}
// END_LATEX
//<img src="TRD/trackletDQDT.gif">
//End_Html
Float_t dq = 0.;
Bool_t hasClusterInChamber = kFALSE;
if(fClusters[ic] && fClusters[ic]->IsInChamber()){
hasClusterInChamber = kTRUE;
dq += TMath::Abs(fClusters[ic]->GetQ());
}
if(fClusters[ic+kNtb] && fClusters[ic+kNtb]->IsInChamber()){
hasClusterInChamber = kTRUE;
dq += TMath::Abs(fClusters[ic+kNtb]->GetQ());
}
if(!hasClusterInChamber) return 0.;
if(dq<1.e-3) return 0.;
Double_t dx = fdX;
if(ic-1>=0 && ic+1<kNtb){
Float_t x2(0.), x1(0.);
if(fClusters[ic-1] && fClusters[ic-1]->IsInChamber()) x2 = fClusters[ic-1]->GetX();
else if(fClusters[ic-1+kNtb] && fClusters[ic-1+kNtb]->IsInChamber()) x2 = fClusters[ic-1+kNtb]->GetX();
else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x2 = fClusters[ic]->GetX()+fdX;
else x2 = fClusters[ic+kNtb]->GetX()+fdX;
if(fClusters[ic+1] && fClusters[ic+1]->IsInChamber()) x1 = fClusters[ic+1]->GetX();
else if(fClusters[ic+1+kNtb] && fClusters[ic+1+kNtb]->IsInChamber()) x1 = fClusters[ic+1+kNtb]->GetX();
else if(fClusters[ic] && fClusters[ic]->IsInChamber()) x1 = fClusters[ic]->GetX()-fdX;
else x1 = fClusters[ic+kNtb]->GetX()-fdX;
dx = .5*(x2 - x1);
}
dx *= TMath::Sqrt(1. + fYfit[1]*fYfit[1] + fZref[1]*fZref[1]);
if(dl) (*dl) = dx;
if(dx>1.e-9) return dq/dx;
else return 0.;
}
Float_t AliTRDseedV1::GetMomentum(Float_t *err) const
{
// p=#frac{1}{1/p_{t}} #sqrt{1+tgl^{2}}
// END_LATEX
// #sigma_{p}^{2} = (#frac{dp}{dp_{t}})^{2} #sigma_{p_{t}}^{2}+(#frac{dp}{dtgl})^{2} #sigma_{tgl}^{2}+2#frac{dp}{dp_{t}}#frac{dp}{dtgl} cov(tgl,1/p_{t})
// END_LATEX
// #sigma_{p}^{2} = p^{2}p_{t}^{4}tgl^{2}#sigma_{tgl}^{2}-2p^{2}p_{t}^{3}tgl cov(tgl,1/p_{t})+p^{2}p_{t}^{2}#sigma_{1/p_{t}}^{2}
// END_LATEX
Double_t p = fPt*TMath::Sqrt(1.+fZref[1]*fZref[1]);
if(err){
Double_t p2 = p*p;
Double_t tgl2 = fZref[1]*fZref[1];
Double_t pt2 = fPt*fPt;
Double_t s2 =
p2*tgl2*pt2*pt2*fRefCov[4]
-2.*p2*fZref[1]*fPt*pt2*fRefCov[5]
+p2*pt2*fRefCov[6];
(*err) = TMath::Sqrt(s2);
}
return p;
}
Int_t AliTRDseedV1::GetTBoccupancy() const
{
Int_t n(0);
for(int ic(0); ic<kNtb; ic++){
if(!fClusters[ic] && !fClusters[ic+kNtb]) continue;
n++;
}
return n;
}
Int_t AliTRDseedV1::GetTBcross() const
{
if(!IsRowCross()) return 0;
Int_t n(0);
for(int ic(0); ic<kNtb; ic++){
if(fClusters[ic] && fClusters[ic+kNtb]) n++;
}
return n;
}
Float_t* AliTRDseedV1::GetProbability(Bool_t force)
{
if(!force) return &fProb[0];
if(!CookPID()) return NULL;
return &fProb[0];
}
Bool_t AliTRDseedV1::CookPID()
{
AliWarning(Form("Obsolete function. Use AliTRDPIDResponse::GetResponse() instead."));
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
AliError("No access to calibration data");
return kFALSE;
}
if (!fkReconstructor) {
AliError("Reconstructor not set.");
return kFALSE;
}
const AliTRDCalPID *pd = calibration->GetPIDObject(fkReconstructor->GetPIDMethod());
if (!pd) {
AliError("No access to AliTRDCalPID object");
return kFALSE;
}
Float_t length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/ TMath::Sqrt((1.0 - GetSnp()*GetSnp()) / (1.0 + GetTgl()*GetTgl()));
CookdEdx(AliTRDCalPID::kNSlicesNN);
AliDebug(4, Form("p=%6.4f[GeV/c] dEdx{%7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f %7.2f} l=%4.2f[cm]", GetMomentum(), fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7], length));
Bool_t kPIDNN(fkReconstructor->GetPIDMethod()==AliTRDpidUtil::kNN);
for(int ispec=0; ispec<AliPID::kSPECIES; ispec++)
fProb[ispec] = pd->GetProbability(ispec, GetMomentum(), &fdEdx[0], length, kPIDNN?GetPlane():fkReconstructor->GetRecoParam()->GetPIDLQslices());
return kTRUE;
}
Float_t AliTRDseedV1::GetQuality(Bool_t kZcorr) const
{
Float_t zcorr = kZcorr ? GetTilt() * (fZfit[0] - fZref[0]) : 0.;
return
.5 * TMath::Abs(18.0 - GetN())
+ 10.* TMath::Abs(fYfit[1] - fYref[1])
+ 5. * TMath::Abs(fYfit[0] - fYref[0] + zcorr)
+ 2. * TMath::Abs(fZfit[0] - fZref[0]) / GetPadLength();
}
void AliTRDseedV1::GetCovAt(Double_t , Double_t *cov) const
{
// Y = T_{x} X^{T}
// END_LATEX
// C_{Y} = T_{x} C_{X} T_{x}^{T}
// END_LATEX
// T_{x} = (1 x); X=(y0 dy/dx); C_{X}=#(){#splitline{Var(y0) Cov(y0, dy/dx)}{Cov(y0, dy/dx) Var(dy/dx)}}
// END_LATEX
// T_{#alpha} = #(){#splitline{cos(#alpha) __ sin(#alpha)}{-sin(#alpha) __ cos(#alpha)}}; X=(y z); C_{X}=#(){#splitline{Var(y) 0}{0 Var(z)}}
// END_LATEX
// C_{Y}=#frac{1}{1+tg^{2}#alpha} #(){#splitline{(#sigma_{y}^{2}+tg^{2}#alpha#sigma_{z}^{2}) __ tg#alpha(#sigma_{z}^{2}-#sigma_{y}^{2})}{tg#alpha(#sigma_{z}^{2}-#sigma_{y}^{2}) __ (#sigma_{z}^{2}+tg^{2}#alpha#sigma_{y}^{2})}}
// END_LATEX
// C_{Y}=#(){#splitline{#sigma_{y}^{2} __ (#sigma_{z}^{2}-#sigma_{y}^{2})tg#alpha}{((#sigma_{z}^{2}-#sigma_{y}^{2})tg#alpha __ #sigma_{z}^{2}}}
// END_LATEX
Double_t sy2 = fCov[0];
Double_t sz2 = fS2Z;
if(fkReconstructor){
Double_t sys[15]; memset(sys, 0, 15*sizeof(Double_t));
fkReconstructor->GetRecoParam()->GetSysCovMatrix(sys);
}
if(!IsRowCross()){
Double_t t2 = GetTilt()*GetTilt();
Double_t correction = 1./(1. + t2);
cov[0] = (sy2+t2*sz2)*correction;
cov[1] = GetTilt()*(sz2 - sy2)*correction;
cov[2] = (t2*sy2+sz2)*correction;
} else {
cov[0] = sy2; cov[1] = 0.; cov[2] = sz2;
}
AliDebug(4, Form("C(%6.1f %+6.3f %6.1f) RC[%c]", 1.e4*TMath::Sqrt(cov[0]), cov[1], 1.e4*TMath::Sqrt(cov[2]), IsRowCross()?'y':'n'));
}
Int_t AliTRDseedV1::GetCovSqrt(const Double_t * const c, Double_t *d)
{
// C^{1/2} = VD^{1/2}V^{-1}
// END_LATEX
const Double_t kZero(1.e-20);
Double_t l[2],
v[3];
Double_t tr = c[0]+c[2],
det = c[0]*c[2]-c[1]*c[1];
if(TMath::Abs(det)<kZero) return 1;
Double_t dd = TMath::Sqrt(tr*tr - 4*det);
l[0] = .5*(tr + dd*(c[0]>c[2]?-1.:1.));
l[1] = .5*(tr + dd*(c[0]>c[2]?1.:-1.));
if(l[0]<kZero || l[1]<kZero) return 2;
Double_t den = (l[0]-c[0])*(l[0]-c[0])+c[1]*c[1];
if(den<kZero){
v[0] = TMath::Sign(0., c[1]);
v[1] = TMath::Sign(1., (l[0]-c[0]));
v[2] = TMath::Sign(0., c[1]*(l[0]-c[0])*(l[1]-c[2]));
} else {
Double_t tmp = 1./TMath::Sqrt(den);
v[0] = c[1]* tmp;
v[1] = (l[0]-c[0])*tmp;
if(TMath::Abs(l[1]-c[2])<kZero) v[2] = TMath::Sign(v[0]*(l[0]-c[0])/kZero, (l[1]-c[2]));
else v[2] = v[0]*(l[0]-c[0])/(l[1]-c[2]);
}
l[0] = TMath::Sqrt(l[0]); l[1] = TMath::Sqrt(l[1]);
d[0] = v[0]*v[0]*l[0]+v[1]*v[1]*l[1];
d[1] = v[0]*v[1]*l[0]+v[1]*v[2]*l[1];
d[2] = v[1]*v[1]*l[0]+v[2]*v[2]*l[1];
return 0;
}
Double_t AliTRDseedV1::GetCovInv(const Double_t * const c, Double_t *d)
{
Double_t det = c[0]*c[2] - c[1]*c[1];
if(TMath::Abs(det)<1.e-20) return 0.;
Double_t invDet = 1./det;
d[0] = c[2]*invDet;
d[1] =-c[1]*invDet;
d[2] = c[0]*invDet;
return det;
}
UShort_t AliTRDseedV1::GetVolumeId() const
{
for(Int_t ic(0);ic<kNclusters; ic++){
if(fClusters[ic]) return fClusters[ic]->GetVolumeId();
}
return 0;
}
void AliTRDseedV1::Calibrate()
{
AliCDBManager *cdb = AliCDBManager::Instance();
if(cdb->GetRun() < 0){
AliError("OCDB manager not properly initialized");
return;
}
AliTRDcalibDB *calib = AliTRDcalibDB::Instance();
AliTRDCalROC *vdROC = calib->GetVdriftROC(fDet),
*t0ROC = calib->GetT0ROC(fDet);;
const AliTRDCalDet *vdDet = calib->GetVdriftDet();
const AliTRDCalDet *t0Det = calib->GetT0Det();
Int_t col = 70, row = 7;
AliTRDcluster **c = &fClusters[0];
if(GetN()){
Int_t ic = 0;
while (ic<kNclusters && !(*c)){ic++; c++;}
if(*c){
col = (*c)->GetPadCol();
row = (*c)->GetPadRow();
}
}
fT0 = (t0Det->GetValue(fDet) + t0ROC->GetValue(col,row)) / AliTRDCommonParam::Instance()->GetSamplingFrequency();
fVD = vdDet->GetValue(fDet) * vdROC->GetValue(col, row);
fS2PRF = calib->GetPRFWidth(fDet, col, row); fS2PRF *= fS2PRF;
fExB = AliTRDCommonParam::Instance()->GetOmegaTau(fVD);
AliTRDCommonParam::Instance()->GetDiffCoeff(fDiffL,
fDiffT, fVD);
AliDebug(4, Form("Calibration params for Det[%3d] Col[%3d] Row[%2d]\n t0[%f] vd[%f] s2PRF[%f] ExB[%f] Dl[%f] Dt[%f]", fDet, col, row, fT0, fVD, fS2PRF, fExB, fDiffL, fDiffT));
SetBit(kCalib, kTRUE);
}
void AliTRDseedV1::SetOwner()
{
if(TestBit(kOwner)) return;
for(int ic=0; ic<kNclusters; ic++){
if(!fClusters[ic]) continue;
fClusters[ic] = new AliTRDcluster(*fClusters[ic]);
}
SetBit(kOwner);
}
void AliTRDseedV1::SetPadPlane(AliTRDpadPlane * const p)
{
fPad[0] = p->GetLengthIPad();
fPad[1] = p->GetWidthIPad();
fPad[2] = TMath::Tan(TMath::DegToRad()*p->GetTiltingAngle());
fPad[3] = p->GetRow0() + p->GetAnodeWireOffset();
}
Bool_t AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_t tilt, Bool_t chgPos, Int_t ev)
{
// r_{y} = 3*#sqrt{12*(#sigma^{2}_{Trk}(y) + #frac{#sigma^{2}_{cl}(y) + tg^{2}(#alpha_{L})#sigma^{2}_{cl}(z)}{1+tg^{2}(#alpha_{L})})}
// r_{z} = 1.5*L_{pad}
// END_LATEX
const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam();
if(!recoParam){
AliError("Tracklets can not be used without a valid RecoParam.");
return kFALSE;
}
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if (!calibration) {
AliError("No access to calibration data");
return kFALSE;
}
const AliTRDCalTrkAttach *attach = calibration->GetAttachObject();
if (!attach) {
AliError("No usable AttachClusters calib object.");
return kFALSE;
}
Int_t t0 = 14;
Int_t kClmin = Int_t(recoParam->GetFindableClusters()*AliTRDtrackerV1::GetNTimeBins());
Int_t kTBmin = 4;
Double_t sysCov[5]; recoParam->GetSysCovMatrix(sysCov);
Double_t s2yTrk= fRefCov[0],
s2yCl = 0.,
s2zCl = GetPadLength()*GetPadLength()/12.,
syRef = TMath::Sqrt(s2yTrk),
t2 = GetTilt()*GetTilt();
const Double_t kroady = 3.;
const Double_t kroadz = GetPadLength() * recoParam->GetRoadzMultiplicator() + 1.;
Short_t sig[] = {0, 0, 10, 30, 10, 0,0};
AliTRDcluster cp(fDet, 6, 75, 0, sig, 0);
if(fkReconstructor->IsHLT()) cp.SetRPhiMethod(AliTRDcluster::kCOG);
if(!IsCalibrated()) Calibrate();
AliDebug(4, Form("\n syTrk[cm]=%4.2f dydxTrk[deg]=%+6.2f Chg[%c] rY[cm]=%4.2f rZ[cm]=%5.2f TC[%c]", syRef, TMath::ATan(fYref[1])*TMath::RadToDeg(), chgPos?'+':'-', kroady, kroadz, tilt?'y':'n'));
Double_t phiTrk(TMath::ATan(fYref[1])),
thtTrk(TMath::ATan(fZref[1]));
const Int_t kNrows = 16;
const Int_t kNcls = 3*kNclusters;
TObjArray clst[kNrows];
Bool_t blst[kNrows][kNcls];
Double_t cond[4],
dx, dy, dz,
yt, zt,
zc[kNrows],
xres[kNrows][kNcls], yres[kNrows][kNcls], zres[kNrows][kNcls], s2y[kNrows][kNcls];
Int_t idxs[kNrows][kNcls], ncl[kNrows], ncls = 0;
memset(ncl, 0, kNrows*sizeof(Int_t));
memset(zc, 0, kNrows*sizeof(Double_t));
memset(idxs, 0, kNrows*kNcls*sizeof(Int_t));
memset(xres, 0, kNrows*kNcls*sizeof(Double_t));
memset(yres, 0, kNrows*kNcls*sizeof(Double_t));
memset(zres, 0, kNrows*kNcls*sizeof(Double_t));
memset(s2y, 0, kNrows*kNcls*sizeof(Double_t));
memset(blst, 0, kNrows*kNcls*sizeof(Bool_t));
Double_t roady(0.), s2Mean(0.); Int_t ns2Mean(0);
AliTRDcluster *c(NULL);
AliTRDchamberTimeBin *layer(NULL);
Bool_t kBUFFER = kFALSE;
for (Int_t it = 0; it < kNtb; it++) {
if(!(layer = chamber->GetTB(it))) continue;
if(!Int_t(*layer)) continue;
dx = fX0 - layer->GetX();
yt = fYref[0] - fYref[1] * dx;
zt = fZref[0] - fZref[1] * dx;
cp.SetLocalTimeBin(it);
cp.SetSigmaY2(0.02, fDiffT, fExB, dx, -1., fYref[1]);
s2yCl = cp.GetSigmaY2() + sysCov[0]; if(!tilt) s2yCl = (s2yCl + t2*s2zCl)/(1.+t2);
if(TMath::Abs(it-12)<7){ s2Mean += cp.GetSigmaY2(); ns2Mean++;}
roady = TMath::Min(3.*TMath::Sqrt(12.*(s2yTrk + s2yCl)), kroady);
AliDebug(5, Form("\n"
" %2d xd[cm]=%6.3f yt[cm]=%7.2f zt[cm]=%8.2f\n"
" syTrk[um]=%6.2f syCl[um]=%6.2f syClTlt[um]=%6.2f\n"
" Ry[mm]=%f"
, it, dx, yt, zt
, 1.e4*TMath::Sqrt(s2yTrk), 1.e4*TMath::Sqrt(cp.GetSigmaY2()+sysCov[0]), 1.e4*TMath::Sqrt(s2yCl)
, 1.e1*roady));
cond[0] = yt; cond[2] = roady;
cond[1] = zt; cond[3] = kroadz;
Int_t n=0, idx[6]; layer->GetClusters(cond, idx, n, 6);
for(Int_t ic = n; ic--;){
c = (*layer)[idx[ic]];
dx = fX0 - c->GetX();
yt = fYref[0] - fYref[1] * dx;
zt = fZref[0] - fZref[1] * dx;
dz = zt - c->GetZ();
dy = yt - (c->GetY() + (tilt ? (GetTilt() * dz) : 0.));
Int_t r = c->GetPadRow();
clst[r].AddAtAndExpand(c, ncl[r]);
blst[r][ncl[r]] = kTRUE;
idxs[r][ncl[r]] = idx[ic];
zres[r][ncl[r]] = dz/GetPadLength();
yres[r][ncl[r]] = dy;
xres[r][ncl[r]] = dx;
zc[r] = c->GetZ();
s2y[r][ncl[r]] = TMath::Min(c->GetSigmaY2()+sysCov[0], 0.025);
AliDebug(5, Form(" -> dy[cm]=%+7.4f yc[cm]=%7.2f row[%d] idx[%2d]", dy, c->GetY(), r, ncl[r]));
ncl[r]++; ncls++;
if(ncl[r] >= kNcls) {
AliWarning(Form("Cluster candidates row[%d] reached buffer limit[%d]. Some may be lost.", r, kNcls));
kBUFFER = kTRUE;
break;
}
}
if(kBUFFER) break;
}
if(ncls<kClmin){
AliDebug(1, Form("CLUSTERS FOUND %d LESS THAN THRESHOLD %d.", ncls, kClmin));
SetErrorMsg(kAttachClFound);
for(Int_t ir(kNrows);ir--;) clst[ir].Clear();
return kFALSE;
}
if(ns2Mean<kTBmin){
AliDebug(1, Form("CLUSTERS IN TimeBins %d LESS THAN THRESHOLD %d.", ns2Mean, kTBmin));
SetErrorMsg(kAttachClFound);
for(Int_t ir(kNrows);ir--;) clst[ir].Clear();
return kFALSE;
}
s2Mean /= ns2Mean;
Int_t idxRow[kNrows], nrc(0); Double_t zresRow[kNrows];
for(Int_t ir(0); ir<kNrows; ir++){
idxRow[ir]=-1; zresRow[ir] = 999.;
if(!ncl[ir]) continue;
dz = 0.; for(Int_t ic = ncl[ir]; ic--;) dz += zres[ir][ic]; dz/=ncl[ir];
idxRow[nrc] = ir; zresRow[nrc] = TMath::Abs(dz); nrc++;
}
AliDebug(4, Form("Found %d clusters in %d rows. Sorting ...", ncls, nrc));
if(nrc>=2){
if(nrc==2){
if(zresRow[0]>zresRow[1]){
Int_t itmp=idxRow[1]; idxRow[1] = idxRow[0]; idxRow[0] = itmp;
Double_t dtmp=zresRow[1]; zresRow[1] = zresRow[0]; zresRow[0] = dtmp;
}
if(TMath::Abs(idxRow[1] - idxRow[0]) != 1){
SetErrorMsg(kAttachRowGap);
AliDebug(2, Form("Rows attached not continuous. Select first candidate.\n"
" row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f",
idxRow[0], ncl[idxRow[0]], zresRow[0], idxRow[1], idxRow[1]<0?0:ncl[idxRow[1]], zresRow[1]));
nrc=1; idxRow[1] = -1; zresRow[1] = 999.;
}
} else {
Int_t idx0[kNrows];
TMath::Sort(nrc, zresRow, idx0, kFALSE);
nrc = 3;
Int_t iatmp[] = {-1, -1, -1}; Double_t datmp[] = {999., 999., 999.};
for(Int_t irc(0); irc<nrc; irc++){
iatmp[irc] = idxRow[idx0[irc]];
datmp[irc] = zresRow[idx0[irc]];
}
idxRow[0] = iatmp[0]; zresRow[0] = datmp[0];
idxRow[1] = iatmp[1]; zresRow[1] = datmp[1];
idxRow[2] = iatmp[2]; zresRow[2] = datmp[2];
if(TMath::Abs(idxRow[1] - idxRow[0]) != 1){
SetErrorMsg(kAttachRowGap);
AliDebug(2, Form("Rows attached not continuous. Turn on selection.\n"
"row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f\n"
"row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f\n"
"row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f",
idxRow[0], ncl[idxRow[0]], zresRow[0],
idxRow[1], ncl[idxRow[1]], zresRow[1],
idxRow[2], ncl[idxRow[2]], zresRow[2]));
if(TMath::Abs(idxRow[0] - idxRow[2]) == 1){
AliDebug(2, "Solved ! Remove second candidate.");
nrc = 2;
idxRow[1] = idxRow[2]; zresRow[1] = zresRow[2];
idxRow[2] = -1; zresRow[2] = 999.;
} else if(TMath::Abs(idxRow[1] - idxRow[2]) == 1){
if(ncl[idxRow[1]]+ncl[idxRow[2]] > ncl[idxRow[0]]){
AliDebug(2, "Solved ! Remove first candidate.");
nrc = 2;
idxRow[0] = idxRow[1]; zresRow[0] = zresRow[1];
idxRow[1] = idxRow[2]; zresRow[1] = zresRow[2];
} else {
AliDebug(2, "Solved ! Remove second and third candidate.");
nrc = 1;
idxRow[1] = -1; zresRow[1] = 999.;
idxRow[2] = -1; zresRow[2] = 999.;
}
} else {
AliDebug(2, "Unsolved !!! Remove second and third candidate.");
nrc = 1;
idxRow[1] = -1; zresRow[1] = 999.;
idxRow[2] = -1; zresRow[2] = 999.;
}
} else {
nrc = 2;
idxRow[2] = -1; zresRow[2] = 999.;
}
}
}
AliDebug(4, Form("Sorted row candidates:\n"
" row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f row[%2d] Ncl[%2d] <dz>[cm]=%+8.2f"
, idxRow[0], ncl[idxRow[0]], zresRow[0], idxRow[1], idxRow[1]<0?0:ncl[idxRow[1]], zresRow[1]));
TTreeSRedirector *pstreamer(NULL);
if((recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming())||
AliTRDReconstructor::GetStreamLevel()>30) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
if(pstreamer){
TVectorD vdy[2], vdx[2], vs2[2];
for(Int_t jr(0); jr<nrc; jr++){
Int_t ir(idxRow[jr]);
vdx[jr].ResizeTo(ncl[ir]); vdy[jr].ResizeTo(ncl[ir]); vs2[jr].ResizeTo(ncl[ir]);
for(Int_t ic(ncl[ir]); ic--;){
vdx[jr](ic) = xres[ir][ic];
vdy[jr](ic) = yres[ir][ic];
vs2[jr](ic) = s2y[ir][ic];
}
}
(*pstreamer) << "AttachClusters4"
<< "r0=" << idxRow[0]
<< "dz0=" << zresRow[0]
<< "dx0=" << &vdx[0]
<< "dy0=" << &vdy[0]
<< "s20=" << &vs2[0]
<< "r1=" << idxRow[1]
<< "dz1=" << zresRow[1]
<< "dx1=" << &vdx[1]
<< "dy1=" << &vdy[1]
<< "s21=" << &vs2[1]
<< "\n";
vdx[0].Clear(); vdy[0].Clear(); vs2[0].Clear();
vdx[1].Clear(); vdy[1].Clear(); vs2[1].Clear();
if(recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 4 ||AliTRDReconstructor::GetStreamLevel()>4){
Int_t idx(idxRow[1]);
if(idx<0){
for(Int_t ir(0); ir<kNrows; ir++){
if(clst[ir].GetEntries()>0) continue;
idx = ir;
break;
}
}
(*pstreamer) << "AttachClusters5"
<< "c0.=" << &clst[idxRow[0]]
<< "c1.=" << &clst[idx]
<< "\n";
}
}
Double_t f[kNcls],
r[2][kNcls],
xm[2][kNcls],
ym[2][kNcls],
sm[2][kNcls],
s[2][kNcls],
p[2][kNcls],
q[2][kNcls];
memset(f, 0, kNcls*sizeof(Double_t));
Int_t index[2][kNcls], n[2][kNcls];
memset(n, 0, 2*kNcls*sizeof(Int_t));
Int_t mts(0), nts[2] = {0, 0};
AliTRDpadPlane *pp(AliTRDtransform::Geometry().GetPadPlane(fDet));
AliTRDtrackletOflHelper helper;
Int_t lyDet(AliTRDgeometry::GetLayer(fDet));
for(Int_t jr(0), n0(0); jr<nrc; jr++){
Int_t ir(idxRow[jr]);
Bool_t kInit(kFALSE);
if(jr==0){
n0 = helper.Init(pp, &clst[ir]); kInit = kTRUE;
if(!n0 || (helper.ClassifyTopology() == AliTRDtrackletOflHelper::kNormal)){
nts[jr] = 1; memset(index[jr], 0, ncl[ir]*sizeof(Int_t));
n[jr][0] = ncl[ir];
}
}
if(!n[jr][0]){
nts[jr] = AliTRDtrackletOflHelper::Segmentation(ncl[ir], xres[ir], yres[ir], index[jr]);
for(Int_t ic(ncl[ir]);ic--;) n[jr][index[jr][ic]]++;
}
mts += nts[jr];
for(Int_t its(0); its<nts[jr]; its++){
if(n[jr][its]<=2) {
xm[jr][its] = 0.;ym[jr][its] = 0.;sm[jr][its] = 0.;
for(Int_t ic(ncl[ir]); ic--;){
if(its != index[jr][ic]) continue;
ym[jr][its] += yres[ir][ic];
xm[jr][its] += xres[ir][ic];
sm[jr][its] += TMath::Sqrt(s2y[ir][ic]);
}
if(n[jr][its]==2){ xm[jr][its] *= 0.5; ym[jr][its] *= 0.5; sm[jr][its] *= 0.5;}
xm[jr][its]= fX0 - xm[jr][its];
r[jr][its] = 0.;
s[jr][its] = 1.e-5;
p[jr][its] = 1.;
q[jr][its] = -1.;
continue;
}
if(!kInit) n0 = helper.Init(pp, &clst[ir], index[jr], its);
Int_t n1 = helper.GetRMS(r[jr][its], ym[jr][its], s[jr][its], fX0);
p[jr][its] = Double_t(n1)/n0;
sm[jr][its] = helper.GetSyMean();
q[jr][its] = helper.GetQ()/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
xm[jr][its] = fX0;
Double_t dxm= fX0 - xm[jr][its];
yt = fYref[0] - fYref[1]*dxm;
zt = fZref[0] - fZref[1]*dxm;
ym[jr][its]+= GetTilt()*(zt - zc[ir]);
r[jr][its] += GetTilt() * fZref[1];
ym[jr][its] = yt - ym[jr][its];
r[jr][its] = (r[jr][its] - fYref[1])/(1+r[jr][its]*fYref[1]);
r[jr][its] = TMath::ATan(r[jr][its]);
if(jr) continue;
f[its] = attach->CookLikelihood(chgPos, lyDet, fPt, phiTrk, n[jr][its], ym[jr][its], r[jr][its]*TMath::RadToDeg(), s[jr][its]/sm[jr][its]);
}
}
AliDebug(4, Form(" Tracklet candidates: row[%2d] = %2d row[%2d] = %2d:", idxRow[0], nts[0], idxRow[1], nts[1]));
if(AliLog::GetDebugLevel("TRD", "AliTRDseedV1")>3){
for(Int_t jr(0); jr<nrc; jr++){
Int_t ir(idxRow[jr]);
for(Int_t its(0); its<nts[jr]; its++){
printf(" segId[%2d] row[%2d] Ncl[%2d] x[cm]=%7.2f dz[pu]=%4.2f dy[mm]=%+7.3f r[deg]=%+6.2f p[%%]=%6.2f s[um]=%7.2f\n",
its, ir, n[jr][its], xm[jr][its], zresRow[jr], 1.e1*ym[jr][its], r[jr][its]*TMath::RadToDeg(), 100.*p[jr][its], 1.e4*s[jr][its]);
}
}
}
if(!pstreamer &&
( (recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 2 && fkReconstructor->IsDebugStreaming()) ||
AliTRDReconstructor::GetStreamLevel()>2 )
) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
if(pstreamer){
TVectorD vidx, vn, vx, vy, vr, vs, vsm, vp, vf;
vidx.ResizeTo(ncl[idxRow[0]]+(idxRow[1]<0?0:ncl[idxRow[1]]));
vn.ResizeTo(mts);
vx.ResizeTo(mts);
vy.ResizeTo(mts);
vr.ResizeTo(mts);
vs.ResizeTo(mts);
vsm.ResizeTo(mts);
vp.ResizeTo(mts);
vf.ResizeTo(mts);
for(Int_t jr(0), jts(0), jc(0); jr<nrc; jr++){
Int_t ir(idxRow[jr]);
for(Int_t its(0); its<nts[jr]; its++, jts++){
vn[jts] = n[jr][its];
vx[jts] = xm[jr][its];
vy[jts] = ym[jr][its];
vr[jts] = r[jr][its];
vs[jts] = s[jr][its];
vsm[jts]= sm[jr][its];
vp[jts] = p[jr][its];
vf[jts] = jr?-1.:f[its];
}
for(Int_t ic(0); ic<ncl[ir]; ic++, jc++) vidx[jc] = index[jr][ic];
}
(*pstreamer) << "AttachClusters3"
<< "idx=" << &vidx
<< "n=" << &vn
<< "x=" << &vx
<< "y=" << &vy
<< "r=" << &vr
<< "s=" << &vs
<< "sm=" << &vsm
<< "p=" << &vp
<< "f=" << &vf
<< "\n";
}
Int_t idx2[kNcls]; memset(idx2, 0, kNcls*sizeof(Int_t));
if(nts[0]>1) TMath::Sort(nts[0], f, idx2);
Int_t is(idx2[0]);
Int_t idxTrklt[kNcls],
kts(0),
nTrklt(n[0][is]);
Double_t fTrklt(f[is]),
rTrklt(r[0][is]),
yTrklt(ym[0][is]),
sTrklt(s[0][is]),
smTrklt(sm[0][is]),
xTrklt(xm[0][is]),
pTrklt(p[0][is]),
qTrklt(q[0][is]);
memset(idxTrklt, 0, kNcls*sizeof(Int_t));
if(f[is]<1.e-2){
AliDebug(1, Form("Seed seg[%d] row[%2d] n[%2d] f[%f]<0.01.", is, idxRow[0], n[0][is], f[is]));
SetErrorMsg(kAttachClAttach);
if(!pstreamer &&
( (recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 1 && fkReconstructor->IsDebugStreaming()) ||
AliTRDReconstructor::GetStreamLevel()>1 )
) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
if(pstreamer){
UChar_t stat(0);
if(IsKink()) SETBIT(stat, 1);
if(IsStandAlone()) SETBIT(stat, 2);
if(IsRowCross()) SETBIT(stat, 3);
SETBIT(stat, 4);
TVectorD vidx; vidx.ResizeTo(1); vidx[0] = is;
(*pstreamer) << "AttachClusters2"
<< "stat=" << stat
<< "ev=" << ev
<< "chg=" << chgPos
<< "det=" << fDet
<< "x0=" << fX0
<< "y0=" << fYref[0]
<< "z0=" << fZref[0]
<< "phi=" << phiTrk
<< "tht=" << thtTrk
<< "pt=" << fPt
<< "s2Trk=" << s2yTrk
<< "s2Cl=" << s2Mean
<< "idx=" << &vidx
<< "n=" << nTrklt
<< "f=" << fTrklt
<< "x=" << xTrklt
<< "y=" << yTrklt
<< "r=" << rTrklt
<< "s=" << sTrklt
<< "sm=" << smTrklt
<< "p=" << pTrklt
<< "q=" << qTrklt
<< "\n";
}
return kFALSE;
}
AliDebug(2, Form("Seed seg[%d] row[%2d] n[%2d] dy[%f] r[%+5.2f] s[%+5.2f] f[%5.3f] q[%6.2f]", is, idxRow[0], n[0][is], ym[0][is], r[0][is]*TMath::RadToDeg(), s[0][is]/sm[0][is], f[is], q[0][is]));
idxTrklt[kts++] = is;
helper.Init(pp, &clst[idxRow[0]], index[0], is);
AliTRDtrackletOflHelper test;
Float_t rcLikelihood(0.); SetBit(kRowCross, kFALSE);
Double_t dyRez[kNcls]; Int_t idx3[kNcls];
Int_t kNSgmDy[2]; attach->GetNsgmDy(kNSgmDy[0], kNSgmDy[1]);
Float_t kLikeMinRelDecrease[2]; attach->GetLikeMinRelDecrease(kLikeMinRelDecrease[0], kLikeMinRelDecrease[1]);
Float_t kRClikeLimit(attach->GetRClikeLimit());
if(nts[0]>1){
Int_t jr(0), ir(idxRow[jr]);
memset(dyRez, 0, nts[jr]*sizeof(Double_t));
for(Int_t jts(1); jts<nts[jr]; jts++) {
Int_t its(idx2[jts]);
Double_t rot(TMath::Tan(r[0][is]));
dyRez[its] = TMath::Abs(ym[0][is] - ym[jr][its] + rot*(xm[0][is]-xm[jr][its]));
}
TMath::Sort(nts[jr], dyRez, idx3, kFALSE);
for (Int_t jts(1); jts<nts[jr]; jts++) {
Int_t its(idx3[jts]);
if(dyRez[its] > kNSgmDy[jr]*smTrklt){
AliDebug(2, Form("Reject seg[%d] row[%2d] n[%2d] dy[%f] > %d*s[%f].", its, idxRow[jr], n[jr][its], dyRez[its], kNSgmDy[jr], kNSgmDy[jr]*smTrklt));
continue;
}
test = helper;
Int_t n0 = test.Expand(&clst[ir], index[jr], its);
Double_t rt, dyt, st, xt, smt, pt, qt, ft;
Int_t n1 = test.GetRMS(rt, dyt, st, fX0);
pt = Double_t(n1)/n0;
smt = test.GetSyMean();
qt = test.GetQ()/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
xt = fX0;
Double_t dxm= fX0 - xt;
yt = fYref[0] - fYref[1]*dxm;
zt = fZref[0] - fZref[1]*dxm;
dyt+= GetTilt()*(zt - zc[idxRow[0]]);
rt += GetTilt() * fZref[1];
dyt = yt - dyt;
rt = (rt - fYref[1])/(1+rt*fYref[1]);
rt = TMath::ATan(rt);
ft = (n0>=2) ? attach->CookLikelihood(chgPos, lyDet, fPt, phiTrk, n0, dyt, rt*TMath::RadToDeg(), st/smt) : 0.;
Bool_t kAccept(ft>=fTrklt*(1.-kLikeMinRelDecrease[jr]));
AliDebug(2, Form("%s seg[%d] row[%2d] n[%2d] dy[%f] r[%+5.2f] s[%+5.2f] f[%f] < %4.2f*F[%f].",
(kAccept?"Adding":"Reject"), its, idxRow[jr], n0, dyt, rt*TMath::RadToDeg(), st/smt, ft, 1.-kLikeMinRelDecrease[jr], fTrklt*(1.-kLikeMinRelDecrease[jr])));
if(kAccept){
idxTrklt[kts++] = its;
nTrklt = n0;
fTrklt = ft;
rTrklt = rt;
yTrklt = dyt;
sTrklt = st;
smTrklt= smt;
xTrklt = xt;
pTrklt = pt;
qTrklt = qt;
helper.Expand(&clst[ir], index[jr], its);
}
}
}
if(nts[1] && (rcLikelihood = zresRow[0]/zresRow[1]) > kRClikeLimit){
Int_t jr(1), ir(idxRow[jr]);
memset(dyRez, 0, nts[jr]*sizeof(Double_t));
Double_t rot(TMath::Tan(r[0][is]));
for(Int_t jts(0); jts<nts[jr]; jts++) {
dyRez[jts] = TMath::Abs(ym[0][is] - ym[jr][jts] + rot*(xm[0][is]-xm[jr][jts]));
}
TMath::Sort(nts[jr], dyRez, idx3, kFALSE);
for (Int_t jts(0); jts<nts[jr]; jts++) {
Int_t its(idx3[jts]);
if(dyRez[its] > kNSgmDy[jr]*smTrklt){
AliDebug(2, Form("Reject seg[%d] row[%2d] n[%2d] dy[%f] > %d*s[%f].", its, idxRow[jr], n[jr][its], dyRez[its], kNSgmDy[jr], kNSgmDy[jr]*smTrklt));
continue;
}
test = helper;
Int_t n0 = test.Expand(&clst[ir], index[jr], its);
Double_t rt, dyt, st, xt, smt, pt, qt, ft;
Int_t n1 = test.GetRMS(rt, dyt, st, fX0);
pt = Double_t(n1)/n0;
smt = test.GetSyMean();
qt = test.GetQ()/TMath::Sqrt(1. + fYref[1]*fYref[1] + fZref[1]*fZref[1]);
xt = fX0;
Double_t dxm= fX0 - xt;
yt = fYref[0] - fYref[1]*dxm;
zt = fZref[0] - fZref[1]*dxm;
dyt+= GetTilt()*(zt - zc[idxRow[0]]);
rt += GetTilt() * fZref[1];
dyt = yt - dyt;
rt = (rt - fYref[1])/(1+rt*fYref[1]);
rt = TMath::ATan(rt);
ft = (n0>=2) ? attach->CookLikelihood(chgPos, lyDet, fPt, phiTrk, n0, dyt, rt*TMath::RadToDeg(), st/smt) : 0.;
Bool_t kAccept(ft>=fTrklt*(1.-kLikeMinRelDecrease[jr]));
AliDebug(2, Form("%s seg[%d] row[%2d] n[%2d] dy[%f] r[%+5.2f] s[%+5.2f] f[%f] < %4.2f*F[%f].",
(kAccept?"Adding":"Reject"), its, idxRow[jr], n0, dyt, rt*TMath::RadToDeg(), st/smt, ft, 1.-kLikeMinRelDecrease[jr], fTrklt*(1.-kLikeMinRelDecrease[jr])));
if(kAccept){
idxTrklt[kts++] = its;
nTrklt = n0;
fTrklt = ft;
rTrklt = rt;
yTrklt = dyt;
sTrklt = st;
smTrklt= smt;
xTrklt = xt;
pTrklt = pt;
qTrklt = qt;
helper.Expand(&clst[ir], index[jr], its);
SetBit(kRowCross, kTRUE);
}
}
}
for(Int_t ir(0); ir<kNrows; ir++) clst[ir].Clear();
if(!pstreamer &&
((recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 1 && fkReconstructor->IsDebugStreaming()) ||
AliTRDReconstructor::GetStreamLevel()>1 )
) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
if(pstreamer){
UChar_t stat(0);
if(IsKink()) SETBIT(stat, 1);
if(IsStandAlone()) SETBIT(stat, 2);
if(IsRowCross()) SETBIT(stat, 3);
TVectorD vidx; vidx.ResizeTo(kts);
for(Int_t its(0); its<kts; its++) vidx[its] = idxTrklt[its];
(*pstreamer) << "AttachClusters2"
<< "stat=" << stat
<< "ev=" << ev
<< "chg=" << chgPos
<< "det=" << fDet
<< "x0=" << fX0
<< "y0=" << fYref[0]
<< "z0=" << fZref[0]
<< "phi=" << phiTrk
<< "tht=" << thtTrk
<< "pt=" << fPt
<< "s2Trk=" << s2yTrk
<< "s2Cl=" << s2Mean
<< "idx=" << &vidx
<< "n=" << nTrklt
<< "q=" << qTrklt
<< "f=" << fTrklt
<< "x=" << xTrklt
<< "y=" << yTrklt
<< "r=" << rTrklt
<< "s=" << sTrklt
<< "sm=" << smTrklt
<< "p=" << pTrklt
<< "\n";
}
Int_t nselected(0), nc(0);
TObjArray *selected(helper.GetClusters());
if(!selected || !(nselected = selected->GetEntriesFast())){
AliError("Cluster candidates missing !!!");
SetErrorMsg(kAttachClAttach);
return kFALSE;
}
for(Int_t ic(0); ic<nselected; ic++){
if(!(c = (AliTRDcluster*)selected->At(ic))) continue;
Int_t it(c->GetPadTime()),
jr(Int_t(helper.GetRow() != c->GetPadRow())),
idx(it+kNtb*jr);
if(fClusters[idx]){
AliDebug(1, Form("Multiple clusters/tb for D[%03d] Tb[%02d] Row[%2d]", fDet, it, c->GetPadRow()));
continue;
}
fIndexes[idx] = chamber->GetTB(it)->GetGlobalIndex(idxs[idxRow[jr]][ic]);
fClusters[idx] = c;
nc++;
}
AliDebug(2, Form("Clusters Found[%2d] Attached[%2d] RC[%c]", nselected, nc, IsRowCross()?'y':'n'));
if (nc < kClmin){
AliDebug(1, Form("NOT ENOUGH CLUSTERS %d ATTACHED TO THE TRACKLET [min %d] FROM FOUND %d.", nc, kClmin, ncls));
SetErrorMsg(kAttachClAttach);
return kFALSE;
}
SetN(nc);
Float_t x[2] = {0.,0.}; Int_t tb[2]={0,0};
for (Int_t it = t0, irp=0; irp<2 && it < AliTRDtrackerV1::GetNTimeBins(); it++) {
if(!fClusters[it]) continue;
x[irp] = fClusters[it]->GetX();
tb[irp] = fClusters[it]->GetLocalTimeBin();
irp++;
}
Int_t dtb = tb[1] - tb[0];
fdX = dtb ? (x[0] - x[1]) / dtb : 0.15;
return kTRUE;
}
void AliTRDseedV1::Bootstrap(const AliTRDReconstructor *rec)
{
fkReconstructor = rec;
AliTRDgeometry g;
SetPadPlane(g.GetPadPlane(fDet));
Int_t n = 0, nshare = 0, nused = 0;
AliTRDcluster **cit = &fClusters[0];
for(Int_t ic = kNclusters; ic--; cit++){
if(!(*cit)) return;
n++;
if((*cit)->IsShared()) nshare++;
if((*cit)->IsUsed()) nused++;
}
SetN(n); SetNUsed(nused); SetNShared(nshare);
Fit();
CookLabels();
GetProbability();
}
Bool_t AliTRDseedV1::Fit(UChar_t opt)
{
// #sigma^{2}_{y} = #sigma^{2}_{PRF} + #frac{x#delta_{t}^{2}}{(1+tg(#alpha_{L}))^{2}} + #frac{x^{2}tg^{2}(#phi-#alpha_{L})tg^{2}(#alpha_{L})}{12}
// END_LATEX
// #sigma_{x|y} = tg(#phi) #sigma_{x}
// END_LATEX
// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
// END_LATEX
// #sigma_{y} = #sigma^{2}_{y_{0}} + 2xcov(y_{0}, dy/dx) + #sigma^{2}_{dy/dx}
// END_LATEX
// x = - cov(y_{0}, dy/dx)/#sigma^{2}_{dy/dx}
// END_LATEX
// #sigma_{y} = #sigma^{2}_{y_{0}} - 2cov^{2}(y_{0}, dy/dx)/#sigma^{2}_{dy/dx} + #sigma^{2}_{dy/dx}
// #sigma_{z} = Pad_{length}/12
// END_LATEX
// #sigma_{z} = tg(#theta) #Delta x_{x_{0}}/12
// END_LATEX
if(!fkReconstructor){
AliError("The tracklet needs the reconstruction setup. Please initialize by SetReconstructor().");
return kFALSE;
}
if(!IsCalibrated()) Calibrate();
if(opt>2){
AliWarning(Form("Option [%d] outside range [0, 2]. Using default",opt));
opt=0;
}
const Int_t kClmin = 8;
const Float_t kScalePulls = 10.;
Double_t y0 = fYref[0];
Double_t dydx = fYref[1];
Double_t z0 = fZref[0];
Double_t dzdx = fZref[1];
AliTRDtrackerV1::AliTRDLeastSquare fitterY;
AliTRDtrackerV1::AliTRDLeastSquare fitterZ;
Double_t qc[kNclusters], xc[kNclusters], yc[kNclusters], zc[kNclusters], sy[kNclusters];
Bool_t tilt(opt==1)
,pseudo(opt==2)
,rc(IsRowCross())
,kDZDX(IsPrimary());
Int_t n(0);
AliTRDcluster *c(NULL), *cc(NULL), **jc = &fClusters[0];
const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam();
const Char_t *tcName[]={"NONE", "FULL", "HALF"};
AliDebug(2, Form("Options : TC[%s] dzdx[%c]", tcName[opt], kDZDX?'Y':'N'));
for (Int_t ic=0; ic<kNclusters; ic++, ++jc) {
xc[ic] = -1.; yc[ic] = 999.; zc[ic] = 999.; sy[ic] = 0.;
if(!(c = (*jc))) continue;
if(!c->IsInChamber()) continue;
if(kDZDX){
fZfit[0] = c->GetZ();
if(rc){
for(Int_t kc=AliTRDseedV1::kNtb; kc<AliTRDseedV1::kNclusters; kc++){
if(!(cc=fClusters[kc])) continue;
if(!cc->IsInChamber()) continue;
fZfit[0] += cc->GetZ(); fZfit[0] *= 0.5;
break;
}
}
fZfit[1] = fZfit[0]/fX0;
if(rc){
fZfit[0] += fZfit[1]*0.5*AliTRDgeometry::CdrHght();
fZfit[1] = fZfit[0]/fX0;
}
kDZDX=kFALSE;
}
qc[n] = TMath::Abs(c->GetQ());
xc[n] = fX0 - c->GetX();
c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], -1., dydx);
c->SetSigmaZ2(fPad[0]*fPad[0]/12.);
sy[n] = TMath::Sqrt(c->GetSigmaY2());
yc[n] = recoParam->UseGAUS() ?
c->GetYloc(y0, sy[n], GetPadWidth()): c->GetY();
zc[n] = c->GetZ();
Float_t correction(0.);
if(tilt) correction = fPad[2]*(xc[n]*dzdx + zc[n] - z0);
else if(pseudo) correction = fPad[2]*(xc[n]*fZfit[1] + zc[n]-fZfit[0]);
yc[n]-=correction;
AliDebug(5, Form(" tb[%2d] dx[%6.3f] y[%6.2f+-%6.3f]", c->GetLocalTimeBin(), xc[n], yc[n], sy[n]));
fitterY.AddPoint(&xc[n], yc[n], sy[n]);
if(rc) fitterZ.AddPoint(&xc[n], qc[n]*(ic<kNtb?1.:-1.), 1.);
n++;
}
if (n < kClmin){
AliDebug(1, Form("Not enough clusters to fit. Clusters: Attached[%d] Fit[%d].", GetN(), n));
SetErrorMsg(kFitCl);
return kFALSE;
}
if(!fitterY.Eval()){
AliDebug(1, "Fit Y failed.");
SetErrorMsg(kFitFailedY);
return kFALSE;
}
fYfit[0] = fitterY.GetFunctionParameter(0);
fYfit[1] = -fitterY.GetFunctionParameter(1);
Double_t p[3];
fitterY.GetCovarianceMatrix(p);
fCov[0] = kScalePulls*p[1];
fCov[1] = kScalePulls*p[2];
fCov[2] = kScalePulls*p[0];
fX = -fCov[1]/fCov[2];
fS2Y = fCov[0] +2.*fX*fCov[1] + fX*fX*fCov[2];
Float_t xs=fX+.5*AliTRDgeometry::CamHght();
if(xs < 0. || xs > AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()){
AliDebug(1, Form("Ref radial position ouside chamber x[%5.2f].", fX));
SetErrorMsg(kFitFailedY);
return kFALSE;
}
if(opt!=1 && IsRowCross()){
if(!fitterZ.Eval()) SetErrorMsg(kFitFailedZ);
if(!HasError(kFitFailedZ) && TMath::Abs(fitterZ.GetFunctionParameter(1))>1.e-10){
}
fS2Z = 0.05+0.4*TMath::Abs(fZref[1]); fS2Z *= fS2Z;
} else {
fS2Z = GetPadLength()*GetPadLength()/12.;
}
return kTRUE;
}
Bool_t AliTRDseedV1::FitRobust(AliTRDpadPlane *pp, TGeoHMatrix *mDet, Float_t bz, Int_t chg, Int_t opt)
{
TTreeSRedirector *pstreamer(NULL);
const AliTRDrecoParam* const recoParam = fkReconstructor->GetRecoParam();
if( (recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming()) ||
AliTRDReconstructor::GetStreamLevel()>3 ) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
Float_t kScalePulls = 1.;
AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
if(!calibration){
AliWarning("No access to calibration data");
} else {
const AliTRDCalTrkAttach *attach = calibration->GetAttachObject();
if(!attach){
AliWarning("No usable AttachClusters calib object.");
} else {
}
SetChmbGood(calibration->IsChamberGood(fDet));
if(!IsChmbGood()) kScalePulls*=10.;
}
AliTRDCommonParam *cp = AliTRDCommonParam::Instance();
Double_t freq(cp?cp->GetSamplingFrequency():10.);
if(EstimatedCrossPoint(pp, bz)<0.) return kFALSE;
Double_t
z0(0.5 * (pp->GetRow0() + pp->GetRowEnd()) + fZfit[0]),
DZ(pp->GetRow0() - pp->GetRowEnd() - pp->GetAnodeWireOffset() + fZfit[0]),
z, d(-1.);
Double_t xc[kNclusters], yc[kNclusters], dz(0.), dzdx(0.),
s2dz(0.), s2dzdx(0.), sy[kNclusters],
s2x((8.33e-2/freq/freq+1.56e-2)*fVD*fVD),
t2(fPad[2]*fPad[2]), loc[3]={0.};
Int_t n(0),
row[]={-1, -1};
Double_t ycorr(UnbiasY(IsRowCross(), bz)),
kS2Ycorr(recoParam->GetS2Ycorr(IsRowCross(), chg>0));
AliTRDcluster *c(NULL), **jc = &fClusters[0];
for(Int_t ic=0; ic<kNtb; ic++, ++jc) {
if(!(c = (*jc))) continue;
if(!c->IsInChamber()) continue;
if(row[0]<0){
row[0] = c->GetPadRow();
z = pp->GetRowPos(row[0]) - 0.5*pp->GetRowSize(row[0]);
switch(opt){
case 0:
dzdx = IsRowCross()?fZfit[1]:0.;
s2dzdx= IsRowCross()?GetS2DZDX(dzdx):0.;
dz = IsRowCross()?(z - z0):0.;
s2dz = IsRowCross()?fS2Z:0.;
break;
case 1:
dzdx = fZref[1];
dz = IsRowCross()?(z - z0):0.;
break;
case 2:
dzdx = fZref[1];
dz = c->GetZ()-fZref[0];
break;
default:
AliError(Form("Wrong option fit %d !", opt));
break;
}
}
Double_t trk[] = {c->GetX(), c->GetY(), c->GetZ()};
mDet->MasterToLocal(trk, loc);
xc[n] = AliTRDgeometry::AnodePos()-loc[0];
yc[n] = loc[1];
yc[n]-= fPad[2]*(dz+xc[n]*dzdx);
yc[n]-= ycorr;
if(IsRowCross()){
d = DZ-xc[n]*dzdx;
d -= ((Int_t)(2 * d)) / 2.0;
if (d > 0.25) d = 0.5 - d;
}
c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], d, fYref[1]);
s2x = c->GetSX(c->GetLocalTimeBin(), d); s2x*=s2x;
sy[n] = c->GetSigmaY2()>0?(TMath::Min(Double_t(c->GetSigmaY2()), 6.4e-3)):6.4e-3;
sy[n]+= t2*(s2dz+xc[n]*xc[n]*s2dzdx+dzdx*dzdx*s2x);
sy[n] = TMath::Sqrt(sy[n]);
n++;
}
for(Int_t ic=kNtb; ic<kNclusters; ic++, ++jc) {
if(!(c = (*jc))) continue;
if(!c->IsInChamber()) continue;
if(row[1]<0){
row[1] = c->GetPadRow();
z = pp->GetRowPos(row[1]) - 0.5*pp->GetRowSize(row[1]);
switch(opt){
case 0:
dz = z - z0;
break;
case 1:
dz = z - z0;
break;
case 2:
dz = c->GetZ()-fZref[0];
break;
default:
AliError(Form("Wrong option fit %d !", opt));
break;
}
}
Double_t trk[] = {c->GetX(), c->GetY(), c->GetZ()};
mDet->MasterToLocal(trk, loc);
xc[n] = AliTRDgeometry::AnodePos()-loc[0];
yc[n] = loc[1];
yc[n]-= fPad[2]*(dz+xc[n]*dzdx);
yc[n]-= ycorr;
d = DZ-xc[n]*dzdx;
d -= ((Int_t)(2 * d)) / 2.0;
if (d > 0.25) d = 0.5 - d;
c->SetSigmaY2(fS2PRF, fDiffT, fExB, xc[n], d, fYref[1]);
s2x = c->GetSX(c->GetLocalTimeBin(), d); s2x*=s2x;
sy[n] = c->GetSigmaY2()>0?(TMath::Min(Double_t(c->GetSigmaY2()), 6.4e-3)):6.4e-3;
sy[n]+= t2*(s2dz+xc[n]*xc[n]*s2dzdx+dzdx*dzdx*s2x);
sy[n] = TMath::Sqrt(sy[n]);
n++;
}
UChar_t status(0);
fX = 0.;
Double_t par[3] = {0.,0.,fX}, cov[3];
if(!AliTRDtrackletOflHelper::Fit(n, xc, yc, sy, par, 1.5, cov)){
AliDebug(1, Form("Tracklet fit failed D[%03d].", fDet));
SetErrorMsg(kFitCl);
return kFALSE;
}
fYfit[0] = par[0] - fX * par[1];
fYfit[1] = -par[1];
fCov[0] = kS2Ycorr*cov[0];
fCov[1] = kScalePulls*cov[2];
fCov[2] = kScalePulls*cov[1];
Float_t xs=fX+.5*AliTRDgeometry::CamHght();
if(xs < 0. || xs > AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()){
AliDebug(1, Form("Ref radial position x[%5.2f] ouside D[%3d].", fX, fDet));
SetErrorMsg(kFitFailedY);
return kFALSE;
}
if(!IsRowCross()){
Double_t padEffLength(fPad[0] - TMath::Abs(dzdx));
fS2Z = padEffLength*padEffLength/12.;
}
AliDebug(2, Form("[I] x[cm]=%6.2f y[cm]=%+5.2f z[cm]=%+6.2f dydx[deg]=%+5.2f", GetX(), GetY(), GetZ(), TMath::ATan(fYfit[1])*TMath::RadToDeg()));
if(pstreamer){
Float_t x= fX0 -fX,
y = GetY(),
yt = fYref[0]-fX*fYref[1];
SETBIT(status, 2);
TVectorD vcov(3); vcov[0]=cov[0];vcov[1]=cov[1];vcov[2]=cov[2];
Double_t sm(0.), chi2(0.), tmp, dy[kNclusters];
for(Int_t ic(0); ic<n; ic++){
sm += sy[ic];
dy[ic] = yc[ic]-(fYfit[0]+(xc[ic]-fX0)*fYfit[1]); tmp = dy[ic]/sy[ic];
chi2 += tmp*tmp;
}
sm /= n; chi2 = TMath::Sqrt(chi2);
Double_t m(0.), s(0.);
AliMathBase::EvaluateUni(n, dy, m, s, 0);
(*pstreamer) << "FitRobust4"
<< "stat=" << status
<< "opt=" << opt
<< "ncl=" << n
<< "det=" << fDet
<< "x0=" << fX0
<< "y0=" << fYfit[0]
<< "x=" << x
<< "y=" << y
<< "dydx=" << fYfit[1]
<< "pt=" << fPt
<< "yt=" << yt
<< "dydxt="<< fYref[1]
<< "cov=" << &vcov
<< "chi2=" << chi2
<< "sm=" << sm
<< "ss=" << s
<< "\n";
}
return kTRUE;
}
void AliTRDseedV1::SetXYZ(TGeoHMatrix *mDet)
{
Double_t loc[] = {AliTRDgeometry::AnodePos(), GetLocalY(), fZfit[0]}, trk[3]={0.};
mDet->LocalToMaster(loc, trk);
fX0 = trk[0];
fY = trk[1];
fZ = trk[2];
return;
}
void AliTRDseedV1::Print(Option_t *o) const
{
AliInfo(Form("Det[%3d] X0[%7.2f] Pad{L[%5.2f] W[%5.2f] Tilt[%+6.2f]}", fDet, fX0, GetPadLength(), GetPadWidth(), GetTilt()));
AliInfo(Form("N[%2d] Nused[%2d] Nshared[%2d] [%d]", GetN(), GetNUsed(), GetNShared(), fN));
AliInfo(Form("FLAGS : RC[%c] Kink[%c] SA[%c]", IsRowCross()?'y':'n', IsKink()?'y':'n', IsStandAlone()?'y':'n'));
AliInfo(Form("CALIB PARAMS : T0[%5.2f] Vd[%5.2f] s2PRF[%5.2f] ExB[%5.2f] Dl[%5.2f] Dt[%5.2f]", fT0, fVD, fS2PRF, fExB, fDiffL, fDiffT));
Double_t cov[3], x=GetX();
GetCovAt(x, cov);
AliInfo(" | x[cm] | y[cm] | z[cm] | dydx | dzdx |");
AliInfo(Form("Fit | %7.2f | %7.2f+-%7.2f | %7.2f+-%7.2f| %5.2f | ----- |", x, GetY(), TMath::Sqrt(cov[0]), GetZ(), TMath::Sqrt(cov[2]), fYfit[1]));
AliInfo(Form("Ref | %7.2f | %7.2f+-%7.2f | %7.2f+-%7.2f| %5.2f | %5.2f |", x, fYref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[0]), fZref[0]-fX*fYref[1], TMath::Sqrt(fRefCov[2]), fYref[1], fZref[1]));
AliInfo(Form("P / Pt [GeV/c] = %f / %f", GetMomentum(), fPt));
if(IsStandAlone()) AliInfo(Form("C Rieman / Vertex [1/cm] = %f / %f", fC[0], fC[1]));
AliInfo(Form("dEdx [a.u.] = %f / %f / %f / %f / %f/ %f / %f / %f", fdEdx[0], fdEdx[1], fdEdx[2], fdEdx[3], fdEdx[4], fdEdx[5], fdEdx[6], fdEdx[7]));
AliInfo(Form("PID = %5.3f / %5.3f / %5.3f / %5.3f / %5.3f", fProb[0], fProb[1], fProb[2], fProb[3], fProb[4]));
if(strcmp(o, "a")!=0) return;
AliTRDcluster* const* jc = &fClusters[0];
for(int ic=0; ic<kNclusters; ic++, jc++) {
if(!(*jc)) continue;
(*jc)->Print(o);
}
}
Bool_t AliTRDseedV1::IsEqual(const TObject *o) const
{
if(!o) return kFALSE;
const AliTRDseedV1 *inTracklet = dynamic_cast<const AliTRDseedV1*>(o);
if(!inTracklet) return kFALSE;
for (Int_t i = 0; i < 2; i++){
if ( fYref[i] != inTracklet->fYref[i] ) return kFALSE;
if ( fZref[i] != inTracklet->fZref[i] ) return kFALSE;
}
if ( TMath::Abs(fS2Y - inTracklet->fS2Y)>1.e-10 ) return kFALSE;
if ( TMath::Abs(GetTilt() - inTracklet->GetTilt())>1.e-10 ) return kFALSE;
if ( TMath::Abs(GetPadLength() - inTracklet->GetPadLength())>1.e-10 ) return kFALSE;
for (Int_t i = 0; i < kNclusters; i++){
if ( fIndexes[i] != inTracklet->fIndexes[i] ) return kFALSE;
}
for (Int_t i=0; i < 2; i++){
if ( fYfit[i] != inTracklet->fYfit[i] ) return kFALSE;
if ( fZfit[i] != inTracklet->fZfit[i] ) return kFALSE;
if ( fLabels[i] != inTracklet->fLabels[i] ) return kFALSE;
}
if ( fN != inTracklet->fN ) return kFALSE;
if ( TMath::Abs(fC[0] - inTracklet->fC[0])>1.e-10 ) return kFALSE;
if ( TMath::Abs(fChi2 - inTracklet->fChi2)>1.e-10 ) return kFALSE;
if ( fDet != inTracklet->fDet ) return kFALSE;
if ( TMath::Abs(fPt - inTracklet->fPt)>1.e-10 ) return kFALSE;
if ( TMath::Abs(fdX - inTracklet->fdX)>1.e-10 ) return kFALSE;
for (Int_t iCluster = 0; iCluster < kNclusters; iCluster++){
AliTRDcluster *curCluster = fClusters[iCluster];
AliTRDcluster *inCluster = inTracklet->fClusters[iCluster];
if (curCluster && inCluster){
if (! curCluster->IsEqual(inCluster) ) {
curCluster->Print();
inCluster->Print();
return kFALSE;
}
} else {
if(curCluster || inCluster) return kFALSE;
}
}
return kTRUE;
}