15 #include <svd/modules/svdPerformance/SVDPerformanceTTreeModule.h>
16 #include <framework/datastore/StoreArray.h>
17 #include <framework/gearbox/Unit.h>
18 #include <framework/core/Environment.h>
19 #include <svd/dataobjects/SVDCluster.h>
20 #include <svd/dataobjects/SVDRecoDigit.h>
21 #include <svd/dataobjects/SVDTrueHit.h>
22 #include <svd/dataobjects/SVDEventInfo.h>
23 #include <vxd/dataobjects/VxdID.h>
24 #include <vxd/geometry/GeoCache.h>
25 #include <vxd/geometry/SensorInfoBase.h>
26 #include <tracking/dataobjects/RecoTrack.h>
27 #include <tracking/dataobjects/RecoHitInformation.h>
28 #include <genfit/TrackPoint.h>
32 #include <mdst/dataobjects/Track.h>
33 #include <mdst/dataobjects/TrackFitResult.h>
52 setDescription(
"The module is used to create a TTree to study SVD clusters, genfit unbiased residuals and many other properties related to the track they belong to.");
54 addParam(
"outputFileName", m_rootFileName,
"Name of output root file.", std::string(
"SVDPerformanceTTree.root"));
55 addParam(
"recoTracksStoreArrayName", m_recoTracksStoreArrayName,
"StoreArray name of the input and output RecoTracks.",
56 m_recoTracksStoreArrayName);
62 recoTracks.isOptional();
64 m_rootFilePtr =
new TFile(m_rootFileName.c_str(),
"RECREATE");
67 m_t_U =
new TTree(
"t_U",
"Tree for SVD u-clusters");
68 m_t_U->Branch(
"svdClSNR", &m_svdClSNR,
"svdClSNR/F");
69 m_t_U->Branch(
"svdClCharge", &m_svdClCharge,
"svdClCharge/F");
70 m_t_U->Branch(
"svdStripCharge", &m_svdStripCharge);
71 m_t_U->Branch(
"svdClTime", &m_svdClTime,
"svdClTime/F");
72 m_t_U->Branch(
"svdStripTime", &m_svdStripTime);
73 m_t_U->Branch(
"svdRes", &m_svdRes,
"svdRes/F");
74 m_t_U->Branch(
"svdClIntStrPos", &m_svdClIntStrPos,
"svdClIntStrPos/F");
75 m_t_U->Branch(
"svdClPos", &m_svdClPos,
"svdClPos/F");
76 m_t_U->Branch(
"svdClPosErr", &m_svdClPosErr,
"svdClPosErr/F");
77 m_t_U->Branch(
"svdTruePos", &m_svdTruePos,
"svdTruePos/F");
78 m_t_U->Branch(
"svdClPhi", &m_svdClPhi,
"svdClPhi/F");
79 m_t_U->Branch(
"svdClZ", &m_svdClZ,
"svdClZ/F");
80 m_t_U->Branch(
"svdTrkd0", &m_svdTrkd0,
"svdTrkd0/F");
81 m_t_U->Branch(
"svdTrkz0", &m_svdTrkz0,
"svdTrkz0/F");
82 m_t_U->Branch(
"svdTrkpT", &m_svdTrkpT,
"svdTrkpT/F");
83 m_t_U->Branch(
"svdTrkpCM", &m_svdTrkpCM,
"svdTrkpCM/F");
84 m_t_U->Branch(
"svdTrkTraversedLength", &m_svdTrkTraversedLength,
"svdTrkTraversedLength/F");
85 m_t_U->Branch(
"svdTrkPXDHits", &m_svdTrkPXDHits,
"svdTrkPXDHits/i");
86 m_t_U->Branch(
"svdTrkSVDHits", &m_svdTrkSVDHits,
"svdTrkSVDHits/i");
87 m_t_U->Branch(
"svdTrkCDCHits", &m_svdTrkCDCHits,
"svdTrkCDCHits/i");
88 m_t_U->Branch(
"svdTrkPos", &m_svdTrkPos,
"svdTrkPos/F");
89 m_t_U->Branch(
"svdTrkPosOS", &m_svdTrkPosOS,
"svdTrkPosOS/F");
90 m_t_U->Branch(
"svdTrkPosErr", &m_svdTrkPosErr,
"svdTrkPosErr/F");
91 m_t_U->Branch(
"svdTrkPosErrOS", &m_svdTrkPosErrOS,
"svdTrkPosErrOS/F");
92 m_t_U->Branch(
"svdTrkQoP", &m_svdTrkQoP,
"svdTrkQoP/F");
93 m_t_U->Branch(
"svdTrkPrime", &m_svdTrkPrime,
"svdTrkPrime/F");
94 m_t_U->Branch(
"svdTrkPrimeOS", &m_svdTrkPrimeOS,
"svdTrkPrimeOS/F");
95 m_t_U->Branch(
"svdTrkPosUnbiased", &m_svdTrkPosUnbiased,
"svdTrkPosUnbiased/F");
96 m_t_U->Branch(
"svdTrkPosErrUnbiased", &m_svdTrkPosErrUnbiased,
"svdTrkPosErrUnbiased/F");
97 m_t_U->Branch(
"svdTrkQoPUnbiased", &m_svdTrkQoPUnbiased,
"svdTrkQoPUnbiased/F");
98 m_t_U->Branch(
"svdTrkPrimeUnbiased", &m_svdTrkPrimeUnbiased,
"svdTrkPrimeUnbiased/F");
99 m_t_U->Branch(
"svdLayer", &m_svdLayer,
"svdLayer/i");
100 m_t_U->Branch(
"svdLadder", &m_svdLadder,
"svdLadder/i");
101 m_t_U->Branch(
"svdSensor", &m_svdSensor,
"svdSensor/i");
102 m_t_U->Branch(
"svdSize", &m_svdSize,
"svdSize/i");
103 m_t_U->Branch(
"svdTB", &m_svdTB,
"svdTB/i");
105 m_t_V =
new TTree(
"t_V",
"Tree for M_SVD v-clusters");
106 m_t_V->Branch(
"svdClSNR", &m_svdClSNR,
"svdClSNR/F");
107 m_t_V->Branch(
"svdClCharge", &m_svdClCharge,
"svdClCharge/F");
108 m_t_V->Branch(
"svdStripCharge", &m_svdStripCharge);
109 m_t_V->Branch(
"svdClTime", &m_svdClTime,
"svdClTime/F");
110 m_t_V->Branch(
"svdStripTime", &m_svdStripTime);
111 m_t_V->Branch(
"svdRes", &m_svdRes,
"svdRes/F");
112 m_t_V->Branch(
"svdClIntStrPos", &m_svdClIntStrPos,
"svdClIntStrPos/F");
113 m_t_V->Branch(
"svdClPos", &m_svdClPos,
"svdClPos/F");
114 m_t_V->Branch(
"svdClPosErr", &m_svdClPosErr,
"svdClPosErr/F");
115 m_t_V->Branch(
"svdTruePos", &m_svdTruePos,
"svdTruePos/F");
116 m_t_V->Branch(
"svdClPhi", &m_svdClPhi,
"svdClPhi/F");
117 m_t_V->Branch(
"svdClZ", &m_svdClZ,
"svdClZ/F");
118 m_t_V->Branch(
"svdTrkd0", &m_svdTrkd0,
"svdTrkd0/F");
119 m_t_V->Branch(
"svdTrkz0", &m_svdTrkz0,
"svdTrkz0/F");
120 m_t_V->Branch(
"svdTrkpT", &m_svdTrkpT,
"svdTrkpT/F");
121 m_t_V->Branch(
"svdTrkpCM", &m_svdTrkpCM,
"svdTrkpCM/F");
122 m_t_V->Branch(
"svdTrkTraversedLength", &m_svdTrkTraversedLength,
"svdTrkTraversedLength/F");
123 m_t_V->Branch(
"svdTrkPXDHits", &m_svdTrkPXDHits,
"svdTrkPXDHits/i");
124 m_t_V->Branch(
"svdTrkSVDHits", &m_svdTrkSVDHits,
"svdTrkSVDHits/i");
125 m_t_V->Branch(
"svdTrkCDCHits", &m_svdTrkCDCHits,
"svdTrkCDCHits/i");
126 m_t_V->Branch(
"svdTrkPos", &m_svdTrkPos,
"svdTrkPos/F");
127 m_t_V->Branch(
"svdTrkPosOS", &m_svdTrkPosOS,
"svdTrkPosOS/F");
128 m_t_V->Branch(
"svdTrkPosErr", &m_svdTrkPosErr,
"svdTrkPosErr/F");
129 m_t_V->Branch(
"svdTrkPosErrOS", &m_svdTrkPosErrOS,
"svdTrkPosErrOS/F");
130 m_t_V->Branch(
"svdTrkQoP", &m_svdTrkQoP,
"svdTrkQoP/F");
131 m_t_V->Branch(
"svdTrkPrime", &m_svdTrkPrime,
"svdTrkPrime/F");
132 m_t_V->Branch(
"svdTrkPrimeOS", &m_svdTrkPrimeOS,
"svdTrkPrimeOS/F");
133 m_t_V->Branch(
"svdTrkPosUnbiased", &m_svdTrkPosUnbiased,
"svdTrkPosUnbiased/F");
134 m_t_V->Branch(
"svdTrkPosErrUnbiased", &m_svdTrkPosErrUnbiased,
"svdTrkPosErrUnbiased/F");
135 m_t_V->Branch(
"svdTrkQoPUnbiased", &m_svdTrkQoPUnbiased,
"svdTrkQoPUnbiased/F");
136 m_t_V->Branch(
"svdTrkPrimeUnbiased", &m_svdTrkPrimeUnbiased,
"svdTrkPrimeUnbiased/F");
137 m_t_V->Branch(
"svdLayer", &m_svdLayer,
"svdLayer/i");
138 m_t_V->Branch(
"svdLadder", &m_svdLadder,
"svdLadder/i");
139 m_t_V->Branch(
"svdSensor", &m_svdSensor,
"svdSensor/i");
140 m_t_V->Branch(
"svdSize", &m_svdSize,
"svdSize/i");
141 m_t_V->Branch(
"svdTB", &m_svdTB,
"svdTB/i");
150 std::string m_svdEventInfoName =
"SVDEventInfo";
152 m_svdEventInfoName =
"SVDEventInfoSim";
154 if (!eventinfo) B2ERROR(
"No SVDEventInfo!");
155 m_svdTB = eventinfo->getModeByte().getTriggerBin();
161 for (
const auto& trk : recoTracks) {
162 if (! trk.wasFitSuccessful()) {
169 if (theTK.
size() == 0) {
176 m_svdTrkd0 = tfr->
getD0();
177 m_svdTrkz0 = tfr->
getZ0();
180 pStar.Boost(0, 0, 3. / 11);
181 m_svdTrkpCM = pStar.P();
185 const vector<SVDCluster* > svdClusters = trk.getSVDHitList();
186 B2DEBUG(40,
"FITTED TRACK: NUMBER OF SVD HITS = " << svdClusters.size());
188 m_svdTrkPXDHits = (trk.getPXDHitList()).size();
189 m_svdTrkSVDHits = (trk.getSVDHitList()).size();
190 m_svdTrkCDCHits = (trk.getCDCHitList()).size();
192 for (
unsigned int i = 0; i < svdClusters.size(); i++) {
203 const auto* hitTrackPoint_1 = trk.getCreatedTrackPoint(infoSVD_1);
204 const auto* fittedResult_1 = hitTrackPoint_1->getFitterInfo();
205 if (!fittedResult_1) {
214 const TVectorD resUnBias_1 = fittedResult_1->getResidual(0,
false).getState();
216 const TVectorD& svd_predIntersect_unbiased = state_unbiased.getState();
217 const TMatrixDSym& covMatrix_unbiased = state_unbiased.getCov();
219 const TVectorD& svd_predIntersect_1 = state_1.getState();
220 const TMatrixDSym& covMatrix_1 = state_1.getCov();
224 const int strips_1 = svd_1->
getSize();
227 const TVector3 svdLocal_1(svd_1->
getPosition(), svd_predIntersect_1[4], 0.);
229 const TVector3& svdGlobal_1 = svdSensor_1.
pointToGlobal(svdLocal_1);
230 double svdPhi_1 = atan2(svdGlobal_1(1), svdGlobal_1(0));
231 double svdZ_1 = svdGlobal_1(2);
236 m_svdClSNR = svd_1->
getSNR();
240 if (isMC && trueHit_1.
size() > 0)
241 m_svdTruePos = trueHit_1[0]->getU();
244 m_svdClPhi = svdPhi_1;
246 m_svdTrkPos = svd_predIntersect_1[3];
247 m_svdTrkPosOS = svd_predIntersect_1[4];
248 m_svdTrkPosErr = sqrt(covMatrix_1[3][3]);
249 m_svdTrkPosErrOS = sqrt(covMatrix_1[4][4]);
250 m_svdTrkQoP = svd_predIntersect_1[0];
251 m_svdTrkPrime = svd_predIntersect_1[1];
252 m_svdTrkPrimeOS = svd_predIntersect_1[2];
253 m_svdTrkTraversedLength = svdSensor_1.
getThickness() * sqrt(1 + m_svdTrkPrimeOS * m_svdTrkPrimeOS + m_svdTrkPrime * m_svdTrkPrime);
254 m_svdTrkPosUnbiased = svd_predIntersect_unbiased[3];
255 m_svdTrkPosErrUnbiased = sqrt(covMatrix_unbiased[3][3]);
256 m_svdTrkQoPUnbiased = svd_predIntersect_unbiased[0];
257 m_svdTrkPrimeUnbiased = svd_predIntersect_unbiased[1];
258 m_svdLayer = svd_Layer_1;
259 m_svdLadder = svd_Ladder_1;
260 m_svdSensor = svd_Sensor_1;
261 m_svdSize = strips_1;
264 float halfLength = 1.92;
265 if (m_svdLayer > 3) {
269 m_svdClIntStrPos = fmod(m_svdClPos + halfLength, pitch) / pitch;
271 m_svdStripCharge.clear();
272 m_svdStripTime.clear();
275 if ((theRecoDigits.
size() != m_svdSize) && (m_svdSize != 128))
276 B2ERROR(
" Inconsistency with cluster size! # recoDigits = " << theRecoDigits.
size() <<
" != " << m_svdSize <<
" cluster size");
278 if (!(m_svdSize == 128 && theRecoDigits.
size() == 0))
279 for (
unsigned int d = 0; d < m_svdSize; d++) {
280 m_svdStripCharge.push_back(theRecoDigits[d]->getCharge());
281 m_svdStripTime.push_back(theRecoDigits[d]->getTime());
287 const int strips_1 = svd_1->
getSize();
289 const TVector3 svdLocal_1(svd_predIntersect_1[3], svd_1->
getPosition(), 0.);
291 const TVector3& svdGlobal_1 = svdSensor_1.
pointToGlobal(svdLocal_1);
292 double svdPhi_1 = atan2(svdGlobal_1(1), svdGlobal_1(0));
293 double svdZ_1 = svdGlobal_1(2);
297 m_svdClSNR = svd_1->
getSNR();
301 if (isMC && trueHit_1.
size() > 0)
302 m_svdTruePos = trueHit_1[0]->getV();
305 m_svdClPhi = svdPhi_1;
307 m_svdTrkPos = svd_predIntersect_1[4];
308 m_svdTrkPosOS = svd_predIntersect_1[3];
309 m_svdTrkPosErr = sqrt(covMatrix_1[4][4]);
310 m_svdTrkPosErrOS = sqrt(covMatrix_1[3][3]);
311 m_svdTrkQoP = svd_predIntersect_1[0];
312 m_svdTrkPrime = svd_predIntersect_1[2];
313 m_svdTrkPrimeOS = svd_predIntersect_1[1];
314 m_svdTrkTraversedLength = svdSensor_1.
getThickness() * sqrt(1 + m_svdTrkPrimeOS * m_svdTrkPrimeOS + m_svdTrkPrime * m_svdTrkPrime);
315 m_svdTrkPosUnbiased = svd_predIntersect_unbiased[4];
316 m_svdTrkPosErrUnbiased = sqrt(covMatrix_unbiased[4][4]);
317 m_svdTrkQoPUnbiased = svd_predIntersect_unbiased[0];
318 m_svdTrkPrimeUnbiased = svd_predIntersect_unbiased[2];
319 m_svdLayer = svd_Layer_1;
320 m_svdLadder = svd_Ladder_1;
321 m_svdSensor = svd_Sensor_1;
322 m_svdSize = strips_1;
324 float pitch = 160e-4;
325 float halfLength = 6.144;
326 if (m_svdLayer > 3) {
329 m_svdClIntStrPos = fmod(m_svdClPos + halfLength, pitch) / pitch;
331 m_svdStripCharge.clear();
332 m_svdStripTime.clear();
337 if ((theRecoDigits.
size() != m_svdSize) && (m_svdSize != 128))
338 B2ERROR(
" Inconsistency with cluster size! # recoDigits = " << theRecoDigits.
size() <<
" != " << m_svdSize <<
" cluster size");
340 if (!(m_svdSize == 128 && theRecoDigits.
size() == 0))
341 for (
unsigned int d = 0; d < m_svdSize; d++) {
342 m_svdStripCharge.push_back(theRecoDigits[d]->getCharge());
343 m_svdStripTime.push_back(theRecoDigits[d]->getTime());
349 B2INFO(
"oops...something went wrong in getting the unbiased state, skipping this cluster.");
361 if (m_rootFilePtr != NULL) {
365 m_rootFilePtr->Close();