Belle II Software development
CDCT0CalibrationCollector.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/CDCT0CalibrationCollector.h"
10
11#include <framework/datastore/RelationArray.h>
12
13#include <genfit/TrackPoint.h>
14#include <genfit/KalmanFitterInfo.h>
15#include <genfit/MeasuredStateOnPlane.h>
16#include <cdc/dataobjects/WireID.h>
17#include <cdc/geometry/CDCGeometryPar.h>
18
19#include <Math/ProbFuncMathCore.h>
20#include <TH1F.h>
21
22using namespace std;
23using namespace Belle2;
24using namespace CDC;
25using namespace genfit;
26
27REG_MODULE(CDCT0CalibrationCollector);
28
30{
31 setDescription("Collector module for cdc T0 calibration");
33 addParam("RecoTracksColName", m_recoTrackArrayName, "Name of collection hold genfit::Track", std::string(""));
34 addParam("BField", m_BField, "If true -> #Params ==5 else #params ==4 for calculate P-Val", false);
35 addParam("EventT0Extraction", m_EventT0Extraction, "use event t0 extract t0 or not", false);
36 addParam("MinimumPt", m_MinimumPt, "Tracks whose Pt lower than this value will not be recoreded", 0.);
37 addParam("PvalCut", m_PvalCut, "Tracks which Pval small than this will not be recoreded", 0.);
38 addParam("NDFCut", m_ndfCut, "Tracks which NDF small than this will not be recoreded", 0.);
39 addParam("Xmin", m_xmin, "minimum dift length", 0.1);
40}
41
45
47{
48 m_Tracks.isRequired(m_trackArrayName);
51 m_CDCHits.isRequired(m_cdcHitArrayName);
53 //Store names to speed up creation later
54 m_relRecoTrackTrackName = relRecoTrackTrack.getName();
55
56 auto m_hNDF = new TH1D("hNDF", "NDF of fitted track;NDF;Tracks", 71, -1, 70);
57 auto m_hPval = new TH1D("hPval", "p-values of tracks;pVal;Tracks", 1000, 0, 1);
58
59 registerObject<TH1D>("hNDF", m_hNDF);
60 registerObject<TH1D>("hPval", m_hPval);
61
63 for (int i = 0; i < 56; ++i) {
64 double R = cdcgeo.senseWireR(i);
65 double nW = cdcgeo.nWiresInLayer(i);
66 halfCSize[i] = M_PI * R / nW;
67 }
68
69 auto m_hTotal = new TH1F("hTotal", "hTotal", 30, -15, 15);
70 registerObject<TH1F>("hTotal", m_hTotal);
71 //for each channel
72 TH1F* m_h1[56][385];
73 TH1F* m_hT0b[300];
74 for (int il = 0; il < 56; il++) {
75 const int nW = cdcgeo.nWiresInLayer(il);
76 for (int ic = 0; ic < nW; ++ic) {
77 m_h1[il][ic] = new TH1F(Form("hdT_lay%d_cell%d", il, ic), Form("Lay%d_cell%d", il, ic), 30, -15, 15);
78 registerObject<TH1F>(Form("hdT_lay%d_cell%d", il, ic), m_h1[il][ic]);
79 }
80 }
81 //for each board
82 for (int ib = 0; ib < 300; ++ib) {
83 m_hT0b[ib] = new TH1F(Form("hT0b%d", ib), Form("boardID_%d", ib), 100, -20, 20);
84 registerObject<TH1F>(Form("hT0b%d", ib), m_hT0b[ib]);
85 }
86
87}
88
90{
92
93 /* CDCHit distribution */
94 // make evt t0 in case we don't use evt t0
95 double evtT0 = 0;
96 const int nTr = m_RecoTracks.getEntries();
97
98 for (int i = 0; i < nTr; ++i) {
99 RecoTrack* track = m_RecoTracks[i];
100 if (track->getDirtyFlag()) continue;
101 const genfit::FitStatus* fs = track->getTrackFitStatus();
102 if (!fs || !fs->isFitted()) continue;
103 if (!fs->isFitConverged()) continue; //not fully converged
104
105 const Belle2::Track* b2track = track->getRelatedFrom<Belle2::Track>();
106 if (!b2track) {
107 B2DEBUG(99, "No relation found");
108 continue;
109 }
110
111 const Belle2::TrackFitResult* fitresult = b2track->getTrackFitResult(Const::muon);
112
113 if (!fitresult) {
114 B2WARNING("track was fitted but Relation not found");
115 continue;
116 }
117
118 double ndf = 0;
119 double Pval = 0;
120 if (m_BField) {
121 ndf = fs->getNdf();
122 Pval = fs->getPVal();
123 } else {
124 ndf = fs->getNdf() + 1;
125 double Chi2 = fs->getChi2();
126 Pval = std::max(0., ROOT::Math::chisquared_cdf_c(Chi2, ndf));
127 }
128
129 getObjectPtr<TH1D>("hPval")->Fill(Pval);
130 getObjectPtr<TH1D>("hNDF")->Fill(ndf);
131
132 if (Pval < m_PvalCut || ndf < m_ndfCut) continue;
133 //cut at Pt
134 if (fitresult->getMomentum().Rho() < m_MinimumPt) continue;
135 //reject events don't have eventT0
137 // event with is fail to extract t0 will be exclude from analysis
138 if (m_eventTimeStoreObject.isValid() && m_eventTimeStoreObject->hasEventT0()) {
139 evtT0 = m_eventTimeStoreObject->getEventT0();
140 } else {
141 continue;
142 }
143 }
144
145 B2DEBUG(99, "start collect hit");
146 static CDCGeometryPar& cdcgeo = CDCGeometryPar::Instance();
147
148 for (const RecoHitInformation::UsedCDCHit* hit : track->getCDCHitList()) {
149 const genfit::TrackPoint* tp = track->getCreatedTrackPoint(track->getRecoHitInformation(hit));
150 int lay = hit->getICLayer();
151 int IWire = hit->getIWire();
152 unsigned short tdc = hit->getTDCCount();
153 unsigned short adc = hit->getADCCount();
154 WireID wireid(lay, IWire);
155 const genfit::KalmanFitterInfo* kfi = tp->getKalmanFitterInfo();
156 if (!kfi) {B2DEBUG(199, "No Fitter Info: Layer " << hit->getICLayer()); continue;}
157 for (unsigned int iMeas = 0; iMeas < kfi->getNumMeasurements(); ++iMeas) {
158 if ((kfi->getWeights().at(iMeas)) > 0.5) {
159 int boardID = cdcgeo.getBoardID(WireID(lay, IWire));
160 const genfit::MeasuredStateOnPlane& mop = kfi->getFittedState();
161 const B2Vector3D pocaOnWire = mop.getPlane()->getO();//Local wire position
162 // const B2Vector3D pocaOnTrack = mop.getPlane()->getU();//residual direction
163 const B2Vector3D pocaMom = mop.getMom();
164 double alpha = cdcgeo.getAlpha(pocaOnWire, pocaMom) ;
165 double theta = cdcgeo.getTheta(pocaMom);
166 double x_u = kfi->getFittedState(false).getState()(3);//x fit unbiased
167 // double weight = kfi->getWeights().at(iMeas);
168
169 double xmax = halfCSize[lay] - 0.1;
170 if ((fabs(x_u) < m_xmin) || (fabs(x_u) > xmax)) continue;
171
172 int lr;
173 if (x_u > 0) lr = 1;
174 else lr = 0;
175
176 //Convert to outgoing
177 if (fabs(alpha) > M_PI / 2) {
178 x_u *= -1;
179 }
180
181 lr = cdcgeo.getOutgoingLR(lr, alpha);
182 theta = cdcgeo.getOutgoingTheta(alpha, theta);
183 alpha = cdcgeo.getOutgoingAlpha(alpha);
184
185 double t_fit = cdcgeo.getDriftTime(std::abs(x_u), lay, lr, alpha, theta);
186 //estimate drift time
187 double dt_flight;
188 double dt_prop;
189 double t = cdcgeo.getT0(wireid) - tdc * cdcgeo.getTdcBinWidth(); // - dt_flight - dt_prop;
190 dt_flight = mop.getTime();
191 if (dt_flight < 50) {
192 t -= dt_flight;
193 } else {
194 continue;
195 }
196
197 double z = pocaOnWire.Z();
198 B2Vector3D m_backWirePos = cdcgeo.wireBackwardPosition(wireid, CDCGeometryPar::c_Aligned);
199 double z_prop = z - m_backWirePos.Z();
200 B2DEBUG(99, "z_prop = " << z_prop << " |z " << z << " |back wire poss: " << m_backWirePos.Z());
201 dt_prop = z_prop * cdcgeo.getPropSpeedInv(lay);
202 if (z_prop < 240) {
203 t -= dt_prop;
204 } else {
205 continue;
206 }
207 // Time Walk
208 t -= cdcgeo.getTimeWalk(wireid, adc);
209 // subtract event t0;
210 t -= evtT0;
211
212 getObjectPtr<TH1F>(Form("hdT_lay%d_cell%d", lay, IWire))->Fill(t - t_fit);
213 getObjectPtr<TH1F>(Form("hT0b%d", boardID))->Fill(t - t_fit);
214 getObjectPtr<TH1F>("hTotal")->Fill(t - t_fit);
215 } //NDF
216 // }//end of if isU
217 }//end of for
218 }//end of for tp
219 }//end of func
220}
221
225
226
double R
typedef autogenerated by FFTW
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition B2Vector3.h:435
The Class for CDC Geometry Parameters.
double getTimeWalk(const WireID &wID, unsigned short adcCount) const
Returns time-walk.
double getTheta(const B2Vector3D &momentum) const
Returns track incident angle (theta in rad.).
unsigned short getBoardID(const WireID &wID) const
Returns frontend board id. corresponding to the wire id.
double getAlpha(const B2Vector3D &posOnWire, const B2Vector3D &momentum) const
Returns track incident angle in rphi plane (alpha in rad.).
const B2Vector3D wireBackwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the backward position of the input sense wire.
double getTdcBinWidth() const
Return TDC bin width (nsec).
double getOutgoingAlpha(const double alpha) const
Converts incoming- to outgoing-alpha.
float getT0(const WireID &wireID) const
Returns t0 parameter of the specified sense wire.
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.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
double getPropSpeedInv(const unsigned int layerID) const
Get the inversel of propagation speed in the sense wire.
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.
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
bool m_BField
fit in case no magnetic field or not, if false, NDF=4 in cal P-value
StoreObjPtr< EventT0 > m_eventTimeStoreObject
Event t0 object.
std::string m_recoTrackArrayName
Belle2::RecoTrack StoreArray name.e.
double m_MinimumPt
minimum pt required for track
StoreArray< TrackFitResult > m_TrackFitResults
Track fit results.
std::string m_cdcHitArrayName
Belle2::CDCHit StoreArray name.
std::string m_relRecoTrackTrackName
Relation between RecoTrack and Belle2:Track.
void collect() override
Event action, collect information for calibration.
std::string m_trackArrayName
Belle2::Track StoreArray name.
void prepare() override
Initializes the Module.
std::string m_trackFitResultArrayName
Belle2::TrackFitResult StoreArray name.
bool m_EventT0Extraction
use Event T0 extraction or not.
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
Low-level class to create/modify relations between StoreArrays.
const std::string & getName() const
Return name under which the object is saved in the DataStore.
Values of the result of a track fit with a given particle hypothesis.
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Class that bundles various TrackFitResults.
Definition Track.h:25
const TrackFitResult * getTrackFitResult(const Const::ChargedStable &chargedStable) const
Default Access to TrackFitResults.
Definition Track.cc:30
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
Abstract base class for different kinds of events.
STL namespace.