diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 176671c5fe9a37d6b1b34165c56d2b050d95b12a..4abd25516b3383424982bf4a3876750fae1b86d2 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -645,7 +645,7 @@ InitStatus CbmL1::Init()
if (L1Algo::TrackingMode::kGlobal == fTrackingMode) { //SGtrd2D!!
trdFrontSigma = 1.1;
trdBackSigma = 1.1;
- stationInfo.SetTimeResolution(1.e10);
+ // stationInfo.SetTimeResolution(1.e10);
stationInfo.SetTimeInfo(false);
}
stationInfo.SetFrontBackStripsGeometry(trdFrontPhi, trdFrontSigma, trdBackPhi, trdBackSigma);
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 6246b84b0c6b753cb066222312a896d94093f060..b9e2aaea24054de7567a9261a91be23ee1986c65 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -76,7 +76,7 @@ struct TmpHit {
int iStation; ///< index of station
int ExtIndex; ///< index of hit in the external TClonesArray array (NOTE: negative for MVD)
double u_front, u_back; ///< positions of strips
- double x, y, z; ///< position of hit (Cortesian coordinates)
+ double x, y, z; ///< position of hit (Cartesian coordinates)
double xmc, ymc, p; ///< mc position of hit, total momentum
double tx, ty; ///< mc slopes of the mc point
double dx, dy, dxy; ///< hit position errors in Cortesian coordinates
@@ -108,46 +108,41 @@ struct TmpHit {
/// \param ip
/// \param NStrips
/// \param st reference to the station info object
- // TODO: Probably, L1Station& st parameter should be constant. Do we really want to modify a station here? (S.Zharko)
- void CreateHitFromPoint(const CbmL1MCPoint& point, int det, int nTmpHits, int nStripF, int ip, int& NStrips,
- const L1Station& st)
+ void CreateHitFromPoint(const CbmL1MCPoint& point, int ip, int det, int nTmpHits, int& NStrips, const L1Station& st)
{
ExtIndex = 0;
Det = det;
- id = nTmpHits; //tmpHits.size();
+ id = nTmpHits;
iStation = point.iStation;
dt = st.dt[0];
time = point.time + gRandom->Gaus(0, dt);
- iStripF = nStripF + ip; //firstDetStrip + ip;
+ iStripF = NStrips;
iStripB = iStripF;
- if (NStrips <= iStripF) { NStrips = iStripF + 1; }
+ NStrips++;
+
+ dx = sqrt(st.XYInfo.C00[0]);
+ dy = sqrt(st.XYInfo.C11[0]);
+ dxy = st.XYInfo.C10[0];
+
+ du = sqrt(st.frontInfo.sigma2[0]);
+ dv = sqrt(st.backInfo.sigma2[0]);
- float f_sigma = 0.1; // sqrt(st.frontInfo.sigma2[0]);
- //1.;
- float b_sigma = 0.1; // sqrt(st.backInfo.sigma2[0]);//1.;
+ std::tie(u_front, u_back) = st.ConvXYtoUV(point.x, point.y);
- dx = f_sigma;
- dy = b_sigma;
- dxy = 0;
- du = f_sigma;
- dv = b_sigma;
+ u_front += gRandom->Gaus(0, du);
+ u_back += gRandom->Gaus(0, dv);
- x = point.x; // + gRandom->Gaus(0, dx);
- y = point.y; // + gRandom->Gaus(0, dy);
- z = point.z;
+ std::tie(x, y) = st.ConvUVtoXY(u_front, u_back);
+ z = point.z;
xmc = point.x;
ymc = point.y;
track = point.ID;
-
- p = point.p;
-
- u_front = x * st.frontInfo.cos_phi[0] + y * st.frontInfo.sin_phi[0];
- u_back = x * st.backInfo.cos_phi[0] + y * st.backInfo.sin_phi[0];
- iMC = ip;
+ p = point.p;
+ iMC = ip;
}
/// Sets randomized position and time of the hit
@@ -264,9 +259,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
for (Int_t iMC = 0; iMC < nMvdPointsInEvent; iMC++) {
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 1)) {
- MC.iStation = -1;
+ MC.iStation = -1;
const L1Station* sta = algo->GetParameters()->GetStations().begin();
- double bestDist = 1.e20;
+ double bestDist = 1.e20;
for (Int_t iSt = 0; iSt < NMvdStationsGeom; iSt++) {
// use z_in since z_out is sometimes very wrong
// due to a problem in transport
@@ -303,9 +298,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
for (Int_t iMC = 0; iMC < nMC; iMC++) {
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 0)) {
- MC.iStation = -1;
+ MC.iStation = -1;
const L1Station* sta = algo->GetParameters()->GetStations().begin();
- double bestDist = 1.e20;
+ double bestDist = 1.e20;
for (Int_t iSt = 0; iSt < NStsStationsGeom; iSt++) {
int iStActive = algo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kSts);
if (iStActive == -1) { continue; }
@@ -341,7 +336,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
for (Int_t iMC = 0; iMC < nMC; iMC++) {
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 2)) {
- MC.iStation = -1;
+ MC.iStation = -1;
const L1Station* sta = algo->GetParameters()->GetStations().begin();
for (Int_t iSt = 0; iSt < NMuchStationsGeom; iSt++) {
int iStActive = algo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kMuch);
@@ -369,7 +364,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
for (Int_t iMC = 0; iMC < fTrdPoints->Size(iFile, iEvent); iMC++) {
CbmL1MCPoint MC;
if (!ReadMCPoint(&MC, iMC, iFile, iEvent, 3)) {
- MC.iStation = -1;
+ MC.iStation = -1;
const L1Station* sta = algo->GetParameters()->GetStations().begin();
for (Int_t iSt = 0; iSt < NTrdStationsGeom; iSt++) {
int iStActive = algo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kTrd);
@@ -421,7 +416,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
if (trk_it == dFEI2vMCTracks.end()) continue;
Int_t IND_Track = trk_it->second;
- MC.iStation = -1;
+ MC.iStation = -1;
const L1Station* sta = algo->GetParameters()->GetStations().begin();
for (Int_t iSt = 0; iSt < NTOFStationGeom; iSt++) {
int iStActive = algo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kTof);
@@ -449,7 +444,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
if (!ReadMCPoint(&MC, TofPointToTrack[iTofSta][iMC], iFile, iEvent, 4)) {
- MC.iStation = -1;
+ MC.iStation = -1;
const L1Station* sta = algo->GetParameters()->GetStations().begin();
for (Int_t iSt = 0; iSt < NTOFStation; iSt++) {
int iStActive = algo->GetParameters()->GetStationIndexActive(iSt, L1DetectorID::kTof);
@@ -542,8 +537,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.z = pos.Z();
const L1Station& st = algo->GetParameters()->GetStation(th.iStation);
- th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
+ th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
+ th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
}
th.Det = 0;
th.iMC = -1;
@@ -570,7 +565,6 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
if (fUseSTS && (2 == fStsUseMcHit)) { // create hits from points
- int firstDetStrip = NStrips;
for (int ip = firstStsPoint; ip < firstStsPoint + nStsPoints; ip++) {
const CbmL1MCPoint& p = vMCPoints[ip];
// int mcTrack = p.ID;
@@ -582,8 +576,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
// << p.time << " mc " << p.ID << " p " << p.p << endl;
TmpHit th;
int DetId = 1;
- th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips,
- algo->GetParameters()->GetStation(p.iStation));
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, algo->GetParameters()->GetStation(p.iStation));
tmpHits.push_back(th);
nStsHits++;
}
@@ -661,8 +654,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.dv = mh->GetDv();
const L1Station& st = algo->GetParameters()->GetStation(th.iStation);
- th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
+ th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
+ th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
}
th.iMC = -1;
@@ -731,7 +724,6 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
// Get MuCh hits
//
if ((2 == fMuchUseMcHit) && fUseMUCH) { // create hits from points
- int firstDetStrip = NStrips;
for (int ip = firstMuchPoint; ip < firstMuchPoint + nMuchPoints; ip++) {
const CbmL1MCPoint& p = vMCPoints[ip];
@@ -744,8 +736,7 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
TmpHit th;
int DetId = 2;
- th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips,
- algo->GetParameters()->GetStation(p.iStation));
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, algo->GetParameters()->GetStation(p.iStation));
tmpHits.push_back(th);
nMuchHits++;
@@ -801,8 +792,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
const L1Station& st = algo->GetParameters()->GetStation(th.iStation);
- th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
+ th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
+ th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
}
th.iMC = -1;
int iMC = -1;
@@ -836,151 +827,146 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
nMuchHits++;
} // for j
- } // if listMvdHits
-
- if ((2 == fTrdUseMcHit) && fTrdHitMatches && fUseTRD) { // create hits from points
- int firstDetStrip = NStrips;
-
- for (int ip = firstTrdPoint; ip < firstTrdPoint + nTrdPoints; ip++) {
- const CbmL1MCPoint& p = vMCPoints[ip];
-
- // int mcTrack = p.ID;
- // if (mcTrack < 0) continue;
- // const CbmL1MCTrack& t = vMCTracks[mcTrack];
-
- TmpHit th;
- int DetId = 3;
- th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips,
- algo->GetParameters()->GetStation(p.iStation));
- tmpHits.push_back(th);
- nTrdHits++;
- }
- }
+ } // if listMuchHits
//
// Get TRD hits
//
- if (fUseTRD && listTrdHits && (2 != fTrdUseMcHit)) {
- int firstDetStrip = NStrips;
- vector<bool> mcUsed(nTrdPoints, 0);
- for (int iHit = 0; iHit < listTrdHits->GetEntriesFast(); iHit++) {
- TmpHit th;
+ if (fUseTRD) {
+
+ if (2 == fTrdUseMcHit) { // create hits from MC points
+ for (int ip = firstTrdPoint; ip < firstTrdPoint + nTrdPoints; ip++) {
+ const CbmL1MCPoint& p = vMCPoints[ip];
+ // int mcTrack = p.ID;
+ // if (mcTrack < 0) continue;
+ // const CbmL1MCTrack& t = vMCTracks[mcTrack];
+ TmpHit th;
+ int DetId = 3;
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, algo->GetParameters()->GetStation(p.iStation));
+ tmpHits.push_back(th);
+ nTrdHits++;
+ }
+ }
+ else { // read hits
- CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(listTrdHits->At(iHit));
- // SGtrd2d!! skip TRD-1D hit
- if ((L1Algo::TrackingMode::kGlobal == fTrackingMode) && (int) mh->GetClassType() != 1) { continue; }
+ assert(listTrdHits);
- th.ExtIndex = iHit;
- th.Det = 3;
+ vector<bool> mcUsed(nTrdPoints, 0);
- th.id = tmpHits.size();
+ for (int iHit = 0; iHit < listTrdHits->GetEntriesFast(); iHit++) {
+ TmpHit th;
- int stIdx = algo->GetParameters()->GetStationIndexActive(mh->GetPlaneId(), L1DetectorID::kTrd);
- if (stIdx == -1) continue;
+ CbmTrdHit* mh = L1_DYNAMIC_CAST<CbmTrdHit*>(listTrdHits->At(iHit));
- th.iStation = stIdx;
-
- // if (mh->GetPlaneId()==0) continue;
- // if (mh->GetPlaneId()==1) continue;
- // if (mh->GetPlaneId()==2) continue;
- // if (mh->GetPlaneId()==3) continue;
-
- th.time = mh->GetTime();
- th.dt = mh->GetTimeError();
+ if ((L1Algo::TrackingMode::kGlobal == fTrackingMode) && (int) mh->GetClassType() != 1) {
+ // SGtrd2d!! skip TRD-1D hit
+ continue;
+ }
- // th.iSector = 0;
- th.iStripF = firstDetStrip + iHit;
- th.iStripB = th.iStripF; //TrdHitsOnStationIndex[num+1][k];
- if (NStrips <= th.iStripF) { NStrips = th.iStripF + 1; }
+ th.ExtIndex = iHit;
+ th.Det = 3;
+ th.id = tmpHits.size();
- TVector3 pos, err;
- mh->Position(pos);
- mh->PositionError(err);
+ int stIdx = algo->GetParameters()->GetStationIndexActive(mh->GetPlaneId(), L1DetectorID::kTrd);
+ if (stIdx == -1) continue;
- th.x = pos.X();
- th.y = pos.Y();
- th.z = pos.Z();
+ th.iStation = stIdx;
+ // if (mh->GetPlaneId()==0) continue;
+ // if (mh->GetPlaneId()==1) continue;
+ // if (mh->GetPlaneId()==2) continue;
+ // if (mh->GetPlaneId()==3) continue;
- th.dx = fabs(mh->GetDx());
- th.dy = fabs(mh->GetDy());
- th.dxy = 0;
+ th.time = mh->GetTime();
+ th.dt = mh->GetTimeError();
- th.du = fabs(mh->GetDx());
- th.dv = fabs(mh->GetDy());
+ // th.iSector = 0;
+ th.iStripF = NStrips;
+ th.iStripB = th.iStripF; //TrdHitsOnStationIndex[num+1][k];
+ NStrips++;
- const L1Station& st = algo->GetParameters()->GetStation(th.iStation);
- th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
+ TVector3 pos, err;
+ mh->Position(pos);
+ mh->PositionError(err);
- th.iMC = -1;
+ th.x = pos.X();
+ th.y = pos.Y();
+ th.z = pos.Z();
- if (fPerformance && fTrdHitMatches) {
+ th.dx = fabs(mh->GetDx());
+ th.dy = fabs(mh->GetDy());
+ th.dxy = 0;
- CbmMatch* trdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(iHit));
+ th.du = fabs(mh->GetDx());
+ th.dv = fabs(mh->GetDy());
- if (trdHitMatch->GetNofLinks() > 0) {
+ const L1Station& st = algo->GetParameters()->GetStation(th.iStation);
- int iMC = trdHitMatch->GetLink(0).GetIndex();
+ std::tie(th.u_front, th.u_back) = st.ConvXYtoUV(th.x, th.y);
+
+ th.iMC = -1;
+ th.track = -1;
+ int iMcTrd = -1;
+ if (fPerformance && fTrdHitMatches) {
+ CbmMatch* trdHitMatch = L1_DYNAMIC_CAST<CbmMatch*>(fTrdHitMatches->At(iHit));
+ if (trdHitMatch->GetNofLinks() > 0) {
+ iMcTrd = trdHitMatch->GetLink(0).GetIndex();
+ assert(iMcTrd >= 0 && iMcTrd < nTrdPoints);
+ th.iMC = iMcTrd + nMvdPoints + nStsPoints + nMuchPoints;
+ th.track = vMCPoints[th.iMC].ID;
+ }
+ }
- assert(iMC >= 0 && iMC < nTrdPoints);
+ if (1 == fTrdUseMcHit) { // replace hit by MC points
- th.iMC = iMC + nMvdPoints + nStsPoints + nMuchPoints;
- th.track = vMCPoints[th.iMC].ID;
+ assert(fPerformance && fTrdHitMatches);
- if (1 == fTrdUseMcHit) { // replace hits by MC points
+ if (th.iMC < 0) continue; // skip noise hits
+ if (mcUsed[iMcTrd]) continue; // one hit per MC point
- if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
- // SGtrd2d!! don't smear
- th.dx = 0.;
- th.dy = 0.;
- th.dt = 0.;
- }
+ if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
+ // SGtrd2d!! don't smear
+ TmpHit thsave = th;
+ th.dx = 0.;
+ th.dy = 0.;
+ th.dt = 0.;
th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetParameters()->GetStation(th.iStation));
- /*
- CbmTrdPoint* pt =
- (CbmTrdPoint*) fTrdPoints->Get(trdHitMatch->GetLink(0).GetFile(), trdHitMatch->GetLink(0).GetEntry(),
- trdHitMatch->GetLink(0).GetIndex());
-
- cout << " dx " << th.x - 0.5 * (pt->GetXOut() + pt->GetXIn()) << " dy "
- << th.y - 0.5 * (pt->GetYOut() + pt->GetYIn()) << " dz "
- << th.z - 0.5 * (pt->GetZOut() + pt->GetZIn()) << " z " << th.z << endl;
- */
- if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
- // SGtrd2d!! enlarge errors
- th.dx = 1.1;
- th.du = 1.1;
- th.dy = 1.1;
- th.dv = 1.1;
- th.dt = 10000.;
- }
- if (th.iMC < 0) continue;
- if (mcUsed[iMC]) continue;
- if (0) { // select specific tracks only
- int mcTrack = -1;
- float mcP = 0;
- if (th.iMC >= 0) {
- mcTrack = vMCPoints[th.iMC].ID;
- if (mcTrack >= 0) { mcP = vMCTracks[mcTrack].p; }
- }
- if (mcP < 1.) continue;
- }
- mcUsed[iMC] = 1;
+ th.dx = thsave.dx;
+ th.dy = thsave.dy;
+ th.dt = thsave.dt;
+ // SGtrd2d!! artificial errors
+ th.dx = 1.1;
+ th.du = 1.1;
+ th.dy = 1.1;
+ th.dv = 1.1;
+ }
+ else {
+ th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetParameters()->GetStation(th.iStation));
+ }
+ mcUsed[iMcTrd] = 1;
+ } // replace hit by MC points
+
+ if (0) { // select specific tracks
+ int mcTrack = -1;
+ float mcP = 0;
+ if (th.iMC >= 0) {
+ mcTrack = vMCPoints[th.iMC].ID;
+ if (mcTrack >= 0) { mcP = vMCTracks[mcTrack].p; }
}
+ if (mcP < 1.) continue;
}
- }
- tmpHits.push_back(th);
- nTrdHits++;
- } // for listTrdHits
- } // if listTrdHits
+
+ tmpHits.push_back(th);
+ nTrdHits++;
+ } // for listTrdHits
+ }
+ } // read TRD hits
//
// Get ToF hits
//
if ((2 == fTofUseMcHit) && fUseTOF) { // create hits from points
- int firstDetStrip = NStrips;
-
for (int ip = firstTofPoint; ip < firstTofPoint + nTofPoints; ip++) {
const CbmL1MCPoint& p = vMCPoints[ip];
@@ -991,12 +977,9 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
//if (t.Points.size() > 4) continue;
TmpHit th;
-
int DetId = 4;
- th.CreateHitFromPoint(p, DetId, tmpHits.size(), firstDetStrip, ip, NStrips,
- algo->GetParameters()->GetStation(p.iStation));
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, algo->GetParameters()->GetStation(p.iStation));
tmpHits.push_back(th);
-
nTofHits++;
}
}
@@ -1065,8 +1048,8 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
th.iStation = stIdx;
const L1Station& st = algo->GetParameters()->GetStation(th.iStation);
- th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
- th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
+ th.u_front = th.x * st.frontInfo.cos_phi[0] + th.y * st.frontInfo.sin_phi[0];
+ th.u_back = th.x * st.backInfo.cos_phi[0] + th.y * st.backInfo.sin_phi[0];
th.iMC = -1;
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index cbc8feb44ddaa9e184fb9b478685969afb75b1ee..44667d32ec2110a4b15407703036901b1b26d0dc 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -98,7 +98,6 @@ void L1Algo::Init(const bool UseHitErrors, const TrackingMode mode, const bool M
fTrackingLevel = fInitManager.GetTrackingLevel();
fGhostSuppression = fInitManager.GetGhostSuppression();
fMomentumCutOff = fInitManager.GetMomentumCutOff();
-
}
@@ -168,10 +167,10 @@ void L1Algo::SetData(L1Vector<L1Hit>& Hits_, int nStrips_, L1Vector<unsigned cha
void L1Algo::GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, char iS)
{
const L1Station& sta = fParameters.GetStation(int(iS));
- fscal u = _h.u;
- fscal v = _h.v;
- _x = (sta.xInfo.sin_phi[0] * u + sta.xInfo.cos_phi[0] * v) / (_h.z - fParameters.GetTargetPositionZ()[0]);
- _y = (sta.yInfo.cos_phi[0] * u + sta.yInfo.sin_phi[0] * v) / (_h.z - fParameters.GetTargetPositionZ()[0]);
+ std::tie(_x, _y) = sta.ConvUVtoXY(_h.u, _h.v);
+ float dz = _h.z - fParameters.GetTargetPositionZ()[0];
+ _x /= dz;
+ _y /= dz;
}
void L1Algo::GetHitCoor(const L1Hit& _h, fscal& _x, fscal& _y, fscal& _z, const L1Station& sta)
@@ -189,13 +188,10 @@ void L1Algo::StripsToCoor(
const fscal& u, const fscal& v, fscal& _x, fscal& _y,
const L1Station& sta) const // TODO: Actually sta.yInfo.sin_phi is same for all stations, so ...
{
- // fvec x,y;
- // StripsToCoor(u,v,x,y,sta);
- // _x = x[0];
- // _y = y[0];
- _x = sta.xInfo.sin_phi[0] * u + sta.xInfo.cos_phi[0] * v;
- _y = sta.yInfo.cos_phi[0] * u + sta.yInfo.sin_phi[0] * v;
+ std::tie(_x, _y) = sta.ConvUVtoXY(u, v);
}
+
+
/// convert strip positions to coordinates
void L1Algo::StripsToCoor(
const fscal& u, const fscal& v, fvec& _x, fvec& _y,
@@ -205,7 +201,7 @@ void L1Algo::StripsToCoor(
// StripsToCoor(u,v,x,y,sta);
// _x = x[0];
// _y = y[0];
- _x = sta.xInfo.sin_phi * u + sta.xInfo.cos_phi * v;
+ _x = sta.xInfo.cos_phi * u + sta.xInfo.sin_phi * v;
_y = sta.yInfo.cos_phi * u + sta.yInfo.sin_phi * v;
}
@@ -216,12 +212,12 @@ void L1Algo::dUdV_to_dY(const fvec& u, const fvec& v, fvec& _y, const L1Station&
void L1Algo::dUdV_to_dX(const fvec& u, const fvec& v, fvec& _x, const L1Station& sta)
{
- _x = sqrt((sta.xInfo.sin_phi * u) * (sta.xInfo.sin_phi * u) + (sta.xInfo.cos_phi * v) * (sta.xInfo.cos_phi * v));
+ _x = sqrt((sta.xInfo.cos_phi * u) * (sta.xInfo.cos_phi * u) + (sta.xInfo.sin_phi * v) * (sta.xInfo.sin_phi * v));
}
void L1Algo::dUdV_to_dXdY(const fvec& u, const fvec& v, fvec& _xy, const L1Station& sta)
{
- _xy = ((sta.xInfo.sin_phi * u) * (sta.yInfo.cos_phi * u) + (sta.xInfo.cos_phi * v) * (sta.yInfo.sin_phi * v));
+ _xy = ((sta.xInfo.cos_phi * u) * (sta.yInfo.cos_phi * u) + (sta.xInfo.sin_phi * v) * (sta.yInfo.sin_phi * v));
}
void L1Algo::StripsToCoor(
@@ -231,7 +227,7 @@ void L1Algo::StripsToCoor(
// only for same-z
// x = u;
// y = (sta.yInfo.cos_phi*u + sta.yInfo.sin_phi*v);
- x = sta.xInfo.sin_phi * u + sta.xInfo.cos_phi * v;
+ x = sta.xInfo.cos_phi * u + sta.xInfo.sin_phi * v;
y = sta.yInfo.cos_phi * u + sta.yInfo.sin_phi * v;
}
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx
index 8b70495d31ea239177661be2b8cbc07611af6726..a978aaaa8c47bd0e8dba87b35df150780714dda1 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.cxx
+++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx
@@ -309,18 +309,21 @@ void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double front
fL1Station.backInfo.sin_phi = sBack;
fL1Station.backInfo.sigma2 = backSigma * backSigma;
- double det = (fabs(backPhi - frontPhi) < 0.1) ? cos(backPhi - frontPhi) : (cFront * sBack - sFront * cBack);
- det *= det;
- fL1Station.XYInfo.C00 = (sBack * sBack * frontSigma * frontSigma + sFront * sFront * backSigma * backSigma) / det;
- fL1Station.XYInfo.C10 = -(sBack * cBack * frontSigma * frontSigma + sFront * cFront * backSigma * backSigma) / det;
- fL1Station.XYInfo.C11 = (cBack * cBack * frontSigma * frontSigma + cFront * cFront * backSigma * backSigma) / det;
-
- fL1Station.xInfo.cos_phi = -sFront / (cFront * sBack - cBack * sFront);
- fL1Station.xInfo.sin_phi = sBack / (cFront * sBack - cBack * sFront);
+ double det = cFront * sBack - sFront * cBack;
+
+ assert(fabs(det) > 1.e-2);
+ double det2 = det * det;
+
+ fL1Station.XYInfo.C00 = (sBack * sBack * frontSigma * frontSigma + sFront * sFront * backSigma * backSigma) / det2;
+ fL1Station.XYInfo.C10 = -(sBack * cBack * frontSigma * frontSigma + sFront * cFront * backSigma * backSigma) / det2;
+ fL1Station.XYInfo.C11 = (cBack * cBack * frontSigma * frontSigma + cFront * cFront * backSigma * backSigma) / det2;
+
+ fL1Station.xInfo.cos_phi = sBack / det;
+ fL1Station.xInfo.sin_phi = -sFront / det;
fL1Station.xInfo.sigma2 = fL1Station.XYInfo.C00;
- fL1Station.yInfo.cos_phi = cBack / (cBack * sFront - cFront * sBack);
- fL1Station.yInfo.sin_phi = -cFront / (cBack * sFront - cFront * sBack);
+ fL1Station.yInfo.cos_phi = -cBack / det;
+ fL1Station.yInfo.sin_phi = cFront / det;
fL1Station.yInfo.sigma2 = fL1Station.XYInfo.C11;
//-------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index 271a9cd96d4d8ea611ef2d64ce2eec070e8357e8..8732fa5c4e48e291a87c73510eff1e120ef4e47b 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -232,16 +232,16 @@ private:
/// Forms parameters container
void FormParametersContainer();
- /* Basic fields */
+ // * Basic fields *
InitController_t fInitController {}; ///< Initialization flags
L1DetectorIDSet_t fActiveDetectorIDs {}; ///< Set of tracking detectors, active during this analysis session
- /* Target fields */
+ // * Target fields *
double fTargetZ {0.}; ///< Target position z component in double precision
- /* Stations related fields */
+ // * Stations related fields *
std::set<L1BaseStationInfo> fStationsInfo {}; ///< Set of L1BaseStationInfo objects
@@ -257,7 +257,7 @@ private:
L1FieldFunction_t fFieldFunction {[](const double (&)[3], double (&)[3]) {}};
// NOTE: Stations of daetectors which will not be assigned as active, will not be included in the tracking!
- /* CA track finder iterations related */
+ // * CA track finder iterations related *
int fCAIterationsNumberCrosscheck {-1}; ///< Number of iterations to be passed (must be used for cross-checks)
diff --git a/reco/L1/L1Algo/L1Station.cxx b/reco/L1/L1Algo/L1Station.cxx
index 18b30f67eb963bef461c2efaf9493466a07e6b47..5aa0c4c5446826232796b9bab40d6f8f0e3b7897 100644
--- a/reco/L1/L1Algo/L1Station.cxx
+++ b/reco/L1/L1Algo/L1Station.cxx
@@ -83,6 +83,24 @@ void L1Station::CheckConsistency() const
xInfo.CheckConsistency();
yInfo.CheckConsistency();
XYInfo.CheckConsistency();
+
+ // Check consistency of coordinate transformations
+ for (float x = -1.f; x < 1.1f; x += 0.2) {
+ for (float y = -1.f; y < 1.1f; y += 0.2) {
+ float u, v;
+ std::tie(u, v) = ConvXYtoUV(x, y);
+ float x1, y1;
+ std::tie(x1, y1) = ConvUVtoXY(u, v);
+
+ if (fabs(x1 - x) > 1.e-6 || fabs(y1 - y) > 1.e-6) {
+ std::stringstream msg;
+ msg << "L1Station: " << this->ToString() << " makes wrong coordinate convertion: "
+ << "x,y " << x << "," << y << " -> u,v " << u << "," << v << " -> x,y " << x1 << "," << y1
+ << " are not the same: dx,dy " << x1 - x << "," << y1 - y;
+ throw std::logic_error(msg.str());
+ }
+ }
+ }
}
// TODO: Improve log style (S.Zh.)
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index e21d21cfd97d6ade96e12bfb14f596e538a32b5c..01b82ae659273771dbe69b5bfdf491519d6f1f14 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -29,8 +29,8 @@ struct L1Station {
L1FieldSlice fieldSlice {};
L1UMeasurementInfo frontInfo {};
L1UMeasurementInfo backInfo {};
- L1UMeasurementInfo xInfo {};
- L1UMeasurementInfo yInfo {};
+ L1UMeasurementInfo xInfo {}; ///< x axis in front,back coordinates
+ L1UMeasurementInfo yInfo {}; ///< y axis in front,back coordinates
L1XYMeasurementInfo XYInfo {};
/// Prints object fields
@@ -47,6 +47,19 @@ struct L1Station {
/// initialization
void CheckConsistency() const;
+ /// convert x,y to u,v
+ std::pair<float, float> ConvXYtoUV(float x, float y) const
+ {
+ return std::make_pair(x * frontInfo.cos_phi[0] + y * frontInfo.sin_phi[0],
+ x * backInfo.cos_phi[0] + y * backInfo.sin_phi[0]);
+ }
+
+ /// convert u,v to x,y
+ std::pair<float, float> ConvUVtoXY(float u, float v) const
+ {
+ return std::make_pair(u * xInfo.cos_phi[0] + v * xInfo.sin_phi[0], u * yInfo.cos_phi[0] + v * yInfo.sin_phi[0]);
+ }
+
} _fvecalignment;
#endif // L1Station_h