Belle II Software development
CDCCosmicSelectorAfterFullSimModule.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/cdcCosmicSelectorAfterFullSim/CDCCosmicSelectorAfterFullSimModule.h>
10#include <framework/logging/Logger.h>
11#include <framework/datastore/RelationArray.h>
12#include <framework/datastore/RelationIndex.h>
13
14using namespace std;
15using namespace Belle2;
16
17// register module
18REG_MODULE(CDCCosmicSelectorAfterFullSim);
20{
21 // Set description
22 setDescription("Modify CDCSimHits and MCParticles for cosmics so that the global time is (approximately) zero at y=0 using FullSim output (=CDCSimHits). And select cosmics passing through a user-specified region at y=0. This module works only for the event with the no. of primary charged MC particles=1. Please place this module after FullSim and before CDCDigitizer. This module is not needed for normal MC; only needed when you want to evaluate performance of T0 extraction module.");
24
25 addParam("xOfRegion", m_xOfRegion, "x of center position of region at y=0 (cm)", 0.);
26 addParam("zOfRegion", m_zOfRegion, "z of center position of region at y=0 (cm)", 0.);
27 addParam("wOfRegion", m_wOfRegion, "full-width (x-direction) of region at y=0 (cm)", 10.);
28 addParam("lOfRegion", m_lOfRegion, "full-length (z-direction) of region at y=0 (cm)", 100.);
29}
30
32{
33 m_mcParticles.isRequired();
34 m_simHits.isRequired();
35}
36
37
39{
40 // Get relation array
41 RelationArray mcParticlesToCDCSimHits(m_mcParticles, m_simHits);
42
43 bool returnVal = false;
44
45 // const double c = 29.9792458; //light speed (cm/ns)
46 // double mass = 0.1056583715;
47
48 int nHits = m_simHits.getEntries();
49 if (nHits <= 1) {
50 B2DEBUG(250, "No .of CDCSimHits <= 1 in the current event.");
51 // B2WARNING("No .of CDCSimHits=0 in the current event= " << nHits);
52 setReturnValue(returnVal);
53 return;
54 }
55
56 int nMCPs = m_mcParticles.getEntries();
57 if (nMCPs == 0) {
58 B2ERROR("No. of MCParticles=0 in the current event.");
59 setReturnValue(returnVal);
60 return;
61 }
62
63 // Loop over prim. charged MC particle
64 unsigned nPrimChgds = 0;
65 for (int iMCP = 0; iMCP < nMCPs; ++iMCP) {
66 MCParticle* m_P = m_mcParticles[iMCP];
67
68 // B2INFO("isPrimaryParticle= " << m_P->isPrimaryParticle());
69 if (!m_P->isPrimaryParticle()) continue;
70
71 unsigned pid = abs(m_P->getPDG());
72 if (pid == 13) {
73 ++nPrimChgds;
74 } else if (pid == 11) {
75 ++nPrimChgds;
76 // mass = 0.000510998928;
77 } else if (pid == 211) {
78 ++nPrimChgds;
79 // mass = 0.13957018;
80 } else if (pid == 321) {
81 ++nPrimChgds;
82 // mass = 0.493677;
83 } else if (pid == 2212) {
84 ++nPrimChgds;
85 // mass = 0.938272046;
86 } else {
87 continue;
88 }
89
90 // B2INFO("No .of prim. charged MC particles= " << nPrimChgds);
91 if (nPrimChgds > 1) continue;
92 /*
93 const B2Vector3D vertex = m_P->getProductionVertex();
94 const double vX0 = vertex.X();
95 const double vY0 = vertex.Y();
96 const double vZ0 = vertex.Z();
97 const B2Vector3D momentum = m_P->getMomentum();
98 const double pX0 = momentum.X();
99 const double pY0 = momentum.Y();
100 const double pZ0 = momentum.Z();
101 const double p = momentum.Mag();
102 const double beta = p / sqrt(p * p + mass * mass);
103
104 //calculate an intersection with y=0 plane
105 double xi4cry = -999.;
106 double yi4cry = 0.;
107 double zi4cry = -999.;
108 if (pY0 != 0.) {
109 xi4cry = (0. - vY0) * (pX0 / pY0) + vX0;
110 zi4cry = (0. - vY0) * (pZ0 / pY0) + vZ0;
111 } else {
112 xi4cry = 0.;
113 zi4cry = -vX0 * (pZ0 / pX0) + vZ0;
114 }
115
116 const double fl4cry = sqrt((xi4cry - vX0) * (xi4cry - vX0) + (yi4cry - vY0) * (yi4cry - vY0) + (zi4cry - vZ0) *
117 (zi4cry - vZ0)); // fl to y=0 plane
118 */
119
120 //try to calculate cdc hit with tof_up and tof_dn
121 double tof_up = -9999.;
122 double tof_dn = 9999.;
123 double ihit_up = -1;
124 double ihit_dn = -1;
126 if (!mcp_to_hit) B2FATAL("No MCParticle->CDCSimHit relation found!");
128 // std::cout <<" "<< std::endl;
129
130 unsigned ntry = 0;
131 double yold = 0.;
132 double tofold = 0.;
133 unsigned iHitold = 0 ;
134 bool crossfind = false;
135
136 for (int iHit = 0; iHit < nHits; ++iHit) {
137 if (crossfind) break;
138 for (const RelationElement& rel : mcp_to_hit.getElementsTo(m_simHits[iHit])) {
139 // std::cout <<"iHit,iMCP,rfromindex= " << iHit <<" "<< iMCP <<" "<< rel.from->getIndex()-1 << std::endl;
140 if ((rel.from->getIndex() - 1) != iMCP) continue;
141 // std::cout << "weight,pdginsimhit= " << rel.weight <<" "<< simHits[iHit]->getPDGCode() << std::endl;
142 if (rel.weight < 0.) continue; //reject 2ndary particle
143 const double y = m_simHits[iHit]->getPosTrack().Y();
144 const double tof = m_simHits[iHit]->getFlightTime();
145 // const double py = simHits[iHit]->getMomentum().Y();
146
147 ++ntry;
148
149 if (ntry == 1) {
150 yold = y;
151 tofold = tof;
152 iHitold = iHit;
153 } else {
154 // std::cout <<"yold,y= " << yold<<" "<< y << std::endl;
155 if (y * yold >= 0.) {
156 yold = y;
157 tofold = tof;
158 iHitold = iHit;
159 } else {
160 tof_up = tofold;
161 ihit_up = iHitold;
162 tof_dn = tof;
163 ihit_dn = iHit;
164 crossfind = true;
165 break;
166 }
167 }
168 /*
169 if (y > 0. && py < 0. && tof > tof_up) {
170 tof_up = tof;
171 ihit_up = iHit;
172 }
173 if (y < 0. && py < 0. && tof < tof_dn) {
174 tof_dn = tof;
175 ihit_dn = iHit;
176 }
177 */
178 }
179 }
180
181 //calculate flight time from y_up to y=0 plane (linear approx.)
182 // std::cout <<"ihit_up,dn= " << ihit_up <<" "<< ihit_dn << std::endl;
183 if (ihit_up < 0 || ihit_dn < 0) continue;
184 const B2Vector3D pos_up = m_simHits[ihit_up]->getPosTrack();
185 const B2Vector3D pos_dn = m_simHits[ihit_dn]->getPosTrack();
186 // const B2Vector3D mom_up = m_simHits[ihit_up]->getMomentum();
187 // const B2Vector3D mom_dn = m_simHits[ihit_dn]->getMomentum();
188 if (tof_up > tof_dn) B2WARNING("tof_up > tof_dn " << tof_up << " " << tof_dn);
189 // std::cout <<"tof_up,dn= " << tof_up <<" "<< tof_dn << std::endl;
190
191 // const B2Vector3D mom_mean = 0.5 * (mom_up + mom_dn);
192 const B2Vector3D dpos = pos_dn - pos_up;
193 double cx = dpos.X() / dpos.Mag();
194 double cy = dpos.Y() / dpos.Mag();
195 double cz = dpos.Z() / dpos.Mag();
196 double vx = pos_up.X();
197 double vy = pos_up.Y();
198 double vz = pos_up.Z();
199 if (tof_up > tof_dn) {
200 cx *= -1.;
201 cy *= -1.;
202 cz *= -1.;
203 vx = pos_dn.X();
204 vy = pos_dn.Y();
205 vz = pos_dn.Z();
206 }
207
208 const double yi = 0.;
209 double xi = 0.;
210 double zi = 0.;
211 if (cy != 0.) {
212 xi = (yi - vy) * (cx / cy) + vx;
213 zi = (yi - vy) * (cz / cy) + vz;
214 } else {
215 zi = -vx * (cz / cx) + vz;
216 }
217
218 bool hitRegion = (fabs(xi - m_xOfRegion) < 0.5 * m_wOfRegion && fabs(zi - m_zOfRegion) < 0.5 * m_lOfRegion) ? true : false;
219
220 // std::cout <<"xi,zi= " << xi <<" "<< zi << std::endl;
221
222 if (hitRegion) {
223 returnVal = true;
224 double fl = (xi - vx) * (xi - vx) + (yi - vy) * (yi - vy) + (zi - vz) * (zi - vz);
225 fl = sqrt(fl);
226 // const double beta_mean = mom_mean.Mag() / sqrt(mom_mean.Mag2() + mass * mass);
227 // const double tofToxyPlane = fl / (c * beta_mean);
228
229 // if ((tof_up + tofToxyPlane) > tof_dn) {
230 // B2WARNING("TOF corr. stragne: " << tof_up << " " << tof_dn << " " << tofToxyPlane);
231 // continue;
232 // }
233
234 // const double deltaT4cry = fl4cry / (c * beta);
235 // double deltaT = tof_up + deltaT4cry + tofToxyPlane;
236 // const double dT = tof_up + tofToxyPlane;
237 double dT = tof_up + (tof_dn - tof_up) * fl / dpos.Mag();
238 if (tof_up > tof_dn) dT = tof_dn + (tof_up - tof_dn) * fl / dpos.Mag();
239
240 // const double dT = tof_up + (tof_dn - tof_up) * fl / dpos.Mag();
241 // std::cout <<"deltaT4cry,tof_up,tofToxyPlane,corrt= "<< deltaT4cry <<" "<< tof_up <<" "<< tofToxyPlane <<" "<< deltaT4cry - deltaT << std::endl;
242 const double pTime = m_P->getProductionTime();
243 m_P->setProductionTime(pTime - dT);
244
245 for (int iHit = 0; iHit < nHits; ++iHit) {
246 // for (const RelationElement& rel : mcp_to_hit.getElementsTo(simHits[iHit])) {
247 // if ((rel.from->getIndex() - 1) != iMCP) continue;
248 const double oldgtime = m_simHits[iHit]->getGlobalTime();
249 m_simHits[iHit]->setFlightTime(oldgtime - dT);
250 m_simHits[iHit]->setGlobalTime(oldgtime - dT);
251 // }
252 }
253 } //end of hitRegion
254 } //end of mcp loop
255
256 setReturnValue(returnVal);
257
258}
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition B2Vector3.h:435
DataType X() const
access variable X (= .at(0) without boundary check)
Definition B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition B2Vector3.h:433
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition B2Vector3.h:159
void initialize() override
Initialize variables, print info, and start CPU clock.
void event() override
Actual digitization of all hits in the CDC.
StoreArray< MCParticle > m_mcParticles
array of MCParticle
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
int getPDG() const
Return PDG code of particle.
Definition MCParticle.h:101
float getProductionTime() const
Return production time in ns.
Definition MCParticle.h:148
void setProductionTime(float time)
Set production time.
Definition MCParticle.h:373
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
Module()
Constructor.
Definition Module.cc:30
void setReturnValue(int value)
Sets the return value for this module as integer.
Definition Module.cc:220
@ 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
Low-level class to create/modify relations between StoreArrays.
Class to store a single element of a relation.
Provides access to fast ( O(log n) ) bi-directional lookups on a specified relation.
RelationIndexContainer< FROM, TO >::Element Element
Struct representing a single element in the index.
range_to getElementsTo(const TO *to) const
Return a range of all elements pointing to the given object.
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
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
bool isPrimaryParticle() const
Check if particle is a primary particle which was created by the generator (and not,...
Definition MCParticle.h:585
Abstract base class for different kinds of events.
STL namespace.