11 #include <beast/microtpc/modules/TpcDigitizerModule.h>
12 #include <beast/microtpc/dataobjects/MicrotpcSimHit.h>
14 #include <mdst/dataobjects/MCParticle.h>
15 #include <framework/logging/Logger.h>
16 #include <framework/gearbox/GearDir.h>
17 #include <framework/core/RandomNumbers.h>
21 #include <boost/foreach.hpp>
25 using namespace microtpc;
40 setDescription(
"Microtpc digitizer module");
43 addParam(
"LowerTimingCut", m_lowerTimingCut,
"Lower timing cut", 0.);
44 addParam(
"UpperTimingCut", m_upperTimingCut,
"Upper timing cut", 1000000.);
47 TpcDigitizerModule::~TpcDigitizerModule()
53 B2INFO(
"TpcDigitizer: Initializing");
54 m_microtpcHit.registerInDataStore();
60 fctToT_Calib1 =
new TF1(
"fctToT_Calib1",
"[0]*(x/[3]+[1])/(x/[3]+[2])", 0., 100000.);
63 fctToT_Calib2 =
new TF1(
"fctToT_Calib2",
"[0]*(x/[3]+[1])/(x/[3]+[2])", 0., 100000.);
94 std::vector<double> T0(
m_nTPC,
98 for (
const auto& microtpcSimHit : microtpcSimHits) {
99 const int detNb = microtpcSimHit.getdetNb();
100 const TVector3 simHitPosition = microtpcSimHit.gettkPos();
101 TVector3 chipPosition(
102 simHitPosition.X() / 100. -
m_TPCCenter[detNb].X(),
103 simHitPosition.Y() / 100. -
m_TPCCenter[detNb].Y(),
106 chipPosition.RotateZ(-
m_TPCAngleZ[detNb] * TMath::DegToRad());
107 chipPosition.RotateX(-
m_TPCAngleX[detNb] * TMath::DegToRad());
108 const double T = chipPosition.Z() +
m_z_DG / 2.;
114 for (
auto& val : T0) {
124 for (
const auto& microtpcSimHit : microtpcSimHits) {
126 const int PDGid = microtpcSimHit.gettkPDG();
128 if (PDGid != 1000020040 || PDGid != 1000060120 || PDGid != 1000080160 || PDGid != 2212) {
133 const int detNb = microtpcSimHit.getdetNb();
134 const double edep = microtpcSimHit.getEnergyDep();
135 const double niel = microtpcSimHit.getEnergyNiel();
136 const int pdg = microtpcSimHit.gettkPDG();
137 const int trkID = microtpcSimHit.gettkID();
139 const TVector3 simHitPosition = microtpcSimHit.gettkPos();
140 TVector3 chipPosition(
141 simHitPosition.X() / 100. -
m_TPCCenter[detNb].X(),
142 simHitPosition.Y() / 100. -
m_TPCCenter[detNb].Y(),
145 chipPosition.RotateZ(-
m_TPCAngleZ[detNb] * TMath::DegToRad());
146 chipPosition.RotateX(-
m_TPCAngleX[detNb] * TMath::DegToRad());
148 TVector3 ChipPosition(0, 0, 0);
150 if (detNb == 0 || detNb == 3) {
151 ChipPosition.SetX(-chipPosition.Y());
152 ChipPosition.SetY(-chipPosition.X());
153 ChipPosition.SetZ(chipPosition.Z() +
m_z_DG / 2.);
155 if (detNb == 1 || detNb == 2) {
156 ChipPosition.SetX(chipPosition.Y());
157 ChipPosition.SetY(chipPosition.X());
158 ChipPosition.SetZ(chipPosition.Z() +
m_z_DG / 2.);
162 ChipPosition.SetX(chipPosition.Y());
163 ChipPosition.SetY(chipPosition.X());
164 ChipPosition.SetZ(chipPosition.Z() +
m_z_DG / 2);
183 (0. < ChipPosition.Z() && ChipPosition.Z() <
m_z_DG) &&
188 const double ionEn = (edep - niel) * 1e3;
198 const double meanEl = ionEn * 1e3 /
m_Workfct;
199 const double sigma = sqrt(
m_Fanofac * meanEl);
200 const int NbEle = (int)gRandom->Gaus(meanEl, sigma);
201 const double NbEle_real = NbEle - NbEle *
m_GasAbs * chipPosition.Z();
204 for (
int ie = 0; ie < (int)NbEle_real; ie++) {
208 double x_DG, y_DG, z_DG, t_DG;
209 Drift(ChipPosition.X(),
222 for (
int ig1 = 0; ig1 < (int)GEM_gain1; ig1++) {
225 const TVector2 GEM1(
GEMGeo1(x_DG, y_DG));
228 double x_TG, y_TG, z_TG, t_TG;
233 for (
int ig2 = 0; ig2 < (int)GEM_gain2; ig2++) {
236 const TVector2 GEM2(
GEMGeo2(x_TG, y_TG));
239 double x_CG, y_CG, z_CG, t_CG;
246 int quT = gRandom->Uniform(-1, 1);
247 int bci = (int)((t_DG + t_TG + t_CG - T0[detNb]) / (double)
m_PixelTimeBin) + quT;
254 (0 <= bci && bci < MAXtSIZE)) {
325 double sl,
double vd)
330 x2 = x1 + gRandom->Gaus(0., sqrt(z1) * st);
332 y2 = y1 + gRandom->Gaus(0., sqrt(z1) * st);
334 z2 = z1 + gRandom->Gaus(0., sqrt(z1) * sl);
338 x2 = -1000; y2 = -1000; z2 = -1000; t2 = -1000;
343 static const double sqrt3o4 = std::sqrt(3. / 4.);
346 int yint = (int)(y1 / (sqrt3o4 *
m_GEMpitch) + 0.5);
353 return TVector2(x2, y2);
358 static const double sqrt3o4 = std::sqrt(3. / 4.);
361 int yint = (int)(x1 / (sqrt3o4 *
m_GEMpitch) + 0.5);
368 return TVector2(x2, y2);
379 std::vector<int> col;
380 std::vector<int> row;
381 std::vector<int> ToT;
382 std::vector<int> bci;
391 const auto& key = keyValuePair.first;
393 int i = std::get<0>(key);
395 int j = std::get<1>(key);
397 if (
m_dchip_map[std::tuple<int, int>(i, j)] == 1) {
400 const int quE = gRandom->Uniform(0, 2);
404 for (
auto& keyValuePair2 :
m_dchip) {
405 const auto& key2 = keyValuePair2.first;
406 int k = std::get<2>(key2);
407 if (
m_dchip[std::tuple<int, int, int>(i, j, k)] > thresEl) {
418 for (
auto& keyValuePair2 :
m_dchip) {
419 const auto& key2 = keyValuePair2.first;
420 int k = std::get<2>(key2);
423 NbOfEl +=
m_dchip[std::tuple<int, int, int>(i, j, k)];
427 if (NbOfEl > thresEl && NbOfEl <= 45.*
m_TOTQ1) {
431 }
else if (NbOfEl > 800.*
m_TOTQ1) {
453 if (NbOfEl > thresEl && NbOfEl <= 45.*
m_TOTQ1) {
457 }
else if (NbOfEl > 800.*
m_TOTQ1) {
479 if (bci.size() > 0) {
483 sort(t0.begin(), t0.end());
486 for (
int j = 0; j < (int)bci.size(); j++) {
689 GearDir content =
GearDir(
"/Detector/DetectorComponent[@name=\"MICROTPC\"]/Content/");
692 BOOST_FOREACH(
const GearDir & activeParams, content.getNodes(
"Active")) {
701 m_phase = content.getDouble(
"Phase");
718 m_TOTA1 = content.getDouble(
"TOTA1");
719 m_TOTB1 = content.getDouble(
"TOTB1");
720 m_TOTC1 = content.getDouble(
"TOTC1");
721 m_TOTQ1 = content.getDouble(
"TOTQ1");
722 m_TOTA2 = content.getDouble(
"TOTA2");
723 m_TOTB2 = content.getDouble(
"TOTB2");
724 m_TOTC2 = content.getDouble(
"TOTC2");
725 m_TOTQ2 = content.getDouble(
"TOTQ2");
726 m_z_DG = content.getDouble(
"z_DG");
727 m_z_TG = content.getDouble(
"z_TG");
728 m_z_CG = content.getDouble(
"z_CG");
729 m_Dl_DG = content.getDouble(
"Dl_DG");
730 m_Dl_TG = content.getDouble(
"Dl_TG");
731 m_Dl_CG = content.getDouble(
"Dl_CG");
732 m_Dt_DG = content.getDouble(
"Dt_DG");
733 m_Dt_TG = content.getDouble(
"Dt_TG");
734 m_Dt_CG = content.getDouble(
"Dt_CG");
735 m_v_DG = content.getDouble(
"v_DG");
736 m_v_TG = content.getDouble(
"v_TG");
737 m_v_CG = content.getDouble(
"v_CG");
738 m_Workfct = content.getDouble(
"Workfct");
739 m_Fanofac = content.getDouble(
"Fanofac");
740 m_GasAbs = content.getDouble(
"GasAbs");
742 B2INFO(
"TpcDigitizer: Aquired tpc locations and gas parameters");
743 B2INFO(
" from MICROTPC.xml. There are " <<
m_nTPC <<
" TPCs implemented");