11 #include "cdc/modules/cdcCalibrationCollector/CDCCalibrationCollector.h"
12 #include <cdc/translators/RealisticTDCCountTranslator.h>
13 #include <framework/datastore/StoreObjPtr.h>
14 #include <framework/datastore/StoreArray.h>
15 #include <framework/datastore/RelationArray.h>
17 #include <mdst/dataobjects/TrackFitResult.h>
18 #include <mdst/dataobjects/Track.h>
19 #include <tracking/dataobjects/RecoTrack.h>
20 #include <tracking/trackFindingCDC/rootification/StoreWrappedObjPtr.h>
21 #include <tracking/trackFindingCDC/eventdata/hits/CDCWireHit.h>
22 #include <tracking/trackFindingCDC/eventdata/tracks/CDCTrack.h>
23 #include <tracking/trackFindingCDC/topology/CDCWireTopology.h>
25 #include <genfit/TrackPoint.h>
26 #include <genfit/KalmanFitterInfo.h>
27 #include <genfit/MeasurementOnPlane.h>
28 #include <genfit/MeasuredStateOnPlane.h>
30 #include <Math/ProbFuncMathCore.h>
32 #include <cdc/dataobjects/WireID.h>
33 #include <cdc/geometry/CDCGeometryPar.h>
41 using namespace TrackFindingCDC;
49 setDescription(
"Collector module for cdc calibration");
50 setPropertyFlags(c_ParallelProcessingCertified);
51 addParam(
"recoTracksColName", m_recoTrackArrayName,
"Name of collection hold genfit::Track", std::string(
""));
52 addParam(
"bField", m_bField,
"If true -> #Params ==5 else #params ==4 for calculate P-Val",
false);
53 addParam(
"calExpectedDriftTime", m_calExpectedDriftTime,
"if true module will calculate expected drift time, it take a time",
55 addParam(
"storeTrackParams", m_storeTrackParams,
"Store Track Parameter or not, it will be multicount for each hit",
false);
56 addParam(
"eventT0Extraction", m_eventT0Extraction,
"use event t0 extract t0 or not",
true);
57 addParam(
"minimumPt", m_minimumPt,
"Tracks with tranverse momentum smaller than this value will not used", 0.15);
58 addParam(
"isCosmic", m_isCosmic,
"True when we process cosmic events, else False (collision)", m_isCosmic);
59 addParam(
"effStudy", m_effStudy,
"When true module collects info only necessary for wire eff study",
false);
73 RelationArray relRecoTrackTrack(recoTracks, storeTrack, m_relRecoTrackTrackName);
75 m_recoTrackArrayName = recoTracks.getName();
76 m_trackFitResultArrayName = storeTrackFitResults.getName();
77 m_relRecoTrackTrackName = relRecoTrackTrack.
getName();
80 auto m_tree =
new TTree(m_treeName.c_str(),
"tree for cdc calibration");
81 m_tree->Branch<Float_t>(
"x_mea", &x_mea);
82 m_tree->Branch<Float_t>(
"x_u", &x_u);
83 m_tree->Branch<Float_t>(
"x_b", &x_b);
84 m_tree->Branch<Float_t>(
"alpha", &alpha);
85 m_tree->Branch<Float_t>(
"theta", &theta);
86 m_tree->Branch<Float_t>(
"t", &t);
87 m_tree->Branch<UShort_t>(
"adc", &adc);
89 m_tree->Branch<UChar_t>(
"lay", &lay);
90 m_tree->Branch<Float_t>(
"weight", &weight);
91 m_tree->Branch<UShort_t>(
"IWire", &IWire);
92 m_tree->Branch<Float_t>(
"Pval", &Pval);
93 m_tree->Branch<Float_t>(
"ndf", &ndf);
94 if (m_storeTrackParams) {
95 m_tree->Branch<Float_t>(
"d0", &d0);
96 m_tree->Branch<Float_t>(
"z0", &z0);
97 m_tree->Branch<Float_t>(
"phi0", &phi0);
98 m_tree->Branch<Float_t>(
"tanL", &tanL);
99 m_tree->Branch<Float_t>(
"omega", &omega);
102 if (m_calExpectedDriftTime) {
103 m_tree->Branch<Float_t>(
"t_fit", &t_fit);
106 registerObject<TTree>(
"tree", m_tree);
109 auto m_efftree =
new TTree(m_effTreeName.c_str(),
"tree for wire efficiency");
110 m_efftree->Branch<
unsigned short>(
"layerID", &layerID);
111 m_efftree->Branch<
unsigned short>(
"wireID", &wireID);
112 m_efftree->Branch<
float>(
"z", &z);
113 m_efftree->Branch<
bool>(
"isFound", &isFound);
115 registerObject<TTree>(
"efftree", m_efftree);
118 auto m_hNDF =
new TH1F(
"hNDF",
"NDF of fitted track;NDF;Tracks", 71, -1, 70);
119 auto m_hPval =
new TH1F(
"hPval",
"p-values of tracks;pVal;Tracks", 1000, 0, 1);
120 auto m_hEventT0 =
new TH1F(
"hEventT0",
"Event T0", 1000, -100, 100);
121 auto m_hNTracks =
new TH1F(
"hNTracks",
"Number of tracks", 50, 0, 10);
122 auto m_hOccupancy =
new TH1F(
"hOccupancy",
"occupancy", 100, 0, 1.0);
124 registerObject<TH1F>(
"hNDF", m_hNDF);
125 registerObject<TH1F>(
"hPval", m_hPval);
126 registerObject<TH1F>(
"hEventT0", m_hEventT0);
127 registerObject<TH1F>(
"hNTracks", m_hNTracks);
128 registerObject<TH1F>(
"hOccupancy", m_hOccupancy);
138 const RelationArray relTrackTrack(recoTracks, storeTrack, m_relRecoTrackTrackName);
145 if (m_eventT0Extraction) {
147 if (m_eventTimeStoreObject.isValid() && m_eventTimeStoreObject->hasEventT0()) {
148 evtT0 = m_eventTimeStoreObject->getEventT0();
149 getObjectPtr<TH1F>(
"hEventT0")->Fill(evtT0);
156 std::vector<unsigned short> wiresInCDCTrack;
158 for (CDCTrack& cdcTrack : *cdcTracks) {
159 for (CDCRecoHit3D& cdcHit : cdcTrack) {
160 unsigned short eWireID = cdcHit.getWire().getEWire();
161 wiresInCDCTrack.push_back(eWireID);
170 for (
int i = 0; i < nTr; ++i) {
172 if (track->getDirtyFlag())
continue;
173 if (!track->getTrackFitStatus()->isFitted())
continue;
178 if (!b2track) {B2DEBUG(99,
"No relation found");
continue;}
181 B2WARNING(
"track was fitted but Relation not found");
186 if (fabs(charge) > 0) {
193 getObjectPtr<TH1F>(
"hNTracks")->Fill(nCTracks);
198 const int nWires = 14336;
199 float oc =
static_cast<float>(nHits) /
static_cast<float>(nWires);
200 getObjectPtr<TH1F>(
"hOccupancy")->Fill(oc);
202 for (
int i = 0; i < nTr; ++i) {
204 if (track->getDirtyFlag())
continue;
206 if (!fs || !fs->
isFitted())
continue;
210 if (!b2track) {B2DEBUG(99,
"No relation found");
continue;}
214 B2WARNING(
"track was fitted but Relation not found");
222 getObjectPtr<TH1F>(
"hPval")->Fill(Pval);
223 getObjectPtr<TH1F>(
"hNDF")->Fill(ndf);
224 B2DEBUG(99,
"ndf = " << ndf);
225 B2DEBUG(99,
"Pval = " << Pval);
227 if (ndf < 15)
continue;
229 Pval = std::max(0., ROOT::Math::chisquared_cdf_c(Chi2, ndf));
232 d0 = fitresult->
getD0();
233 z0 = fitresult->
getZ0();
236 phi0 = fitresult->
getPhi0() * 180 / M_PI;
240 if (m_isCosmic ==
true && phi0 > 0.0)
continue;
252 if (fitresult->
getD0() > 2 || fitresult->
getZ0() > 5)
continue;
254 buildEfficiencies(wiresInCDCTrack, helixFit);
266 B2DEBUG(99,
"start collect hit");
271 const genfit::TrackPoint* tp = track->getCreatedTrackPoint(track->getRecoHitInformation(hit));
273 lay = hit->getICLayer();
274 IWire = hit->getIWire();
275 adc = hit->getADCCount();
276 unsigned short tdc = hit->getTDCCount();
277 WireID wireid(lay, IWire);
279 if (!kfi) {B2DEBUG(199,
"No Fitter Info: Layer " << hit->getICLayer());
continue;}
280 for (
unsigned int iMeas = 0; iMeas < kfi->getNumMeasurements(); ++iMeas) {
284 const TVector3 pocaOnWire = mop.getPlane()->getO();
285 const TVector3 pocaOnTrack = mop.getPlane()->getU();
286 const TVector3 pocaMom = mop.getMom();
287 alpha = cdcgeo.
getAlpha(pocaOnWire, pocaMom) ;
289 x_mea = kfi->getMeasurementOnPlane(iMeas)->getState()(0);
299 if (fabs(alpha) > M_PI / 2) {
303 x_mea = copysign(x_mea, x_b);
308 B2DEBUG(99,
"x_unbiased " << x_u <<
" |left_right " << lr);
309 if (m_calExpectedDriftTime) { t_fit = cdcgeo.
getDriftTime(abs(x_u), lay, lr, alpha , theta);}
313 t = tdcTrans->
getDriftTime(tdc, wireid, mop.getTime(), pocaOnWire.Z(), adc);
314 getObjectPtr<TTree>(
"tree")->Fill();
322 const Helix& helixFit)
const
328 const CDCWire& oneWire = layer.getWire(1);
329 double newR = oneWire.getWirePos2DAtZ(xyz.z()).norm();
330 double arcLength = helixFit.getArcLength2DAtCylindricalR(newR);
331 TVector3 xyzOnWire = helixFit.getPositionAtArcLength2D(arcLength);
335 const CDCWire& wire = layer.getClosestWire(crosspoint);
343 for (
const CDCWireLayer& wireLayer : wireTopology.
getWireLayers()) {
344 const double radiusofLayer = wireLayer.getRefCylindricalR();
346 const double arcLength = helixFit.getArcLength2DAtCylindricalR(radiusofLayer);
347 const TVector3 xyz = helixFit.getPositionAtArcLength2D(arcLength);
348 if (!xyz.x())
continue;
349 const CDCWire& wireIntersected = getIntersectingWire(xyz, wireLayer, helixFit);
350 unsigned short crossedWire = wireIntersected.getEWire();
351 unsigned short crossedCWire = wireIntersected.getNeighborCW()->getEWire();
352 unsigned short crossedCCWire = wireIntersected.getNeighborCCW()->getEWire();
354 if (find(wireHits.begin(), wireHits.end(), crossedWire) != wireHits.end()
355 || find(wireHits.begin(), wireHits.end(), crossedCWire) != wireHits.end()
356 || find(wireHits.begin(), wireHits.end(), crossedCCWire) != wireHits.end())
361 wireID = wireIntersected.getIWire();
362 layerID = wireIntersected.getICLayer();
364 getObjectPtr<TTree>(
"efftree")->Fill();