9#include <tracking/trackFindingCDC/sim/CDCSimpleSimulation.h>
11#include <tracking/trackFindingCDC/eventdata/tracks/CDCTrack.h>
13#include <tracking/trackFindingCDC/eventdata/hits/CDCRecoHit3D.h>
14#include <tracking/trackFindingCDC/eventdata/hits/CDCRLWireHit.h>
15#include <tracking/trackFindingCDC/eventdata/hits/CDCWireHit.h>
17#include <tracking/trackFindingCDC/topology/CDCWire.h>
18#include <tracking/trackFindingCDC/topology/CDCWireLayer.h>
19#include <tracking/trackFindingCDC/topology/CDCWireTopology.h>
21#include <tracking/trackFindingCDC/utilities/VectorRange.h>
23#include <framework/gearbox/Const.h>
24#include <framework/logging/Logger.h>
31using namespace TrackFindingCDC;
44 return std::move(
simulate(std::vector<CDCTrajectory3D>(1, trajectory3D)).front());
50 std::vector<SimpleSimHit> simpleSimHits;
51 const size_t nMCTracks = trajectories3D.size();
53 for (
size_t iMCTrack = 0; iMCTrack < nMCTracks; ++iMCTrack) {
59 Helix globalHelix = localHelix;
60 const double arcLength2DOffset = globalHelix.
passiveMoveBy(-localOrigin);
61 std::vector<SimpleSimHit> simpleSimHitsForTrajectory =
createHits(globalHelix, arcLength2DOffset);
63 for (
SimpleSimHit& simpleSimHit : simpleSimHitsForTrajectory) {
64 simpleSimHit.m_iMCTrack = iMCTrack;
65 simpleSimHits.push_back(simpleSimHit);
69 std::vector<CDCTrack> mcTracks =
constructMCTracks(nMCTracks, std::move(simpleSimHits));
72 for (
size_t iMCTrack = 0; iMCTrack < nMCTracks; ++iMCTrack) {
73 CDCTrack& mcTrack = mcTracks[iMCTrack];
75 if (not mcTrack.empty()) {
93 std::stable_sort(simpleSimHits.begin(), simpleSimHits.end(),
95 return lhs.m_wireHit < rhs.m_wireHit;
100 const CDCWire* lastWire =
nullptr;
101 size_t nSameWire = 0;
104 auto exceedsMaxNHitOnWire =
105 [&lastWire, &nSameWire, maxNHitOnWire](
const SimpleSimHit & simpleSimHit) ->
bool {
107 if (&(simpleSimHit.m_wireHit.getWire()) == lastWire)
113 lastWire = &(simpleSimHit.m_wireHit.getWire());
115 return nSameWire > maxNHitOnWire ? true :
false;
118 auto itLast = std::remove_if(simpleSimHits.begin(), simpleSimHits.end(), exceedsMaxNHitOnWire);
119 simpleSimHits.erase(itLast, simpleSimHits.end());
124 std::vector<CDCWireHit> wireHits;
125 wireHits.reserve(simpleSimHits.size());
127 wireHits.push_back(simpleSimHit.m_wireHit);
130 B2ASSERT(
"WireHits should be sorted as a result from sorting the SimpleSimHits. "
131 "Algorithms may relay on the sorting o the WireHits",
132 std::is_sorted(wireHits.begin(), wireHits.end()));
134 m_sharedWireHits.reset(
new const std::vector<CDCWireHit>(std::move(wireHits)));
140 std::vector<CDCTrack> mcTracks;
141 mcTracks.resize(nMCTracks);
143 const size_t nWireHits = wireHits.
size();
145 for (
size_t iWireHit = 0; iWireHit < nWireHits; ++iWireHit) {
146 const CDCWireHit& wireHit = wireHits[iWireHit];
147 const SimpleSimHit& simpleSimHit = simpleSimHits[iWireHit];
153 mcTrack.push_back(recoHit3D);
157 for (
CDCTrack& mcTrack : mcTracks) {
166std::vector<CDCSimpleSimulation::SimpleSimHit>
168 double arcLength2DOffset)
const
171 std::vector<SimpleSimHit> simpleSimHits;
180 const double localArcLength2DToOuterWall = arcLength2DOffset + globalArcLength2DToOuterWall;
182 if (localArcLength2DToOuterWall < 0) {
185 B2WARNING(
"Simple simulation got trajectory outside CDC that moves away from the detector.");
186 return simpleSimHits;
192 const bool isCurler = std::isnan(localArcLength2DToOuterWall);
193 const double perimeterXY = globalHelix.
perimeterXY();
194 const double maxArcLength2D = isCurler ? fabs(perimeterXY) : localArcLength2DToOuterWall;
197 B2INFO(
"Simulating curler");
201 double outerR = wireLayer.getOuterCylindricalR();
202 double innerR = wireLayer.getInnerCylindricalR();
204 if ((maxR < innerR) or (outerR < minR)) {
209 double centerR = (std::min(outerR, maxR) + std::max(innerR, minR)) / 2;
212 double localArcLength2D = arcLength2DOffset + globalArcLength2D;
215 std::vector<SimpleSimHit> simpleSimHitsInLayer;
216 if (localArcLength2D > 0 and localArcLength2D < maxArcLength2D) {
219 const CDCWire& closestWire = wireLayer.getClosestWire(pos3DAtLayer);
221 simpleSimHitsInLayer =
createHitsForLayer(closestWire, globalHelix, arcLength2DOffset);
223 for (
SimpleSimHit& simpleSimHit : simpleSimHitsInLayer) {
224 if (simpleSimHit.m_arcLength2D < maxArcLength2D) {
225 simpleSimHits.push_back(simpleSimHit);
229 B2INFO(
"Arc length to long");
232 bool oneSegment = outerR > maxR or innerR < minR;
233 if (not oneSegment) {
236 double secondGlobalArcLength2D = -globalArcLength2D;
237 double secondArcLength2DOffset = arcLength2DOffset;
238 double secondLocalArcLength2D = secondArcLength2DOffset + secondGlobalArcLength2D;
240 if (isCurler and secondLocalArcLength2D < 0) {
241 secondLocalArcLength2D += perimeterXY;
242 secondArcLength2DOffset += perimeterXY;
243 secondGlobalArcLength2D += perimeterXY;
246 if (secondLocalArcLength2D > 0 and secondLocalArcLength2D < maxArcLength2D) {
248 const CDCWire& closestWire = wireLayer.getClosestWire(pos3DAtLayer);
251 bool wireAlreadyHit =
false;
252 for (
const SimpleSimHit& simpleSimHit : simpleSimHits) {
253 if (simpleSimHit.m_wireHit.isOnWire(closestWire)) {
254 wireAlreadyHit =
true;
257 if (not wireAlreadyHit) {
258 std::vector<SimpleSimHit> secondSimpleSimHitsInLayer =
261 for (
SimpleSimHit& simpleSimHit : secondSimpleSimHitsInLayer) {
262 if (simpleSimHit.m_arcLength2D < maxArcLength2D) {
263 simpleSimHits.push_back(simpleSimHit);
271 return simpleSimHits;
274std::vector<CDCSimpleSimulation::SimpleSimHit>
276 const Helix& globalHelix,
277 double arcLength2DOffset)
const
279 std::vector<SimpleSimHit> result;
283 result.push_back(simpleSimHit);
294 result.push_back(simpleSimHitForWire);
306 result.push_back(simpleSimHitForWire);
316 const Helix& globalHelix,
317 double arcLength2DOffset)
const
320 if ((arcLength2D + arcLength2DOffset) < 0) {
328 double correctedArcLength2D = arcLength2D;
330 for (
int c_Iter = 0; c_Iter < 2; c_Iter++) {
335 if ((correctedArcLength2D + arcLength2DOffset) < 0) {
338 correctedPos3D = globalHelix.
atArcLength2D(correctedArcLength2D);
341 const double trueDriftLength = wire.
getDriftLength(correctedPos3D);
342 const double smearedDriftLength = trueDriftLength + gRandom->Gaus(0,
m_driftLengthSigma);
346 double arcLength3D = hypot2(1, globalHelix.
tanLambda()) * (correctedArcLength2D + arcLength2DOffset);
354 double distanceToBack = (wirePos.
z() - backwardZ) * hypot2(1, wire.
getTanStereoAngle());
359 double measuredDriftLength = smearedDriftLength + delayTime *
m_driftSpeed;
376 correctedArcLength2D,
385 const size_t nMCTracks = 2;
386 std::vector<SimpleSimHit> simpleSimHits;
387 simpleSimHits.reserve(128 + 64);
611 std::vector<CDCTrack> mcTracks =
constructMCTracks(nMCTracks, std::move(simpleSimHits));
static const double speedOfLight
[cm/ns]
Class representing an oriented hit wire including a hypotheses whether the causing track passes left ...
Class representing a three dimensional reconstructed hit.
ConstVectorRange< CDCWireHit > getWireHits() const
Getter for the wire hits created in the simulation.
std::vector< CDCTrack > constructMCTracks(int nMCTracks, std::vector< SimpleSimHit > simpleSimHits)
Creates CDCWireHits and uses them to construct the true CDCTracks.
std::vector< CDCTrack > simulate(const std::vector< CDCTrajectory3D > &trajectories3D)
Propagates the trajectories through the CDC as without energy loss until they first leave the CDC.
bool m_addInWireSignalDelay
Switch to activate the in wire signal delay.
double m_propSpeed
Electrical current propagation speed in the wires.
std::shared_ptr< const std::vector< CDCWireHit > > m_sharedWireHits
Space for the memory of the generated wire hits.
double getEventTime() const
Getter for a global event time offset.
double m_driftLengthVariance
Variance by which the drift length should be smeared.
std::vector< SimpleSimHit > createHitsForLayer(const CDCWire &nearWire, const Helix &globalHelix, double arcLength2DOffset) const
Generate connected hits for wires in the same layer close to the given wire.
bool m_addTOFDelay
Switch to activate the addition of the time of flight.
SimpleSimHit createHitForCell(const CDCWire &wire, const Helix &globalHelix, double arcLength2DOffset) const
Generate a hit for the given wire.
int m_maxNHitOnWire
Maximal number of hits allowed on each wire (0 means all).
double m_driftLengthSigma
Standard deviation by which the drift length should be smeared.
std::vector< CDCTrack > loadPreparedEvent()
Fills the wire hits with a hard coded event from the real simulation.
double m_driftSpeed
Electron drift speed in the cdc gas.
std::vector< SimpleSimHit > createHits(const Helix &globalHelix, double arcLength2DOffset) const
Generate hits for the given helix in starting from the two dimensional arc length.
Class representing a sequence of three dimensional reconstructed hits.
void setStartTrajectory3D(const CDCTrajectory3D &startTrajectory3D)
Setter for the two dimensional trajectory.
void sortByArcLength2D()
Sort the recoHits according to their perpS information.
void setEndTrajectory3D(const CDCTrajectory3D &endTrajectory3D)
Setter for the three dimensional trajectory.
Particle full three dimensional trajectory.
const UncertainHelix & getLocalHelix() const
Getter for the helix in local coordinates.
double setLocalOrigin(const Vector3D &localOrigin)
Setter for the origin of the local coordinate system.
const Vector3D & getLocalOrigin() const
Getter for the origin of the local coordinate system.
Class representing a hit wire in the central drift chamber.
double getRefDriftLength() const
Getter for the drift length at the reference position of the wire.
Class representing a sense wire layer in the central drift chamber.
Class representing the sense wire arrangement in the whole of the central drift chamber.
const std::vector< Belle2::TrackFindingCDC::CDCWireLayer > & getWireLayers() const
Getter for the underlying storing layer vector.
static CDCWireTopology & getInstance()
Getter for the singleton instance of the wire topology.
double getOuterCylindricalR() const
Getter for the outer radius of the outer most wire layer.
Class representing a sense wire in the central drift chamber.
double getDriftLength(const Vector3D &pos3D) const
Calculates the straight drift length from the position to the wire This is essentially the same as th...
const WireID & getWireID() const
Getter for the wire id.
Vector2D getWirePos2DAtZ(const double z) const
Gives the xy projected position of the wire at the given z coordinate.
Vector3D getClosest(const Vector3D &pos3D) const
Calculates the closest approach in the wire to the position.
MayBePtr< const CDCWire > getNeighborCCW() const
Gives the closest neighbor in the counterclockwise direction - always exists.
MayBePtr< const CDCWire > getNeighborCW() const
Gives the closest neighbor in the clockwise direction - always exists.
const Vector2D & getRefPos2D() const
Getter for the wire reference position for 2D tracking Gives the wire's reference position projected ...
double getTanStereoAngle() const
Getter for the tangents of the stereo angle of the wire.
double getBackwardZ() const
Getter for the z coordinate at the backward joint points of the wires.
Extension of the generalized circle also caching the perigee coordinates.
double minimalCylindricalR() const
Gives the minimal cylindrical radius the circle reaches (unsigned)
double perimeterXY() const
Getter for the perimeter of the circle in the xy projection.
double arcLength2DToXY(const Vector2D &point) const
Calculates the two dimensional arc length that is closest to two dimensional point in the xy projecti...
double tanLambda() const
Getter for the proportinality factor from arc length in xy space to z.
double maximalCylindricalR() const
Gives the maximal cylindrical radius the circle reaches.
double arcLength2DToCylindricalR(double cylindricalR) const
Calculates the two dimensional arc length that first reaches a cylindrical radius on the helix Return...
const PerigeeCircle & circleXY() const
Getter for the projection into xy space.
Vector3D atArcLength2D(double s) const
Calculates the point, which lies at the give perpendicular travel distance (counted from the perigee)
double passiveMoveBy(const Vector3D &by)
Moves the coordinates system by the given vector.
ERightLeft isRightOrLeft(const Vector2D &point) const
Indicates if the point is on the right or left side of the circle.
A pair of iterators usable with the range base for loop.
std::size_t size() const
Returns the total number of objects in this range.
A general helix class including a covariance matrix.
A two dimensional vector which is equipped with functions for correct handling of orientation relate...
A three dimensional vector.
double z() const
Getter for the z coordinate.
Class to identify a wire inside the CDC.
ERightLeft
Enumeration to represent the distinct possibilities of the right left passage.
Abstract base class for different kinds of events.
Structure to accommodate information about the individual hits during the simulation.
Vector3D m_pos3D
Memory for the true position on the track closest to the wire.
double m_arcLength2D
Memory for the true two dimensional arc length on the helix to this hit.
CDCWireHit m_wireHit
Memory for the wire hit instance that will be given to the reconstruction.
ERightLeft m_rlInfo
Memory for the true right left passage information.
double m_trueDriftLength
Memory for the true drift length from the true position to the wire.
size_t m_iMCTrack
Memory for the true index of the track this hit is contained in.