Belle II Software development
CDCCosmicSelectorModule.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/cdcCosmicSelector/CDCCosmicSelectorModule.h>
10#include <framework/logging/Logger.h>
11
12using namespace std;
13using namespace Belle2;
14
15// register module
16REG_MODULE(CDCCosmicSelector);
18{
19 // Set description
20 setDescription("Modify MCParticles for cosmics so that the global time is zero at y=0 assuming a cosmic trajectory is a line. And select cosmics passing through a trigger counter. This module works only for the event with the no. of primary charged MC particles=1. Please place this module after the event-generator and before FullSim.");
22
23 addParam("xOfCounter", m_xOfCounter, "x-position of trigger counter (cm)", -0.6);
24 addParam("yOfCounter", m_yOfCounter, "y-position of trigger counter (cm)", -13.25);
25 addParam("zOfCounter", m_zOfCounter, "z-position of trigger counter (cm)", 17.3);
26 addParam("phiOfCounter", m_phiOfCounter, "phi-angle of trigger counter (deg)", 0.);
27 addParam("wOfCounter", m_wOfCounter, "full-width of trigger counter (cm)", 7.0);
28 addParam("lOfCounter", m_lOfCounter, "full-length of trigger counter (cm)", 12.5);
29 addParam("TOF", m_tof, "Tof=1(2): TOF from production point to trigger counter (IP) is subtracted", 1);
30 addParam("TOP", m_top, "TOP from hit point to pmt in trigger counter is subtracted (assuming PMT is put at -z end of counter",
31 false);
32 addParam("propSpeed", m_propSpeed, "Light prop. speed in counter (cm/ns)", 14.4);
33 addParam("cryGenerator", m_cryGenerator, "true: CRY generator; false Cosmics generator", true);
34}
35
37{
38 m_mcParticles.isRequired();
39}
40
41
43{
44 bool returnVal = false;
45
46 int nMCPs = m_mcParticles.getEntries();
47 if (nMCPs == 0) {
48 B2ERROR("No. of MCParticles=0 in the current event.");
49 setReturnValue(returnVal);
50 return;
51 }
52
53 /*
54 const double yOfCounter = -16.25 + 3.0;
55 const double wOfCounter = 3.5;
56 const double zminOfCounter = -1.5 - 15.5;
57 const double zmaxOfCounter = -1.5 + 15.5;
58 */
59
60 // Loop over prim. charged MC particle
61 const double c = 29.9792458; //light speed (cm/ns)
62 double mass = 0.1056583715; //muon mass (GeV)
63 unsigned nPrimChgds = 0;
64
65 for (int iMCPs = 0; iMCPs < nMCPs; ++iMCPs) {
66 MCParticle* m_P = m_mcParticles[iMCPs];
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 ROOT::Math::XYZVector vertex = m_P->getProductionVertex();
94 const double vX0 = vertex.X();
95 const double vY0 = vertex.Y();
96 const double cosphi = cos(m_phiOfCounter * M_PI / 180.);
97 const double sinphi = sin(m_phiOfCounter * M_PI / 180.);
98 //(vx,vy) is prod. vertex in frame rotated in phi wrt lab. frame
99 const double vX = cosphi * vX0 + sinphi * vY0;
100 const double vY = -sinphi * vX0 + cosphi * vY0;
101 const double vZ = vertex.Z();
102
103 // std::cout <<" "<< std::endl;
104 // std::cout <<"vr,vx,vy,yz= "<< sqrt(vX*vX + vY*vY) <<" "<<vX <<" "<< vY <<" "<< vZ << std::endl;
105
106 const ROOT::Math::XYZVector momentum = m_P->getMomentum();
107 //(px,py) is momentum in frame rotated in phi wrt lab. frame
108 const double pX0 = momentum.X();
109 const double pY0 = momentum.Y();
110 const double pX = cosphi * pX0 + sinphi * pY0;
111 const double pY = -sinphi * pX0 + cosphi * pY0;
112 const double pZ = momentum.Z();
113
114 double xi = -999.;
115 double yi = -999.;
116 double zi = -999.;
117
118 //(x...,y...) is counter pos. in frame rotated in phi wrt lab. frame
119 const double xOfCounter = cosphi * m_xOfCounter + sinphi * m_yOfCounter;
120 const double yOfCounter = -sinphi * m_xOfCounter + cosphi * m_yOfCounter;
121
122 if (pY != 0.) {
123 xi = (yOfCounter - vY) * (pX / pY) + vX;
124 yi = yOfCounter;
125 zi = (yOfCounter - vY) * (pZ / pY) + vZ;
126 } else {
127 //todo improve the following 3 lines and add check accept or not also in y-direction
128 xi = 0.;
129 yi = yOfCounter;
130 zi = -vX * (pZ / pX) + vZ;
131 }
132
133 // std::cout << "xi,yi,zi= " << xi <<" "<< yi <<" "<< zi << std::endl;
134
135 bool hitCounter = (fabs(xi - xOfCounter) < 0.5 * m_wOfCounter && fabs(zi - m_zOfCounter) < 0.5 * m_lOfCounter) ? true : false;
136
137 //if hit, re-set the production time to zero at IP
138 if (hitCounter) {
139 // std::cout <<"checkxiyizi " << cosphi*xi-sinphi*yi<<" "<< sinphi*xi+cosphi *yi <<" "<< zi << std::endl;
140 returnVal = true;
141 const double p = sqrt(pX * pX + pY * pY + pZ * pZ);
142 const double cX = pX / p;
143 const double cY = pY / p;
144
145 double fl = 0.;
146 if (m_tof == 1) {
147 fl = sqrt((xi - vX) * (xi - vX) + (yi - vY) * (yi - vY) + (zi - vZ) * (zi - vZ)); // fl to counter
148 // std::cout <<"fl= " << fl << std::endl;
149 } else if (m_tof == 2) {
150 fl = -(cX * vX + cY * vY) / (cX * cX + cY * cY); //fl to origin
151 } else {
152 B2FATAL("invalid tof option !");
153 }
154 // if (fl < 0.) B2FATAL("negative tof !");
155 // std::cout << "fl,flp,dif= " << fl <<" "<< flp <<" "<< fl-flp << std::endl;
156 // double dot = pX*(xi - vX) + pY*(yi - vY);
157 // if (dot < 0.) fl *= -1.;
158
159 const double beta = p / sqrt(p * p + mass * mass);
160 double pTime = m_P->getProductionTime();
161 // const double pTimeOrg = pTime;
162
163 if (m_cryGenerator) {
164 //calculate an intersection with y=0 plane
165 double xi4cry = -999.;
166 double yi4cry = 0.;
167 double zi4cry = -999.;
168 if (pY0 != 0.) {
169 xi4cry = (0. - vY0) * (pX0 / pY0) + vX0;
170 zi4cry = (0. - vY0) * (pZ / pY0) + vZ;
171 } else {
172 xi4cry = 0.;
173 zi4cry = -vX0 * (pZ / pX0) + vZ;
174 }
175
176 const double fl4cry = sqrt((xi4cry - vX0) * (xi4cry - vX0) + (yi4cry - vY0) * (yi4cry - vY0) + (zi4cry - vZ) *
177 (zi4cry - vZ)); // fl to y=0 plane
178
179 //reset production time (to the time which old CRY set) once
180 pTime += fl4cry / (c * beta);
181 }
182
183 const double tofToCounter = fl / (c * beta);
184 const double topToPMT = m_top ? hypot(xi - xOfCounter, zi - (m_zOfCounter - 0.5 * m_lOfCounter)) / m_propSpeed : 0.;
185 /*
186 std::cout << "xi = " << xi << std::endl;
187 std::cout << "xOfCounter = " << xOfCounter << std::endl;
188 std::cout << "zi = " << zi << std::endl;
189 std::cout << "m_zOfCounter= " << m_zOfCounter << std::endl;
190 std::cout << "m_lOfCounter= " << m_lOfCounter << std::endl;
191 std::cout << "m_propSpeed = " << m_propSpeed << std::endl;
192 std::cout << "topToPMT = " << topToPMT << std::endl;
193 */
194 m_P->setProductionTime(pTime - tofToCounter - topToPMT);
195 // std::cout <<"org,mod= " << pTimeOrg << m_P->getProductionTime() << std::endl;
196 //if not hit, reverse the momentum vector so that the particle will not be simulated
197 //better to use condition in basf2...
198 } else {
199 // m_P->setMomentum(-1. * momentum);
200
201 }
202 } // end loop over MCParticles
203 setReturnValue(returnVal);
204}
double m_propSpeed
Light speed in counter (cm/ns)
void initialize() override
Initialize variables, print info, and start CPU clock.
double m_wOfCounter
full-width of counter (cm)
void event() override
Actual digitization of all hits in the CDC.
double m_phiOfCounter
phi-angle of counter (deg)
double m_lOfCounter
full-length of counter (cm)
bool m_cryGenerator
cry or cosmics generator
StoreArray< MCParticle > m_mcParticles
array of MCParticle
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
ROOT::Math::XYZVector getProductionVertex() const
Return production vertex position.
Definition MCParticle.h:178
int getPDG() const
Return PDG code of particle.
Definition MCParticle.h:101
float getProductionTime() const
Return production time in ns.
Definition MCParticle.h:148
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition MCParticle.h:187
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
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
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.