9 #include <ecl/modules/eclTrackBremFinder/ECLTrackBremFinderModule.h>
12 #include <framework/dataobjects/EventMetaData.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/gearbox/Const.h>
15 #include <framework/datastore/StoreArray.h>
18 #include <mdst/dataobjects/ECLCluster.h>
19 #include <mdst/dataobjects/Track.h>
20 #include <mdst/dataobjects/PIDLikelihood.h>
21 #include <mdst/dataobjects/TrackFitResult.h>
24 #include <tracking/dataobjects/RecoTrack.h>
25 #include <tracking/dataobjects/BremHit.h>
28 #include <ecl/modules/eclTrackBremFinder/BestMatchContainer.h>
29 #include <ecl/modules/eclTrackBremFinder/BremFindingMatchCompute.h>
39 setDescription(
"Use Track direction to pick up possible ECL Brem Cluster");
40 setPropertyFlags(c_ParallelProcessingCertified);
42 addParam(
"eclClustersStoreArrayName", m_param_eclClustersStoreArrayName,
"StoreArray name of the ECLClusters for brem matching",
43 m_param_eclClustersStoreArrayName);
45 addParam(
"tracksStoreArrayName", m_param_tracksStoreArrayName,
"StoreArray name of the Tracks for brem matching",
46 m_param_tracksStoreArrayName);
48 addParam(
"recoTracksStoreArrayName", m_param_recoTracksStoreArrayName,
"StoreArray name of the reco tracks used for brem search",
49 m_param_recoTracksStoreArrayName);
51 addParam(
"clusterAcceptanceFactor", m_clusterAcceptanceFactor,
52 "Factor which is multiplied onto the cluster position error to check for matches",
53 m_clusterAcceptanceFactor);
55 addParam(
"virtualHitRadii", m_virtualHitRadii,
"Radii where virtual hits for the extrapolation will be generated",
58 addParam(
"relativeClusterEnergy", m_relativeClusterEnergy,
"Fraction of the tracks energy the ECL cluster has "
59 "to possess to be considered for bremsstrahlung finding",
60 m_relativeClusterEnergy);
62 addParam(
"requestedNumberOfCDCHits", m_requestedNumberOfCDCHits,
"Minimal/Maximal number of CDC hits, the track has to possess "
63 "to be considered for bremsstrahlung finding",
64 m_requestedNumberOfCDCHits);
66 addParam(
"electronProbabilityCut", m_electronProbabilityCut,
"Cut on the electron probability (from pid) of track",
67 m_electronProbabilityCut);
69 addParam(
"clusterDistanceCut", m_clusterDistanceCut,
70 "Cut on the distance between the cluster position angle and the extrapolation angle",
71 m_clusterDistanceCut);
81 const std::string relationName =
"Bremsstrahlung";
108 double highestProb = 0;
110 double probability = pid->getProbability(pdg);
111 if (probability > highestProb) {
112 highestProb = probability;
117 B2DEBUG(20,
"Track is expected not to be from electron");
123 double trackMomentum;
124 if (trackFitResult) {
125 trackMomentum = trackFitResult->
getMomentum().Mag();
127 if (trackMomentum > 5.0) {
128 B2DEBUG(20,
"Track momentum higher than 5GeV! Track is not used for bremsstrahlung finding");
136 B2DEBUG(20,
"Checking track for related ECLCluster");
140 auto relatedClustersToTrack =
143 for (
auto& relatedCluster : relatedClustersToTrack) {
145 primaryClusterOfTrack = &relatedCluster;
155 B2DEBUG(20,
"No RecoTrack for this Track");
162 B2DEBUG(20,
"RecoTrack has not requested number of CDC hits");
167 std::vector<std::pair<float, RecoHitInformation*>> extrapolationParams = {};
168 std::vector<genfit::MeasuredStateOnPlane> extrapolatedStates = {};
172 for (
auto hit : recoTrack->getRecoHitInformations(
true)) {
173 if (hit->useInFit() && recoTrack->hasTrackFitStatus()) {
175 auto measState = recoTrack->getMeasuredStateOnPlaneFromRecoHit(hit);
176 float hitRadius = measState.getPos().Perp();
177 float distance = abs(hitRadius - virtualHitRadius);
180 nearestHitContainer.
add(hit, distance);
182 }
catch (NoTrackFitResult&) {
183 B2DEBUG(29,
"No track fit result available for this hit! Event: " <<
m_evtPtr->getEvent());
186 if (nearestHitContainer.
hasMatch()) {
188 extrapolationParams.push_back({virtualHitRadius, nearestHit});
193 for (
auto param : extrapolationParams) {
194 auto fitted_state = recoTrack->getMeasuredStateOnPlaneFromRecoHit(param.second);
196 fitted_state.extrapolateToCylinder(param.first);
197 extrapolatedStates.push_back(fitted_state);
199 B2DEBUG(20,
"Extrapolation failed!");
203 B2WARNING(
"Exception" << e.what());
211 B2DEBUG(20,
"Cluster has wrong hypothesis!");
216 auto relatedTrack = cluster.getRelatedFrom<
Track>();
218 B2DEBUG(20,
"Cluster already related to track, bailing out");
224 auto relatedBremHit = cluster.getRelated<
BremHit>();
225 if (relatedBremHit) {
226 B2DEBUG(20,
"Cluster already assumed to be bremsstrahlung cluster!");
234 B2DEBUG(20,
"Relative energy of cluster higher than 1 or below threshold!");
238 typedef std::tuple<ECLCluster*, genfit::MeasuredStateOnPlane, double, double > ClusterMSoPPair;
244 for (
auto hit : recoTrack->getRecoHitInformations(
true)) {
245 if (hit->getTrackingDetector() == RecoHitInformation::c_PXD || hit->getTrackingDetector() == RecoHitInformation::c_SVD) {
247 if (!recoTrack->hasTrackFitStatus()) {
250 auto measState = recoTrack->getMeasuredStateOnPlaneFromRecoHit(hit);
252 if (bremFinder.isMatch()) {
253 ClusterMSoPPair match_pair = std::make_tuple(&cluster, measState, bremFinder.getDistanceHitCluster(),
254 bremFinder.getEffAcceptanceFactor());
255 matchContainer.
add(match_pair, bremFinder.getDistanceHitCluster());
257 }
catch (NoTrackFitResult&) {
258 B2DEBUG(29,
"No track fit result available for this hit! Event: " <<
m_evtPtr->getEvent());
260 B2WARNING(
"Exception" << e.what());
266 for (
auto fitted_state : extrapolatedStates) {
268 if (bremFinder.isMatch()) {
269 ClusterMSoPPair match_pair = std::make_tuple(&cluster, fitted_state, bremFinder.getDistanceHitCluster(),
270 bremFinder.getEffAcceptanceFactor());
271 matchContainer.
add(match_pair, bremFinder.getDistanceHitCluster());
279 const auto fitted_state = std::get<1>(matchClustermSoP);
281 const auto fitted_pos = fitted_state.getPos();
282 const auto fitted_mom = fitted_state.getMom();
283 const auto fitted_dir = fitted_state.getDir();
285 const auto hit_theta = fitted_mom.Theta();
286 const auto hit_phi = fitted_mom.Phi();
288 B2DEBUG(20,
"Best Cluster" << std::endl
289 <<
" Cluster Phi=" << std::get<0>(matchClustermSoP)->getPhi() <<
" Theta=" << std::get<0>(matchClustermSoP)->getTheta()
290 <<
" TrackHit Phi=" << hit_phi <<
" Theta=" << hit_theta);
296 double effAcceptanceFactor = std::get<3>(matchClustermSoP);
297 ECLCluster* bremCluster = std::get<0>(matchClustermSoP);
298 double clusterDistance = std::get<2>(matchClustermSoP);
303 clusterDistance, effAcceptanceFactor);
307 bremCluster->
addRelationTo(&track, effAcceptanceFactor,
"Bremsstrahlung");
309 if (primaryClusterOfTrack) {
310 primaryClusterOfTrack->
addRelationTo(bremCluster, effAcceptanceFactor);