36 #include "G4SDManager.hh"
39 #include "G4VProcess.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
62 preStepPoint(nullptr),
63 postStepPoint(nullptr) {
68 SetVerboseLevel(verbn);
72 if (
name ==
"TotemHitsT1") {
74 }
else if (
name ==
"TotemHitsT2Si") {
76 }
else if (
name ==
"TotemHitsT2Gem") {
78 }
else if (
name ==
"TotemHitsRP") {
81 edm::LogWarning(
"ForwardSim") <<
"TotemSD: ReadoutName not supported\n";
107 LogDebug(
"ForwardSim") <<
"TotemSD : Initialize called for " << GetName();
120 int thehc_entries =
theHC->entries();
121 for (
int j = 0;
j < thehc_entries &&
j < 15000;
j++) {
124 LogDebug(
"ForwardSim") <<
"HIT NUMERO " <<
j <<
"unit ID = " << aHit->
getUnitID() <<
"\n"
126 <<
"enrty z " << aHit->
getEntry().z() <<
"\n"
146 LogDebug(
"ForwardSim") <<
"TotemSD: Collection " <<
theHC->GetName();
147 theHC->PrintAllHits();
157 LogDebug(
"ForwardSim") <<
" Dispatched BeginOfEvent for " << GetName() <<
" !";
164 G4ThreeVector localPoint;
165 const G4VTouchable* touch =
preStepPoint->GetTouchable();
166 localPoint = touch->GetHistory()->GetTopTransform().TransformPoint(global);
181 edeposit = aStep->GetTotalEnergyDeposit();
192 Pabs = aStep->GetPreStepPoint()->GetMomentum().mag() /
GeV;
193 Tof = aStep->GetPostStepPoint()->GetGlobalTime() / nanosecond;
195 Eloss = aStep->GetTotalEnergyDeposit() /
GeV;
198 ThetaAtEntry = aStep->GetPreStepPoint()->GetPosition().theta() / deg;
199 PhiAtEntry = aStep->GetPreStepPoint()->GetPosition().phi() / deg;
209 edm::LogWarning(
"ForwardSim") <<
"***** TotemSD error: primaryID = " <<
primaryID <<
" maybe detector name changed";
227 int thehc_entries =
theHC->entries();
228 for (
int j = 0;
j < thehc_entries && !
found;
j++) {
248 LogDebug(
"ForwardSim") <<
"TotemSD CreateNewHit for"
250 <<
" MVid = " <<
currentPV->GetMother()->GetCopyNo() <<
" Unit " <<
unitID <<
"\n"
252 <<
theTrack->GetTrackID() <<
" which is a " <<
theTrack->GetDefinition()->GetParticleName();
260 if (
theTrack->GetCreatorProcess() !=
nullptr)
263 LogDebug(
"ForwardSim") <<
"NO process";
314 G4ThreeVector _PosizioEvo;
329 const G4ThreeVector& Pos,
double vx,
double vy,
double vz,
double pabs,
int& accettanza) {
332 G4ThreeVector PosEvo;
333 double ThetaX = std::atan((Pos.x() - vx) / (Pos.z() - vz));
334 double ThetaY = std::atan((Pos.y() - vy) / (Pos.z() - vz));
335 double X_at_0 = (vx - ((Pos.x() - vx) / (Pos.z() - vz)) * vz) / 1000.;
336 double Y_at_0 = (vy - ((Pos.y() - vy) / (Pos.z() - vz)) * vz) / 1000.;
344 double csi =
std::abs((7000. - pabs) / 7000.);
348 double x_par[no_rp + 1];
349 double y_par[no_rp + 1];
351 double rp[no_rp] = {141., 149., 198., 220.};
353 double leffx[][2] = {{122.5429, -46.9312}, {125.4194, -49.1849}, {152.6, -81.157}, {98.8914, -131.8390}};
355 double leffy[][2] = {{124.2314, -55.4852}, {127.7825, -57.4503}, {179.455, -76.274}, {273.0931, -40.4626}};
357 double avx[][2] = {{0.515483, -1.0123}, {0.494122, -1.0534}, {0.2217, -1.483}, {0.004633, -1.0719}};
359 double avy[][2] = {{0.371418, -1.6327}, {0.349035, -1.6955}, {0.0815, -2.59}, {0.007592, -4.0841}};
361 double ddx[][4] = {{-0.082336, -0.092513, 112.3436, -82.5029},
362 {-0.086927, -0.097670, 114.9513, -82.9835},
363 {-0.092117, -0.0915, 180.6236, -82.443},
364 {-0.050470, 0.058837, 208.1106, 20.8198}};
366 double detlim[][2] = {{0, 0}, {0.0028, 0.0021}, {0, 0}, {0.0008, 0.0013}};
368 double pipelim[][2] = {{0.026, 0.026}, {0.04, 0.04}, {0.0226, 0.0177}, {0.04, 0.04}};
370 for (
int j = 0;
j < no_rp;
j++) {
373 y_par[
j] = ThetaY * (leffy[
j][0] + leffy[
j][1] * csi) + (avy[
j][0] + avy[
j][1] * csi) * Y_at_0;
374 x_par[
j] = ThetaX * (leffx[
j][0] + leffx[
j][1] * csi) + (avx[
j][0] + avx[
j][1] * csi) * X_at_0 -
375 csi * (ddx[
j][0] + (ddx[
j][1] + ddx[
j][2] * ThetaX) * csi + ddx[
j][3] * ThetaX);
379 if (
std::abs(y_par[0]) < pipelim[0][1] &&
sqrt((y_par[0] * y_par[0]) + (x_par[0] * x_par[0])) < pipelim[0][0]) {
381 if ((
sqrt((y_par[1] * y_par[1]) + (x_par[1] * x_par[1])) < pipelim[1][0]) &&
382 (
std::abs(y_par[1]) > detlim[1][1] || x_par[1] > detlim[1][0])) {
388 if (
std::abs(y_par[0]) < pipelim[0][1] &&
sqrt((y_par[0]) * (y_par[0]) + (x_par[0]) * (x_par[0])) < pipelim[0][0]) {
390 if (
std::abs(y_par[2]) < pipelim[2][1] &&
sqrt((y_par[2] * y_par[2]) + (x_par[2] * x_par[2])) < pipelim[2][0]) {
392 if ((
sqrt((y_par[3] * y_par[3]) + (x_par[3] * x_par[3])) < pipelim[3][0]) &&
393 (
std::abs(y_par[3]) > detlim[3][1] || x_par[3] > detlim[3][0])) {
396 PosEvo.setX(1000 * x_par[3]);
397 PosEvo.setY(1000 * y_par[3]);
398 PosEvo.setZ(1000 * rp[3]);
400 PosEvo.setZ(-1000 * rp[3]);
442 if (
hit ==
nullptr) {
443 edm::LogWarning(
"ForwardSim") <<
"TotemSD: hit to be stored is NULL !!";