9#include "cdc/modules/cdcCalibrationCollector/CDCCalibrationCollector.h"
10#include <cdc/translators/RealisticTDCCountTranslator.h>
11#include <framework/datastore/RelationArray.h>
13#include <tracking/trackFindingCDC/topology/CDCWireTopology.h>
15#include <genfit/TrackPoint.h>
16#include <genfit/KalmanFitterInfo.h>
17#include <genfit/MeasurementOnPlane.h>
18#include <genfit/MeasuredStateOnPlane.h>
20#include <Math/ProbFuncMathCore.h>
22#include <cdc/dataobjects/WireID.h>
23#include <cdc/geometry/CDCGeometryPar.h>
32using namespace genfit;
33using namespace TrackFindingCDC;
49 addParam(
"minimumPt",
m_minimumPt,
"Tracks with transverse momentum smaller than this value will not used", 0.15);
50 addParam(
"minimumNDF",
m_minimumNDF,
"Discard tracks whose degree-of-freedom below this value", 5.);
52 addParam(
"effStudy",
m_effStudy,
"When true module collects info only necessary for wire eff study",
false);
71 auto m_tree =
new TTree(
m_treeName.c_str(),
"tree for cdc calibration");
72 m_tree->Branch<Float_t>(
"x_mea", &
x_mea);
73 m_tree->Branch<Float_t>(
"x_u", &
x_u);
74 m_tree->Branch<Float_t>(
"x_b", &
x_b);
75 m_tree->Branch<Float_t>(
"alpha", &
alpha);
76 m_tree->Branch<Float_t>(
"theta", &
theta);
77 m_tree->Branch<Float_t>(
"t", &
t);
78 m_tree->Branch<UShort_t>(
"adc", &
adc);
80 m_tree->Branch<UChar_t>(
"lay", &
lay);
81 m_tree->Branch<Float_t>(
"weight", &
weight);
82 m_tree->Branch<UShort_t>(
"IWire", &
IWire);
83 m_tree->Branch<Float_t>(
"Pval", &
Pval);
84 m_tree->Branch<Float_t>(
"ndf", &
ndf);
86 m_tree->Branch<Float_t>(
"d0", &
d0);
87 m_tree->Branch<Float_t>(
"z0", &
z0);
88 m_tree->Branch<Float_t>(
"phi0", &
phi0);
89 m_tree->Branch<Float_t>(
"tanL", &
tanL);
90 m_tree->Branch<Float_t>(
"omega", &
omega);
94 m_tree->Branch<Float_t>(
"t_fit", &
t_fit);
100 auto m_efftree =
new TTree(
m_effTreeName.c_str(),
"tree for wire efficiency");
101 m_efftree->Branch<
unsigned short>(
"layerID", &
layerID);
102 m_efftree->Branch<
unsigned short>(
"wireID", &
wireID);
103 m_efftree->Branch<
float>(
"z", &
z);
104 m_efftree->Branch<
bool>(
"isFound", &
isFound);
109 auto m_hNDF =
new TH1F(
"hNDF",
"NDF of fitted track;NDF;Tracks", 71, -1, 70);
110 auto m_hPval =
new TH1F(
"hPval",
"p-values of tracks;pVal;Tracks", 1000, 0, 1);
111 auto m_hEventT0 =
new TH1F(
"hEventT0",
"Event T0", 1000, -100, 100);
112 auto m_hNTracks =
new TH1F(
"hNTracks",
"Number of tracks", 50, 0, 10);
113 auto m_hOccupancy =
new TH1F(
"hOccupancy",
"occupancy", 100, 0, 1.0);
121 ROOT::Math::XYZVector pos(0, 0, 0);
123 if (bfield.Z() > 0.5) {
125 B2INFO(
"CDCCalibrationCollector: Magnetic field is ON");
128 B2INFO(
"CDCCalibrationCollector: Magnetic field is OFF");
130 B2INFO(
"BField at (0,0,0) = " << bfield.R());
154 std::vector<unsigned short> wiresInCDCTrack;
158 unsigned short eWireID = cdcHit.getWire().getEWire();
159 wiresInCDCTrack.push_back(eWireID);
164 const int nTr =
m_Tracks.getEntries();
165 const int nHits =
m_CDCHits.getEntries();
166 const int nWires = 14336;
167 float oc =
static_cast<float>(nHits) /
static_cast<float>(nWires);
171 for (
int i = 0; i < nTr; ++i) {
175 B2WARNING(
"No track fit result found.");
181 B2WARNING(
"Can not access RecoTrack of this Belle2::Track");
193 B2DEBUG(99,
"ndf = " <<
ndf);
194 B2DEBUG(99,
"Pval = " <<
Pval);
197 double Chi2 = fs->getChi2();
198 Pval = std::max(0., ROOT::Math::chisquared_cdf_c(Chi2,
ndf));
216 }
catch (
const genfit::Exception& e) {
217 B2ERROR(
"Exception when harvest information from recotrack: " << e.what());
222 if (fitresult->
getD0() > 2 || fitresult->
getZ0() > 5)
continue;
236 B2DEBUG(99,
"start collect hit");
241 const genfit::TrackPoint* tp = track->getCreatedTrackPoint(track->getRecoHitInformation(hit));
243 lay = hit->getICLayer();
244 IWire = hit->getIWire();
245 adc = hit->getADCCount();
246 unsigned short tdc = hit->getTDCCount();
248 const genfit::KalmanFitterInfo* kfi = tp->getKalmanFitterInfo();
249 if (!kfi) {B2DEBUG(199,
"No Fitter Info: Layer " << hit->getICLayer());
continue;}
250 for (
unsigned int iMeas = 0; iMeas < kfi->getNumMeasurements(); ++iMeas) {
251 if ((kfi->getWeights().at(iMeas)) > 0.5) {
253 const genfit::MeasuredStateOnPlane& mop = kfi->getFittedState();
254 const TVector3 pocaOnWire = mop.getPlane()->getO();
255 const TVector3 pocaOnTrack = mop.getPlane()->getU();
256 const TVector3 pocaMom = mop.getMom();
259 x_mea = kfi->getMeasurementOnPlane(iMeas)->getState()(0);
260 x_b = kfi->getFittedState(
true).getState()(3);
261 x_u = kfi->getFittedState(
false).getState()(3);
262 weight = kfi->getWeights().at(iMeas);
269 if (std::abs(
alpha) > M_PI / 2) {
278 B2DEBUG(99,
"x_unbiased " <<
x_u <<
" |left_right " << lr);
283 t = tdcTrans->
getDriftTime(tdc, wireid, mop.getTime(), pocaOnWire.Z(),
adc);
292 const Helix& helixFit)
const
298 const CDCWire& oneWire = layer.getWire(1);
300 double arcLength = helixFit.getArcLength2DAtCylindricalR(newR);
301 ROOT::Math::XYZVector xyzOnWire =
B2Vector3D(helixFit.getPositionAtArcLength2D(arcLength));
305 const CDCWire& wire = layer.getClosestWire(crosspoint);
314 const double radiusofLayer = wireLayer.getRefCylindricalR();
316 const double arcLength = helixFit.getArcLength2DAtCylindricalR(radiusofLayer);
317 const ROOT::Math::XYZVector xyz =
B2Vector3D(helixFit.getPositionAtArcLength2D(arcLength));
318 if (!xyz.X())
continue;
320 unsigned short crossedWire = wireIntersected.
getEWire();
321 unsigned short crossedCWire = wireIntersected.
getNeighborCW()->getEWire();
322 unsigned short crossedCCWire = wireIntersected.
getNeighborCCW()->getEWire();
324 if (find(wireHits.begin(), wireHits.end(), crossedWire) != wireHits.end()
325 || find(wireHits.begin(), wireHits.end(), crossedCWire) != wireHits.end()
326 || find(wireHits.begin(), wireHits.end(), crossedCCWire) != wireHits.end())
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
StoreObjPtr< EventT0 > m_eventTimeStoreObject
Event t0 object.
std::string m_recoTrackArrayName
Belle2::RecoTrack StoreArray name.e.
bool m_calExpectedDriftTime
Calculate expected drift time from x_fit or not.
Float_t phi0
Track Parameter, phi0.
Float_t tanL
Track Parameter, tanL.
double m_minimumNDF
minimum NDF required for track
std::string m_effTreeName
Name of efficiency tree for the output file.
std::string m_cdcTrackVectorName
Belle2::CDCTrack vectorpointer name.
bool m_storeTrackParams
Store Track parameter or not.
Float_t weight
Weight of hit.
TrackFindingCDC::StoreWrappedObjPtr< std::vector< TrackFindingCDC::CDCTrack > > m_CDCTracks
CDC tracks.
StoreArray< TrackFitResult > m_TrackFitResults
Track fit results.
void harvest(Belle2::RecoTrack *track)
collect hit information of fitted track.
CDCCalibrationCollectorModule()
Constructor.
std::string m_cdcHitArrayName
Belle2::CDCHit StoreArray name.
Float_t theta
Entrance Polar angle of hit (degree).
Float_t x_b
X_fit for biased track fit.
Float_t Pval
P-value of fitted track.
std::string m_relRecoTrackTrackName
Relation between RecoTrack and Belle2:Track.
unsigned short wireID
wireID for hit-level wire monitoring
void collect() override
Event action, collect information for calibration.
bool m_isCosmic
true when we process cosmic events, else false (collision).
Float_t t_fit
Drift time calculated from x_fit.
virtual ~CDCCalibrationCollectorModule()
Destructor.
Float_t omega
Track Parameter, omega.
bool m_eventT0Extraction
use Event T0 extract t0 or not.
void buildEfficiencies(std::vector< unsigned short > wireHits, const Helix helixFit)
fills efficiency objects
bool m_bField
fit incase no magnetic Field of not, if false, NDF=4 in cal P-value
std::string m_trackArrayName
Belle2::Track StoreArray name.
void prepare() override
Initializes the Module.
Float_t z0
Track Parameter, z0.
Float_t t
Measurement Drift time.
Float_t x_u
X_fit for unbiased track fit.
void finish() override
Termination action.
Float_t ndf
degree of freedom.
std::string m_trackFitResultArrayName
Belle2::TrackFitResult StoreArray name.
Float_t d0
Track Parameter, d0.
unsigned short layerID
layerID for hit-level wire monitoring
Float_t x_mea
measure drift length (signed by left right).
const TrackFindingCDC::CDCWire & getIntersectingWire(const ROOT::Math::XYZVector &xyz, const TrackFindingCDC::CDCWireLayer &layer, const Helix &helixFit) const
extrapolates the helix fit to a given layer and finds the wire which it would be hitting
std::string m_treeName
Name of tree for the output file.
StoreArray< RecoTrack > m_RecoTracks
Tracks.
StoreArray< Track > m_Tracks
Tracks.
Float_t alpha
Entrance Azimuthal angle of hit (degree).
double m_minimumPt
minimum pt required for track
bool m_effStudy
When true module collects info only necessary for wire eff study.
float z
z of hit for hit-level wire monitoring
bool isFound
flag for a hit that has been found near a track as expected by extrapolation
StoreArray< CDCHit > m_CDCHits
CDC hits.
The Class for CDC Geometry Parameters.
double getTheta(const B2Vector3D &momentum) const
Returns track incident angle (theta in rad.).
double getAlpha(const B2Vector3D &posOnWire, const B2Vector3D &momentum) const
Returns track incident angle in rphi plane (alpha in rad.).
double getOutgoingAlpha(const double alpha) const
Converts incoming- to outgoing-alpha.
unsigned short getOutgoingLR(const unsigned short lr, const double alpha) const
Converts incoming-lr to outgoing-lr.
double getOutgoingTheta(const double alpha, const double theta) const
Converts incoming- to outgoing-theta.
double getDriftTime(double dist, unsigned short layer, unsigned short lr, double alpha, double theta) const
Return the drift time to the sense wire.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Translator mirroring the realistic Digitization.
double getDriftTime(unsigned short tdcCount, const WireID &wireID, double timeOfFlightEstimator, double z, unsigned short adcCount) override
Get Drift time.
void registerObject(std::string name, T *obj)
Register object with a name, takes ownership, do not access the pointer beyond prepare()
CalibrationCollectorModule()
Constructor. Sets the default prefix for calibration dataobjects.
T * getObjectPtr(std::string name)
Calls the CalibObjManager to get the requested stored collector data.
static const ChargedStable muon
muon particle
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
This is the Reconstruction Event-Data Model Track.
const genfit::FitStatus * getTrackFitStatus(const genfit::AbsTrackRep *representation=nullptr) const
Return the track fit status for the given representation or for the cardinal one. You are not allowed...
Low-level class to create/modify relations between StoreArrays.
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
const std::string & getName() const
Return name under which the object is saved in the DataStore.
Class representing a three dimensional reconstructed hit.
Class representing a sequence of three dimensional reconstructed hits.
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.
Class representing a sense wire in the central drift chamber.
Vector2D getWirePos2DAtZ(const double z) const
Gives the xy projected position of the wire at the given z coordinate.
MayBePtr< const CDCWire > getNeighborCCW() const
Gives the closest neighbor in the counterclockwise direction - always exists.
IWire getIWire() const
Getter for the wire id within its layer.
unsigned short getEWire() const
Getter for the encoded wire number.
ILayer getICLayer() const
Getter for the continuous layer id ranging from 0 - 55.
MayBePtr< const CDCWire > getNeighborCW() const
Gives the closest neighbor in the clockwise direction - always exists.
double norm() const
Calculates the length of the vector.
Values of the result of a track fit with a given particle hypothesis.
Helix getHelix() const
Conversion to framework Helix (without covariance).
double getOmega() const
Getter for omega.
double getD0() const
Getter for d0.
double getTransverseMomentum() const
Getter for the absolute value of the transverse momentum at the perigee.
double getTanLambda() const
Getter for tanLambda.
double getZ0() const
Getter for z0.
double getPhi0() const
Getter for phi0.
Class that bundles various TrackFitResults.
const TrackFitResult * getTrackFitResultWithClosestMass(const Const::ChargedStable &requestedType) const
Return the track fit for the fit hypothesis with the closest mass.
Class to identify a wire inside the CDC.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
B2Vector3< double > B2Vector3D
typedef for common usage with double
HepGeom::Vector3D< double > Vector3D
3D Vector
Abstract base class for different kinds of events.