Belle II Software development
CDCCalibrationCollector.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include "cdc/modules/cdcCalibrationCollector/CDCCalibrationCollector.h"
10#include <cdc/translators/RealisticTDCCountTranslator.h>
11#include <framework/datastore/RelationArray.h>
12
13#include <tracking/trackFindingCDC/topology/CDCWireTopology.h>
14
15#include <genfit/TrackPoint.h>
16#include <genfit/KalmanFitterInfo.h>
17#include <genfit/MeasurementOnPlane.h>
18#include <genfit/MeasuredStateOnPlane.h>
19
20#include <Math/ProbFuncMathCore.h>
21
22#include <cdc/dataobjects/WireID.h>
23#include <cdc/geometry/CDCGeometryPar.h>
24
25#include <TH1F.h>
26
27#include <cmath>
28
29//using namespace std;
30using namespace Belle2;
31using namespace CDC;
32using namespace genfit;
33using namespace TrackFindingCDC;
34
35
36REG_MODULE(CDCCalibrationCollector);
37
38
40{
41 setDescription("Collector module for cdc calibration");
42 setPropertyFlags(c_ParallelProcessingCertified); // specify this flag if you need parallel processing
43 addParam("recoTracksColName", m_recoTrackArrayName, "Name of collection hold genfit::Track", std::string(""));
44 //addParam("bField", m_bField, "If true -> #Params ==5 else #params ==4 for calculate P-Val", false);
45 addParam("calExpectedDriftTime", m_calExpectedDriftTime, "if true module will calculate expected drift time, it take a time",
46 true);
47 addParam("storeTrackParams", m_storeTrackParams, "Store Track Parameter or not, it will be multicount for each hit", false);
48 addParam("eventT0Extraction", m_eventT0Extraction, "use event t0 extract t0 or not", true);
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.);
51 addParam("isCosmic", m_isCosmic, "True when we process cosmic events, else False (collision)", m_isCosmic);
52 addParam("effStudy", m_effStudy, "When true module collects info only necessary for wire eff study", false);
53}
54
58
60{
61 m_Tracks.isRequired(m_trackArrayName);
64 m_CDCHits.isRequired(m_cdcHitArrayName);
67 //Store names to speed up creation later
68 m_relRecoTrackTrackName = relRecoTrackTrack.getName();
69
70 if (!m_effStudy) { // by default collects calibration data
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);
79 // m_tree->Branch<int>("boardID", &boardID);
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);
91 }
92
93 if (m_calExpectedDriftTime) { // expected drift time, calculated form xfit
94 m_tree->Branch<Float_t>("t_fit", &t_fit);
95 }
96
97 registerObject<TTree>("tree", m_tree);
98 }
99 if (m_effStudy) { //if m_effStudy is changed to true prepares to only run wire efficiency study
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);
105
106 registerObject<TTree>("efftree", m_efftree);
107 }
108
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);
114
115 registerObject<TH1F>("hNDF", m_hNDF);
116 registerObject<TH1F>("hPval", m_hPval);
117 registerObject<TH1F>("hEventT0", m_hEventT0);
118 registerObject<TH1F>("hNTracks", m_hNTracks);
119 registerObject<TH1F>("hOccupancy", m_hOccupancy);
120
121 ROOT::Math::XYZVector pos(0, 0, 0);
122 ROOT::Math::XYZVector bfield = BFieldManager::getFieldInTesla(pos);
123 if (bfield.Z() > 0.5) {
124 m_bField = true;
125 B2INFO("CDCCalibrationCollector: Magnetic field is ON");
126 } else {
127 m_bField = false;
128 B2INFO("CDCCalibrationCollector: Magnetic field is OFF");
129 }
130 B2INFO("BField at (0,0,0) = " << bfield.R());
131
132}
133
135{
137
138 /* CDCHit distribution */
139 // make evt t0 incase we dont use evt t0
140 evtT0 = 0;
141
142 //reject events don't have eventT0
144 // event with is fail to extract t0 will be exclude from analysis
145 if (m_eventTimeStoreObject.isValid() && m_eventTimeStoreObject->hasEventT0()) {
146 evtT0 = m_eventTimeStoreObject->getEventT0();
147 getObjectPtr<TH1F>("hEventT0")->Fill(evtT0);
148 } else {
149 return;
150 }
151 }
152 // Collects the WireID and Layer of every hit in this event
153 // Used in wire efficiency building
154 std::vector<unsigned short> wiresInCDCTrack;
155
156 for (CDCTrack& cdcTrack : *m_CDCTracks) {
157 for (CDCRecoHit3D& cdcHit : cdcTrack) {
158 unsigned short eWireID = cdcHit.getWire().getEWire();
159 wiresInCDCTrack.push_back(eWireID);
160 }
161 }
162 // WireID collection finished
163
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);
168 getObjectPtr<TH1F>("hNTracks")->Fill(nTr);
169 getObjectPtr<TH1F>("hOccupancy")->Fill(oc);
170
171 for (int i = 0; i < nTr; ++i) {
172 const Belle2::Track* b2track = m_Tracks[i];
174 if (!fitresult) {
175 B2WARNING("No track fit result found.");
176 continue;
177 }
178
180 if (!recoTrack) {
181 B2WARNING("Can not access RecoTrack of this Belle2::Track");
182 continue;
183 }
184 const genfit::FitStatus* fs = recoTrack->getTrackFitStatus();
185 if (!fs) continue;
186 ndf = fs->getNdf();
187 if (!m_bField) {
188 ndf += 1;
189 }
190
191 getObjectPtr<TH1F>("hPval")->Fill(Pval);
192 getObjectPtr<TH1F>("hNDF")->Fill(ndf);
193 B2DEBUG(99, "ndf = " << ndf);
194 B2DEBUG(99, "Pval = " << Pval);
195
196 if (ndf < m_minimumNDF) continue;
197 double Chi2 = fs->getChi2();
198 Pval = std::max(0., ROOT::Math::chisquared_cdf_c(Chi2, ndf));
199 //store track parameters
200
201 d0 = fitresult->getD0();
202 z0 = fitresult->getZ0();
203 tanL = fitresult->getTanLambda();
204 omega = fitresult->getOmega();
205 phi0 = fitresult->getPhi0() * 180 / M_PI;
206
207 // Rejection of suspicious cosmic tracks.
208 // phi0 of cosmic track must be negative in our definition!
209 if (m_isCosmic == true && phi0 > 0.0) continue;
210
211 //cut at Pt
212 if (fitresult->getTransverseMomentum() < m_minimumPt) continue;
213 if (!m_effStudy) { // all harvest to fill the tree if collecting calibration info
214 try {
215 harvest(recoTrack);
216 } catch (const genfit::Exception& e) {
217 B2ERROR("Exception when harvest information from recotrack: " << e.what());
218 }
219 }
220 if (m_effStudy) { // call buildEfficiencies for efficiency study
221 // Request tracks coming from IP
222 if (fitresult->getD0() > 2 || fitresult->getZ0() > 5) continue;
223 const Helix helixFit = fitresult->getHelix();
224 buildEfficiencies(wiresInCDCTrack, helixFit);
225 }
226 }
227
228}
229
233
235{
236 B2DEBUG(99, "start collect hit");
237 static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
239
240 for (const RecoHitInformation::UsedCDCHit* hit : track->getCDCHitList()) {
241 const genfit::TrackPoint* tp = track->getCreatedTrackPoint(track->getRecoHitInformation(hit));
242 if (!tp) continue;
243 lay = hit->getICLayer();
244 IWire = hit->getIWire();
245 adc = hit->getADCCount();
246 unsigned short tdc = hit->getTDCCount();
247 WireID wireid(lay, IWire);
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) {
252 // int boardID = cdcgeo.getBoardID(WireID(lay,IWire));
253 const genfit::MeasuredStateOnPlane& mop = kfi->getFittedState();
254 const TVector3 pocaOnWire = mop.getPlane()->getO();//Local wire position
255 const TVector3 pocaOnTrack = mop.getPlane()->getU();//residual direction
256 const TVector3 pocaMom = mop.getMom();
257 alpha = cdcgeo.getAlpha(pocaOnWire, pocaMom) ;
258 theta = cdcgeo.getTheta(pocaMom);
259 x_mea = kfi->getMeasurementOnPlane(iMeas)->getState()(0);
260 x_b = kfi->getFittedState(true).getState()(3);// x fit biased
261 x_u = kfi->getFittedState(false).getState()(3);//x fit unbiased
262 weight = kfi->getWeights().at(iMeas);
263
264 int lr;
265 if (x_u > 0) lr = 1;
266 else lr = 0;
267
268 //Convert to outgoing
269 if (std::abs(alpha) > M_PI / 2) {
270 x_b *= -1;
271 x_u *= -1;
272 }
273 x_mea = copysign(x_mea, x_b);
274 lr = cdcgeo.getOutgoingLR(lr, alpha);
276 alpha = cdcgeo.getOutgoingAlpha(alpha);
277
278 B2DEBUG(99, "x_unbiased " << x_u << " |left_right " << lr);
279 if (m_calExpectedDriftTime) { t_fit = cdcgeo.getDriftTime(std::abs(x_u), lay, lr, alpha, theta);}
280 alpha *= 180 / M_PI;
281 theta *= 180 / M_PI;
282 //estimate drift time
283 t = tdcTrans->getDriftTime(tdc, wireid, mop.getTime(), pocaOnWire.Z(), adc);
284 getObjectPtr<TTree>("tree")->Fill();
285 } //NDF
286 // }//end of if isU
287 }//end of for
288 }//end of for tp
289}//end of func
290
291const CDCWire& CDCCalibrationCollectorModule::getIntersectingWire(const ROOT::Math::XYZVector& xyz, const CDCWireLayer& layer,
292 const Helix& helixFit) const
293{
294 Vector3D crosspoint;
295 if (layer.isAxial())
296 crosspoint = Vector3D(xyz);
297 else {
298 const CDCWire& oneWire = layer.getWire(1);
299 double newR = oneWire.getWirePos2DAtZ(xyz.Z()).norm();
300 double arcLength = helixFit.getArcLength2DAtCylindricalR(newR);
301 ROOT::Math::XYZVector xyzOnWire = B2Vector3D(helixFit.getPositionAtArcLength2D(arcLength));
302 crosspoint = Vector3D(xyzOnWire);
303 }
304
305 const CDCWire& wire = layer.getClosestWire(crosspoint);
306
307 return wire;
308}
309
310void CDCCalibrationCollectorModule::buildEfficiencies(std::vector<unsigned short> wireHits, const Helix helixFit)
311{
313 for (const CDCWireLayer& wireLayer : wireTopology.getWireLayers()) {
314 const double radiusofLayer = wireLayer.getRefCylindricalR();
315 //simple extrapolation of fit
316 const double arcLength = helixFit.getArcLength2DAtCylindricalR(radiusofLayer);
317 const ROOT::Math::XYZVector xyz = B2Vector3D(helixFit.getPositionAtArcLength2D(arcLength));
318 if (!xyz.X()) continue;
319 const CDCWire& wireIntersected = getIntersectingWire(xyz, wireLayer, helixFit);
320 unsigned short crossedWire = wireIntersected.getEWire();
321 unsigned short crossedCWire = wireIntersected.getNeighborCW()->getEWire();
322 unsigned short crossedCCWire = wireIntersected.getNeighborCCW()->getEWire();
323
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())
327 isFound = true;
328 else
329 isFound = false;
330
331 wireID = wireIntersected.getIWire();
332 layerID = wireIntersected.getICLayer();
333 z = xyz.Z();
334 getObjectPtr<TTree>("efftree")->Fill();
335 }
336}
337
338
339
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.
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.
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.
std::string m_cdcHitArrayName
Belle2::CDCHit StoreArray name.
Float_t theta
Entrance Polar angle of hit (degree).
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.
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 x_u
X_fit for unbiased track fit.
std::string m_trackFitResultArrayName
Belle2::TrackFitResult StoreArray name.
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.
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
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.
Helix parameter class.
Definition Helix.h:48
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
Definition Const.h:660
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
CDCHit UsedCDCHit
Define, use of CDC hits as CDC hits (for symmetry).
This is the Reconstruction Event-Data Model Track.
Definition RecoTrack.h:79
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...
Definition RecoTrack.h:621
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.
Definition CDCTrack.h:41
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.
Definition CDCWire.h:58
Vector2D getWirePos2DAtZ(const double z) const
Gives the xy projected position of the wire at the given z coordinate.
Definition CDCWire.h:192
MayBePtr< const CDCWire > getNeighborCCW() const
Gives the closest neighbor in the counterclockwise direction - always exists.
Definition CDCWire.cc:159
IWire getIWire() const
Getter for the wire id within its layer.
Definition CDCWire.h:146
unsigned short getEWire() const
Getter for the encoded wire number.
Definition CDCWire.h:131
ILayer getICLayer() const
Getter for the continuous layer id ranging from 0 - 55.
Definition CDCWire.h:150
MayBePtr< const CDCWire > getNeighborCW() const
Gives the closest neighbor in the clockwise direction - always exists.
Definition CDCWire.cc:164
double norm() const
Calculates the length of the vector.
Definition Vector2D.h:175
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.
Definition Track.h:25
const TrackFitResult * getTrackFitResultWithClosestMass(const Const::ChargedStable &requestedType) const
Return the track fit for the fit hypothesis with the closest mass.
Definition Track.cc:104
Class to identify a wire inside the CDC.
Definition WireID.h:34
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition B2Vector3.h:516
HepGeom::Vector3D< double > Vector3D
3D Vector
Definition Cell.h:34
Abstract base class for different kinds of events.