Belle II Software  release-05-01-25
CDCCKFEclSeedCreator.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Simon Kurz, Nils Braun *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <tracking/ckf/cdc/findlets/CDCCKFEclSeedCreator.h>
11 
12 #include <tracking/trackFindingCDC/utilities/StringManipulation.h>
13 
14 #include <tracking/dataobjects/RecoTrack.h>
15 #include <tracking/ckf/cdc/entities/CDCCKFState.h>
16 
17 #include <framework/core/ModuleParamList.h>
18 
19 #include <framework/gearbox/Const.h>
20 
21 using namespace Belle2;
22 
23 CDCCKFEclSeedCreator::CDCCKFEclSeedCreator() : Super()
24 {
25 
26 }
27 
28 void CDCCKFEclSeedCreator::exposeParameters(ModuleParamList* moduleParamList, const std::string& prefix)
29 {
30  moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "inputECLshowersStoreArrayName"),
31  m_param_inputEclShowerStoreArrayName,
32  "StoreArray name of the input Shower Store Array.");
33 
34  moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "eclSeedRecoTrackStoreArrayName"),
35  m_param_eclSeedRecoTrackStoreArrayName,
36  "StoreArray name of the output Track Store Array.");
37 
38  moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "minimalEnRequirementCluster"),
39  m_param_minimalEnRequirement,
40  "Minimal energy requirement for the input clusters",
41  m_param_minimalEnRequirement);
42 
43  moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "restrictToForwardSeeds"),
44  m_param_restrictToForwardSeeds,
45  "Don't do Ecl seeding in central region to save computing time",
46  m_param_restrictToForwardSeeds);
47 
48  moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "tanLambdaForwardNeg"),
49  m_param_tanLambdaForwardNeg,
50  "Up to which (neg) tanLambda value should the seeding be performed",
51  m_param_tanLambdaForwardNeg);
52 
53  moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "tanLambdaForwardPos"),
54  m_param_tanLambdaForwardPos,
55  "Up to which (pos) tanLambda value should the seeding be performed",
56  m_param_tanLambdaForwardPos);
57 
58  moduleParamList->addParameter(TrackFindingCDC::prefixed(prefix, "showerDepth"),
59  m_param_showerDepth,
60  "Don't do Ecl seeding in central region to save computing time",
61  m_param_showerDepth);
62 }
63 
64 void CDCCKFEclSeedCreator::initialize()
65 {
66  Super::initialize();
67 
68  m_inputECLshowers.isRequired(m_param_inputEclShowerStoreArrayName);
69 
70  m_eclSeedRecoTracks.registerInDataStore(m_param_eclSeedRecoTrackStoreArrayName);
71  RecoTrack::registerRequiredRelations(m_eclSeedRecoTracks);
72 
73  m_eclSeedRecoTracks.registerRelationTo(m_inputECLshowers);
74 }
75 
76 void CDCCKFEclSeedCreator::apply(std::vector<CDCCKFPath>& seeds)
77 {
78  // loop over all showers and create seed objects
79  for (auto& shower : m_inputECLshowers) {
80  if (shower.getHypothesisId() != ECLShower::c_nPhotons) {
81  continue;
82  }
83 
84  const double Eclus = shower.getEnergy();
85  if (Eclus < m_param_minimalEnRequirement) {
86  continue;
87  }
88 
89  // Get shower position, momentum
90  const double thetaClus = shower.getTheta();
91  const double phiClus = shower.getPhi();
92  const double rClus = shower.getR();
93 
94  const double sinTheta = sin(thetaClus);
95  const double cosTheta = cos(thetaClus);
96  const double sinPhi = sin(phiClus);
97  const double cosPhi = cos(phiClus);
98 
99  TVector3 pos(rClus * sinTheta * cosPhi, rClus * sinTheta * sinPhi, rClus * cosTheta);
100  const double tanLambda = pos.Z() / pos.Perp();
101 
102  // restrict to forward seeds
103  if (m_param_restrictToForwardSeeds) {
104  if (tanLambda > m_param_tanLambdaForwardNeg && tanLambda < m_param_tanLambdaForwardPos) {
105  continue;
106  }
107  }
108 
109  // Correction if shower is assumed to start in a certain depth
110  pos = pos - m_param_showerDepth / pos.Mag() * pos;
111 
112  TVector3 mom = Eclus / pos.Mag() * pos;
113 
114  // Calculate helix trajectory for negative and positive charge seeds
115  // Find center of circular trajectory
116  // factor 1. for charge // speed of light in [cm per ns] // const magnetic field
117  double rad = std::abs(mom.Pt() / (Const::speedOfLight * 1e-4 * 1. * 1.5));
118 
119  // Particle would not be able to reach ECL (underestimation of shower energy)
120  if (2. * rad < pos.Perp()) {
121  rad = pos.Perp() / 2.0 + 1.0;
122  }
123 
124  // Use pq formula (center of circle has to be on perpendicular line through center of line between (0,0) and seed position)
125  double q = pos.Perp();
126  double y3 = pos.Y() / 2.0;
127  double x3 = pos.X() / 2.0;
128 
129  double basex = sqrt(rad * rad - q * q / 4.0) * (-pos.Y()) / q;
130  double basey = sqrt(rad * rad - q * q / 4.0) * pos.X() / q;
131 
132  double centerx1 = x3 + basex;
133  double centery1 = y3 + basey;
134  double centerx2 = x3 - basex;
135  double centery2 = y3 - basey;
136 
137  // vectors for tangent at seed position (perpendicular to radius at seed position)
138  double momx1 = pos.Y() - centery1;
139  double momy1 = - (pos.X() - centerx1);
140  double momx2 = pos.Y() - centery2;
141  double momy2 = - (pos.X() - centerx2);
142 
143  // two solutions (pointing toward and away from center)
144  // make sure that particle moves inwards
145  if (momx1 * pos.X() + momy1 * pos.Y() < 0) {
146  momx1 = -momx1;
147  momy1 = -momy1;
148  }
149  if (momx2 * pos.X() + momy2 * pos.Y() < 0) {
150  momx2 = -momx2;
151  momy2 = -momy2;
152  }
153 
154  // scale to unit length
155  double mom1abs = sqrt(momx1 * momx1 + momy1 * momy1);
156  double mom2abs = sqrt(momx2 * momx2 + momy2 * momy2);
157  momx1 = momx1 / mom1abs;
158  momy1 = momy1 / mom1abs;
159  momx2 = momx2 / mom2abs;
160  momy2 = momy2 / mom2abs;
161 
162  TVector3 mom1(momx1 * mom.Perp(), momy1 * mom.Perp(), mom.Z());
163  TVector3 mom2(momx2 * mom.Perp(), momy2 * mom.Perp(), mom.Z());
164 
165  // Pick the right momentum for positive/negative charge
166  // Use cross product if momentum vector is (counter)clockwise wrt position vector
167  bool clockwise1 = true;
168  bool clockwise2 = true;
169  if (pos.Y() * mom1.X() - pos.X() * mom1.Y() > 0) {
170  clockwise1 = false;
171  }
172  if (pos.Y() * mom2.X() - pos.X() * mom2.Y() > 0) {
173  clockwise2 = false;
174  }
175 
176  if (clockwise1 == clockwise2) {
177  B2WARNING("Something went wrong during helix extrapolation. Skipping track.");
178  continue;
179  }
180 
181  TVector3 mompos;
182  TVector3 momneg;
183  if (clockwise1) {
184  mompos = mom2;
185  momneg = mom1;
186  } else {
187  mompos = mom1;
188  momneg = mom2;
189  }
190 
191  // electron and positron hypothesis
192  RecoTrack* eclSeedNeg = m_eclSeedRecoTracks.appendNew(pos, momneg, -1);
193  eclSeedNeg->addRelationTo(&shower);
194  RecoTrack* eclSeedPos = m_eclSeedRecoTracks.appendNew(pos, mompos, +1);
195  eclSeedPos->addRelationTo(&shower);
196 
197  // define MeasuredStateOnPlane (use pion hypothesis)
198  genfit::AbsTrackRep* repNeg = RecoTrackGenfitAccess::createOrReturnRKTrackRep(*eclSeedNeg, -Const::pion.getPDGCode());
199  genfit::MeasuredStateOnPlane msopNeg(repNeg);
200  genfit::AbsTrackRep* repPos = RecoTrackGenfitAccess::createOrReturnRKTrackRep(*eclSeedPos, Const::pion.getPDGCode());
201  genfit::MeasuredStateOnPlane msopPos(repPos);
202 
203  // set position, momentum, cov, sharedPlanePtr
204  TMatrixDSym cov(6);
205  double covArray[6];
206  shower.getCovarianceMatrixAsArray(covArray);
207 
208  // Calculate uncertainties on position from ECLShower
209  double dx2 = rClus * cosTheta * cosPhi * rClus * cosTheta * cosPhi * covArray[5]
210  + rClus * sinTheta * sinPhi * rClus * sinTheta * sinPhi * covArray[2];
211  double dy2 = rClus * cosTheta * sinPhi * rClus * cosTheta * sinPhi * covArray[5]
212  + rClus * sinTheta * cosPhi * rClus * sinTheta * cosPhi * covArray[2];
213  double dz2 = rClus * sinTheta * rClus * sinTheta * covArray[5];
214 
215  double dpx2 = std::abs(mom1.X() - mom2.X()) / 4.0 * std::abs(mom1.X() - mom2.X()) / 4.0;
216  double dpy2 = std::abs(mom1.Y() - mom2.Y()) / 4.0 * std::abs(mom1.Y() - mom2.Y()) / 4.0;
217  double dpz2 = 0.25 * 0.25 * mom.Z() * mom.Z();
218 
219  cov(0, 0) = dx2;
220  cov(1, 1) = dy2;
221  cov(2, 2) = dz2;
222  cov(3, 3) = dpx2;
223  cov(4, 4) = dpy2;
224  cov(5, 5) = dpz2;
225 
226  // set properties of genfit objects
227  genfit::SharedPlanePtr planeNeg(new genfit::DetPlane(pos, pos));
228  genfit::SharedPlanePtr planePos(new genfit::DetPlane(pos, pos));
229  msopNeg.setPosMomCov(pos, momneg, cov);
230  msopNeg.setPlane(planeNeg);
231  msopPos.setPosMomCov(pos, mompos, cov);
232  msopPos.setPlane(planePos);
233 
234  // create CDCCKF states
235  CDCCKFState seedStateNeg(eclSeedNeg, msopNeg);
236  seeds.push_back({seedStateNeg});
237  CDCCKFState seedStatePos(eclSeedPos, msopPos);
238  seeds.push_back({seedStatePos});
239  }
240 }
Belle2::CDCCKFEclSeedCreator::m_param_restrictToForwardSeeds
bool m_param_restrictToForwardSeeds
Don't do Ecl seeding in central region to save computing time.
Definition: CDCCKFEclSeedCreator.h:69
Belle2::RecoTrack::registerRequiredRelations
static void registerRequiredRelations(StoreArray< RecoTrack > &recoTracks, std::string const &pxdHitsStoreArrayName="", std::string const &svdHitsStoreArrayName="", std::string const &cdcHitsStoreArrayName="", std::string const &bklmHitsStoreArrayName="", std::string const &eklmHitsStoreArrayName="", std::string const &recoHitInformationStoreArrayName="")
Convenience method which registers all relations required to fully use a RecoTrack.
Definition: RecoTrack.cc:42
Belle2::CDCCKFEclSeedCreator::exposeParameters
void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix) override
Expose the parameters of the sub findlets.
Definition: CDCCKFEclSeedCreator.cc:28
genfit::SharedPlanePtr
std::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
Definition: SharedPlanePtr.h:40
Belle2::CDCCKFEclSeedCreator::m_param_minimalEnRequirement
double m_param_minimalEnRequirement
Minimal pt requirement.
Definition: CDCCKFEclSeedCreator.h:67
genfit::MeasuredStateOnPlane
#StateOnPlane with additional covariance matrix.
Definition: MeasuredStateOnPlane.h:39
Belle2::CDCCKFEclSeedCreator::CDCCKFEclSeedCreator
CDCCKFEclSeedCreator()
Add the subfindlets.
Definition: CDCCKFEclSeedCreator.cc:23
Belle2::CDCCKFEclSeedCreator::m_inputECLshowers
StoreArray< ECLShower > m_inputECLshowers
Input ECL Showers Store Array.
Definition: CDCCKFEclSeedCreator.h:79
Belle2::TrackFindingCDC::Findlet< CDCCKFPath >
Belle2::CDCCKFEclSeedCreator::apply
void apply(std::vector< CDCCKFPath > &seeds) override
Load in the reco tracks and the hits.
Definition: CDCCKFEclSeedCreator.cc:76
Belle2::RelationsInterface::addRelationTo
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
Definition: RelationsObject.h:144
genfit::AbsTrackRep
Abstract base class for a track representation.
Definition: AbsTrackRep.h:66
Belle2::CDCCKFEclSeedCreator::m_param_tanLambdaForwardNeg
double m_param_tanLambdaForwardNeg
Up to which (neg) tanLambda value should the seeding be performed.
Definition: CDCCKFEclSeedCreator.h:73
Belle2::TrackFindingCDC::CompositeProcessingSignalListener::initialize
void initialize() override
Receive and dispatch signal before the start of the event processing.
Definition: CompositeProcessingSignalListener.cc:17
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:568
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:535
Belle2::CDCCKFEclSeedCreator::m_param_tanLambdaForwardPos
double m_param_tanLambdaForwardPos
Up to which (pos) tanLambda value should the seeding be performed.
Definition: CDCCKFEclSeedCreator.h:75
Belle2::CDCCKFEclSeedCreator::initialize
void initialize() override
Create the store arrays.
Definition: CDCCKFEclSeedCreator.cc:64
Belle2::ModuleParamList::addParameter
void addParameter(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module list.
Definition: ModuleParamList.templateDetails.h:38
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDCCKFEclSeedCreator::m_eclSeedRecoTracks
StoreArray< RecoTrack > m_eclSeedRecoTracks
Output Reco Tracks Store Array.
Definition: CDCCKFEclSeedCreator.h:81
Belle2::ECLShower::c_nPhotons
@ c_nPhotons
CR is split into n photons (N1)
Definition: ECLShower.h:54
Belle2::RecoTrackGenfitAccess::createOrReturnRKTrackRep
static genfit::AbsTrackRep * createOrReturnRKTrackRep(RecoTrack &recoTrack, int PDGcode)
Checks if a TrackRap for the PDG id of the RecoTrack (and its charge conjugate) does already exit and...
Definition: RecoTrack.cc:394
Belle2::CDCCKFEclSeedCreator::m_param_showerDepth
double m_param_showerDepth
Correction if the shower is assumed to start in a certain depth.
Definition: CDCCKFEclSeedCreator.h:71
genfit::DetPlane
Detector plane.
Definition: DetPlane.h:59
Belle2::CDCCKFState
Define states for CKF algorithm, which can be seed track or CDC wire hit.
Definition: CDCCKFState.h:37
Belle2::ModuleParamList
The Module parameter list class.
Definition: ModuleParamList.h:46
Belle2::CDCCKFEclSeedCreator::m_param_inputEclShowerStoreArrayName
std::string m_param_inputEclShowerStoreArrayName
StoreArray name of the input Ecl Shower Store Array.
Definition: CDCCKFEclSeedCreator.h:63
Belle2::CDCCKFEclSeedCreator::m_param_eclSeedRecoTrackStoreArrayName
std::string m_param_eclSeedRecoTrackStoreArrayName
StoreArray name of the output Track Store Array.
Definition: CDCCKFEclSeedCreator.h:65