9#include <framework/gearbox/Gearbox.h>
10#include <framework/gearbox/GearDir.h>
11#include <framework/logging/Logger.h>
12#include <framework/utilities/FileSystem.h>
14#include <cdc/geometry/CDCGeometryPar.h>
15#include <cdc/geometry/CDCGeoControlPar.h>
16#include <cdc/simulation/CDCSimControlPar.h>
17#include <cdc/utilities/OpenFile.h>
22#include <boost/iostreams/filtering_stream.hpp>
24#include <Math/ChebyshevPol.h>
45 if ((*m_t0FromDB).isValid()) {
52 if ((*m_badWireFromDB).isValid()) {
58 if ((*m_badBoardsFromDB).isValid()) {
64 if ((*m_propSpeedFromDB).isValid()) {
71 if ((*m_timeWalkFromDB).isValid()) {
78 if ((*m_xtRelFromDB).isValid()) {
85 if ((*m_sResolFromDB).isValid()) {
92 if ((*m_fFactorFromDB).isValid()) {
99 if ((*m_chMapFromDB).isValid()) {
106 if ((*m_displacementFromDB).isValid()) {
113 if ((*m_alignmentFromDB).isValid()) {
121 if ((*m_misalignmentFromDB).isValid()) {
130 if ((*m_eDepToADCConversionsFromDB).isValid()) {
143 B2WARNING(
"CDCGeometryPar: Strange that readFromDB is not called! Please make sure that CDC is included in Geometry.");
182 for (
unsigned i = 0; i < 4; ++i) {
184 for (
unsigned j = 0; j < 2; ++j)
187 for (
unsigned i = 0; i < c_maxNSenseLayers; ++i) {
199 for (
unsigned i = 0; i < c_maxNFieldLayers; ++i) {
205 for (
unsigned L = 0; L < c_maxNSenseLayers; ++L) {
206 for (
unsigned C = 0; C < c_maxNDriftCells; ++C) {
207 for (
unsigned i = 0; i < 3; ++i) {
215 for (
unsigned i = 0; i < 7; ++i) {
225 for (
unsigned L = 0; L < c_maxNSenseLayers; ++L) {
226 for (
unsigned i = 0; i < 2; ++i) {
227 for (
unsigned alpha = 0; alpha < c_maxNAlphaPoints; ++alpha) {
228 for (
unsigned theta = 0; theta < c_maxNThetaPoints; ++theta) {
229 for (
unsigned xtparam = 0; xtparam < c_nXTParams; ++xtparam) {
230 m_XT[L][i][alpha][theta][xtparam] = 0.;
233 for (
unsigned sigmaparam = 0; sigmaparam < c_nSigmaParams; ++sigmaparam) {
234 m_Sigma[L][i][alpha][theta][sigmaparam] = 0.;
241 for (
unsigned board = 0; board < c_nBoards; ++board) {
242 for (
unsigned i = 0; i < 2; ++i) {
245 for (
unsigned channel = 0; channel < 48; ++channel) {
250 for (
unsigned superLayer = 0; superLayer < c_nSuperLayers; ++superLayer) {
251 for (
unsigned layer = 0; layer < 8; ++layer) {
264 m_rWall[0] = geom.getInnerWall(2).getRmin();
265 m_zWall[0][0] = geom.getInnerWall(0).getZbwd();
266 m_zWall[0][1] = geom.getInnerWall(0).getZfwd();
268 m_rWall[1] = geom.getInnerWall(0).getRmax();
269 m_zWall[1][0] = geom.getInnerWall(0).getZbwd();
270 m_zWall[1][1] = geom.getInnerWall(0).getZbwd();
273 m_rWall[2] = geom.getOuterWall(0).getRmin();
274 m_zWall[2][0] = geom.getOuterWall(0).getZbwd();
275 m_zWall[2][1] = geom.getOuterWall(0).getZfwd();
277 m_rWall[3] = geom.getOuterWall(1).getRmax();
278 m_zWall[3][0] = geom.getOuterWall(0).getZbwd();
279 m_zWall[3][1] = geom.getOuterWall(0).getZfwd();
288 B2DEBUG(100,
"CDCGeometryPar: Define a mixture of gases and wires in the tracking volume.");
291 B2FATAL(
"CDCGeometryPar: Materialdefinition=2 is disabled for now.");
293 B2FATAL(
"CDCGeometryPar: Materialdefinition mode you specify is invalid.");
304 GearDir content =
GearDir(
"/Detector/DetectorComponent[@name=\"CDC\"]/Content/");
312 for (
const auto& sense : geom.getSenseLayers()) {
313 uint layerId = sense.getId();
318 m_nWires[layerId] = sense.getNWires();
320 m_offSet[layerId] = sense.getOffset();
333 B2FATAL(
"CDCGeometryPar: invalid wire z definition mode specified");
337 const int nWires =
m_nWires[layerId];
338 for (
int iCell = 0; iCell < nWires; ++iCell) {
345 for (
const auto& field : geom.getFieldLayers()) {
346 uint layerId = field.getId();
378 B2FATAL(
"HardwareClockSettings payloads are not valid.");
379 const double officialClockFreq4TDC = 2 *
m_clockSettings->getAcceleratorRF();
381 B2WARNING(
"ClockFreq4TDC changed from cdclocal " << scientific << setprecision(6) <<
m_clockFreq4TDC <<
" to official " <<
382 officialClockFreq4TDC <<
" (GHz) (difference larger than 0.01%)");
385 B2DEBUG(100,
"CDCGeometryPar: Clock freq. for TDC= " <<
m_clockFreq4TDC <<
" (GHz).");
387 B2DEBUG(100,
"CDCGeometryPar: TDC bin width= " <<
m_tdcBinWidth <<
" (ns).");
398 B2DEBUG(100,
"CDCGeometryPar: Load displacement params. (=1); not load (=0):" <<
m_displacement);
401 B2DEBUG(100,
"CDCGeometryPar: Read displacement from DB");
410 B2DEBUG(100,
"CDCGeometryPar: Load alignment params. (=1); not load (=0):" <<
414 B2DEBUG(100,
"CDCGeometryPar: Read alignment from DB");
423 B2DEBUG(100,
"CDCGeometryPar: Load misalignment params. (=1); not load (=0):" <<
427 B2DEBUG(100,
"CDCGeometryPar: Read misalignment from DB");
440 B2FATAL(
"ModifiedLeftRightFlag = true is disabled for now; need to update a G4-related code in framework...");
449 B2DEBUG(100,
"CDCGeometryPar: Read xt from DB");
456 B2DEBUG(100,
"CDCGeometryPar: Read sigma from DB");
463 B2DEBUG(100,
"CDCGeometryPar: Read fudge factors from DB");
470 B2DEBUG(100,
"CDCGeometryPar: Read prop-speed from DB");
477 B2DEBUG(100,
"CDCGeometryPar: Read t0 from DB");
484 B2DEBUG(100,
"CDCGeometryPar: Read badwire from DB");
488 B2FATAL(
"Text file input mode for bdwires is disabled now!");
492 B2DEBUG(100,
"CDCGeometryPar: Read ch-map from DB");
499 B2DEBUG(100,
"CDCGeometryPar: Read time-walk from DB");
504 B2DEBUG(100,
"CDCGeometryPar: Time-walk param. mode= " <<
m_twParamMode);
507 B2DEBUG(29,
"CDCGeometryPar: Read EDepToADC from DB");
508 if ((*m_eDepToADCConversionsFromDB).isValid()) {
580 std::string fileName0;
585 }
else if (set == c_Misaligned) {
587 }
else if (set == c_Aligned) {
593 }
else if (set == c_Misaligned) {
595 }
else if (set == c_Aligned) {
602 boost::iostreams::filtering_istream ifs;
607 double back[np], fwrd[np], tension;
612 for (
int i = 0; i < np; ++i) {
615 for (
int i = 0; i < np; ++i) {
621 if (ifs.eof())
break;
629 for (
int i = 0; i < np; ++i) {
633 }
else if (set == c_Misaligned) {
636 }
else if (set == c_Aligned) {
648 }
else if (set == c_Misaligned) {
653 }
else if (set == c_Aligned) {
672 if (nRead !=
m_nSenseWires) B2FATAL(
"CDCGeometryPar::readWirePositionParams: #lines read-in (=" << nRead <<
673 ") is inconsistent with total #sense wires (=" <<
m_nSenseWires <<
") !");
676 boost::iostreams::close(ifs);
684 for (
unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
691 auto layerID =
WireID(iL, 511);
702 for (
unsigned iC = 0; iC <
m_nWires[iL]; ++iC) {
714 m_BWirPosAlign[iL][iC][0] = d_layerXbwd + cos(d_layerPhiBwd) * wireXbwd + sin(d_layerPhiBwd) * wireYbwd;
715 m_BWirPosAlign[iL][iC][1] = d_layerYbwd - sin(d_layerPhiBwd) * wireXbwd + cos(d_layerPhiBwd) * wireYbwd;
718 m_FWirPosAlign[iL][iC][0] = d_layerXfwd + cos(d_layerPhiFwd) * wireXfwd + sin(d_layerPhiFwd) * wireYfwd;
719 m_FWirPosAlign[iL][iC][1] = d_layerYfwd - sin(d_layerPhiFwd) * wireXfwd + cos(d_layerPhiFwd) * wireYfwd;
725 double back[np], fwrd[np];
727 for (
unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
733 for (
unsigned iC = 0; iC <
m_nWires[iL]; ++iC) {
744 for (
int i = 0; i < np; ++i) {
768 double back[np], fwrd[np];
770 for (
unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
776 for (
unsigned iC = 0; iC <
m_nWires[iL]; ++iC) {
787 for (
int i = 0; i < np; ++i) {
822 fileName0 = gbxParams.
getString(
"xt4ReconFileName");
825 boost::iostreams::filtering_istream ifs;
849 unsigned short nAlphaBins = 0;
850 if (ifs >> nAlphaBins) {
851 if (nAlphaBins == 0 || nAlphaBins > c_maxNAlphaPoints) B2FATAL(
"Fail to read alpha bins !");
853 B2FATAL(
"Fail to read alpha bins !");
856 double alpha0, alpha1, alpha2;
857 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
858 ifs >> alpha0 >> alpha1 >> alpha2;
863 unsigned short nThetaBins = 0;
864 if (ifs >> nThetaBins) {
865 if (nThetaBins == 0 || nThetaBins > c_maxNThetaPoints) B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
867 B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
870 double theta0, theta1, theta2;
872 for (
unsigned short i = 0; i < nThetaBins; ++i) {
873 ifs >> theta0 >> theta1 >> theta2;
878 unsigned short iCL, iLR;
879 const unsigned short npx = c_nXTParams - 1;
881 double theta, alpha, dummy1;
886 if (np <= 0 || np > npx) B2FATAL(
"CDCGeometryPar: no. of xt-params. outside limits !");
888 const double epsi = 0.1;
896 ifs >> theta >> alpha >> dummy1 >> iLR;
897 for (
int i = 0; i < np; ++i) {
902 for (
unsigned short i = 0; i < nThetaBins; ++i) {
908 if (itheta < 0) B2FATAL(
"CDCGeometryPar: thetas in xt.dat are inconsistent !");
911 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
917 if (ialpha < 0) B2FATAL(
"CDCGeometryPar: alphas in xt.dat are inconsistent !");
919 for (
int i = 0; i < np; ++i) {
920 m_XT[iCL][iLR][ialpha][itheta][i] = xtc[i];
923 double boundT = xtc[6];
925 m_XT[iCL][iLR][ialpha][itheta][np] = ROOT::Math::Chebyshev5(boundT, xtc[0], xtc[1], xtc[2], xtc[3], xtc[4], xtc[5]);
927 m_XT[iCL][iLR][ialpha][itheta][np] =
938 boost::iostreams::close(ifs);
941 const double degrad = M_PI / 180.;
942 for (
unsigned i = 0; i < nAlphaBins; ++i) {
945 for (
unsigned i = 0; i < nThetaBins; ++i) {
968 fileName0 = gbxParams.
getString(
"sigma4ReconFileName");
976 unsigned short nAlphaBins = 0;
977 if (ifs >> nAlphaBins) {
978 if (nAlphaBins == 0 || nAlphaBins > c_maxNAlphaPoints) B2FATAL(
"Fail to read alpha bins !");
980 B2FATAL(
"Fail to read alpha bins !");
984 double alpha0, alpha1, alpha2;
985 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
986 ifs >> alpha0 >> alpha1 >> alpha2;
992 unsigned short nThetaBins = 0;
993 if (ifs >> nThetaBins) {
994 if (nThetaBins == 0 || nThetaBins > c_maxNThetaPoints) B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
996 B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
1000 double theta0, theta1, theta2;
1002 for (
unsigned short i = 0; i < nThetaBins; ++i) {
1003 ifs >> theta0 >> theta1 >> theta2;
1008 unsigned short np = 0;
1009 unsigned short iCL, iLR;
1010 double sigma[c_nSigmaParams];
1011 double theta, alpha;
1017 if (np > c_nSigmaParams) B2FATAL(
"CDCGeometryPar: no. of sigma-params. outside limits !");
1021 const double epsi = 0.1;
1023 while (ifs >> iCL) {
1029 ifs >> theta >> alpha >> iLR;
1031 for (
int i = 0; i < np; ++i) {
1036 for (
unsigned short i = 0; i < nThetaBins; ++i) {
1042 if (itheta < 0) B2FATAL(
"CDCGeometryPar: thetas in sigma.dat are inconsistent !");
1045 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
1051 if (ialpha < 0) B2FATAL(
"CDCGeometryPar: alphas in sigma.dat are inconsistent !");
1053 for (
int i = 0; i < np; ++i) {
1054 m_Sigma[iCL][iLR][ialpha][itheta][i] = sigma[i];
1061 const double degrad = M_PI / 180.;
1062 for (
unsigned i = 0; i < nAlphaBins; ++i) {
1065 for (
unsigned i = 0; i < nThetaBins; ++i) {
1078 fileName0 = gbxParams.
getString(
"fudgeFactorFileName");
1080 B2WARNING(
"readFFactor is not ready! " << fileName0);
1090 fileName0 = gbxParams.
getString(
"propSpeed4ReconFileName");
1103 if (ifs.eof())
break;
1109 if (
m_debug) B2DEBUG(150, iL <<
" " << speed);
1112 if (nRead != c_maxNSenseLayers) B2FATAL(
"CDCGeometryPar::readPropSpeed: #lines read-in (=" << nRead <<
1113 ") is inconsistent with total #layers (=" << c_maxNSenseLayers <<
") !");
1158 fileName0 = gbxParams.
getString(
"t04ReconFileName");
1170 ifs >> iL >> iC >> t0;
1176 if (ifs.eof())
break;
1183 B2DEBUG(150, iL <<
" " << iC <<
" " << t0);
1187 if (nRead !=
m_nSenseWires) B2FATAL(
"CDCGeometryPar::readT0: #lines read-in (=" << nRead <<
1188 ") is inconsistent with total #sense wires (=" <<
m_nSenseWires <<
") !");
1239 fileName0 = gbxParams.
getString(
"tw4ReconFileName");
1246 unsigned short nPars(0);
1249 B2FATAL(
"CDCGeometryPar::readTW: invalid mode specified!");
1252 B2FATAL(
"CDCGeometryPar::readTW: invalid #params specified!");
1255 unsigned iBoard = 0;
1258 while (ifs >> iBoard) {
1259 for (
unsigned short i = 0; i < nPars; ++i) {
1265 if (nRead != c_nBoards) B2FATAL(
"CDCGeometryPar::readTW: #lines read-in (=" << nRead <<
") is inconsistent with #boards (=" <<
1283 unsigned short iSL, iL, iW, iB, iC;
1288 ifs >> iSL >> iL >> iW >> iB >> iC;
1289 if (ifs.eof())
break;
1294 m_wireToBoard.insert(pair<WireID, unsigned short>(wID, iB));
1297 if (nRead !=
m_nSenseWires) B2FATAL(
"CDCGeometryPar::readChMap: #lines read-in (=" << nRead <<
1298 ") is inconsistent with #sense-wires (="
1311 fileName0 = gbxParams.
getString(
"fudgeFactorFileName");
1316 std::string fileName1 =
"/data/cdc/" + fileName0;
1319 if (fileName ==
"") {
1323 if (fileName ==
"") {
1324 B2FATAL(
"CDC::openFile: " << fileName0 <<
" not exist!");
1327 B2DEBUG(29,
"CDC::openFile: open " << fileName);
1328 ifs.open(fileName.c_str());
1329 if (!ifs) B2FATAL(
"CDC::openFile: cannot open " << fileName <<
" !");
1332 unsigned short paramMode(0), nParams(0);
1333 ifs >> paramMode >> nParams;
1334 if (paramMode > 1) B2FATAL(
"Param mode > 1!");
1335 if (nParams > 7) B2FATAL(
"No. of params. > 7!");
1336 unsigned short groupId(0);
1338 B2DEBUG(29, paramMode <<
" " << nParams <<
" " << groupId);
1339 if (groupId > 0) B2FATAL(
"GgroupId > 0!");
1342 unsigned short cLMin[c_nSuperLayers], cLMax[c_nSuperLayers];
1345 for (
unsigned int sl = 1; sl < c_nSuperLayers; ++sl) {
1346 cLMin[sl] = cLMax[0] + 6 * sl - 5;
1347 cLMax[sl] = cLMax[0] + 6 * sl;
1350 unsigned short id = 0;
1352 unsigned short nRead = 0;
1354 for (
unsigned short i = 0; i < nParams; ++i) {
1356 for (
unsigned short cL = cLMin[
id]; cL <= cLMax[id]; ++cL) {
1357 for (
unsigned short cell = 0; cell <
m_nWires[cL]; ++cell) {
1364 if (nRead > c_nSuperLayers) B2FATAL(
"No. of read in lines > " << c_nSuperLayers <<
" !");
1374 for (
unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1375 for (
unsigned short iW = 0; iW < c_maxNDriftCells; ++iW) {
1380 for (
auto const& ent : (*m_t0FromDB)->getT0s()) {
1383 const unsigned short iW = wid.
getIWire();
1394 double oldMeanT0 = 0;
1395 unsigned short it1 = 0;
1396 for (
unsigned short it = 0; it < maxIt; ++it) {
1398 double effiSum = 0.;
1401 for (
unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1402 for (
unsigned short iW = 0; iW <
m_nWires[iCL]; ++iW) {
1403 if (
m_t0[iCL][iW] < minT0 ||
m_t0[iCL][iW] > maxT0)
continue;
1418 B2DEBUG(29, it <<
" " << effiSum <<
" " <<
m_meanT0 <<
" " << stdvT0);
1419 if (fabs(
m_meanT0 - oldMeanT0) < epsi)
break;
1424 B2FATAL(
"Wire efficiency sum <= 0!");
1427 if (it1 == maxIt - 1) B2WARNING(
"Max. iterations(=" << maxIt <<
") needed to calculate the mean t0. Strange.");
1447 for (
unsigned short iCL = 0; iCL < (*m_propSpeedFromDB)->getEntries(); ++iCL) {
1459 for (
unsigned short iBd = 0; iBd < (*m_timeWalkFromDB)->getEntries(); ++iBd) {
1460 int np = ((*m_timeWalkFromDB)->getTimeWalkParams(iBd)).size();
1461 for (
int i = 0; i < np; ++i) {
1462 m_timeWalkCoef[iBd][i] = ((*m_timeWalkFromDB)->getTimeWalkParams(iBd))[i];
1488 for (
unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1495 for (
unsigned short iLR = 0; iLR < 2; ++iLR) {
1498 const std::vector<float> params = (*m_xtRelFromDB)->getXtParams(iCL, iLR, iA, iT);
1499 unsigned short np = params.size();
1501 for (
unsigned short i = 0; i < np; ++i) {
1502 m_XT[iCL][iLR][iA][iT][i] = params[i];
1505 double boundT =
m_XT[iCL][iLR][iA][iT][6];
1507 m_XT[iCL][iLR][iA][iT][np] = ROOT::Math::Chebyshev5(boundT,
m_XT[iCL][iLR][iA][iT][0],
m_XT[iCL][iLR][iA][iT][1],
1508 m_XT[iCL][iLR][iA][iT][2],
m_XT[iCL][iLR][iA][iT][3],
m_XT[iCL][iLR][iA][iT][4],
m_XT[iCL][iLR][iA][iT][5]);
1510 m_XT[iCL][iLR][iA][iT][np] =
1511 m_XT[iCL][iLR][iA][iT][0] + boundT
1512 * (
m_XT[iCL][iLR][iA][iT][1] + boundT
1513 * (
m_XT[iCL][iLR][iA][iT][2] + boundT
1514 * (
m_XT[iCL][iLR][iA][iT][3] + boundT
1515 * (
m_XT[iCL][iLR][iA][iT][4] + boundT
1516 * (
m_XT[iCL][iLR][iA][iT][5])))));
1551 for (
unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1552 for (
unsigned short iLR = 0; iLR < 2; ++iLR) {
1555 const std::vector<float> params = (*m_sResolFromDB)->getSigmaParams(iCL, iLR, iA, iT);
1556 unsigned short np = params.size();
1558 for (
unsigned short i = 0; i < np; ++i) {
1572 unsigned short groupId = (*m_fFactorFromDB)->getGroupID();
1573 unsigned short nEnt = (*m_fFactorFromDB)->getEntries();
1574 B2DEBUG(29,
"setFFactor called: groupId,nEnt= " << groupId <<
" " << nEnt);
1578 B2FATAL(
"CDCGeometryPar:: Invalid group-id " << groupId <<
" specified!");
1581 for (
unsigned short id = 0;
id < nEnt; ++id) {
1582 unsigned short np = ((*m_fFactorFromDB)->getFactors(
id)).size();
1583 if (np != 3) B2FATAL(
"CDCGeometryPar:: No. of fudge factors != 3!");
1584 for (
unsigned short i = 0; i < np; ++i) {
1601 const unsigned short isl = cm.getISuperLayer();
1603 const uint il = cm.getILayer();
1604 const int iw = cm.getIWire();
1605 const int iBd = cm.getBoardID();
1606 const WireID wID(isl, il, iw);
1607 m_wireToBoard.insert(pair<WireID, unsigned short>(wID, iBd));
1608 const int iCh = cm.getBoardChannel();
1617 unsigned short groupId = (*m_eDepToADCConversionsFromDB)->getGroupID();
1618 unsigned short nEnt = (*m_eDepToADCConversionsFromDB)->getEntries();
1620 if (nEnt > c_nSuperLayers) B2FATAL(
"CDCGeometryPar:: group-id " << groupId <<
" and #entries " << nEnt <<
" are inconsistent!");
1621 }
else if (groupId == 1) {
1622 if (nEnt > c_maxNSenseLayers) B2FATAL(
"CDCGeometryPar:: group-id " << groupId <<
" and #entries " << nEnt <<
" are inconsistent!");
1624 B2FATAL(
"CDCGeometryPar:: Invalid group-id " << groupId <<
" specified !");
1627 unsigned short cLMin[c_nSuperLayers], cLMax[c_nSuperLayers];
1630 for (
unsigned int sl = 1; sl < c_nSuperLayers; ++sl) {
1631 cLMin[sl] = cLMax[0] + 6 * sl - 5;
1632 cLMax[sl] = cLMax[0] + 6 * sl;
1635 for (
unsigned short id = 0;
id < nEnt; ++id) {
1636 unsigned short np = ((*m_eDepToADCConversionsFromDB)->getParams(
id)).size();
1637 if (np > 7) B2FATAL(
"CDCGeometryPar:: No. of edep-to-ADC conversion params. > 7");
1639 for (
unsigned short cL = cLMin[
id]; cL <= cLMax[id]; ++cL) {
1640 for (
unsigned short cell = 0; cell <
m_nWires[cL]; ++cell) {
1641 for (
unsigned short i = 0; i < np; ++i) {
1646 }
else if (groupId == 1) {
1647 for (
unsigned short cell = 0; cell <
m_nWires[id]; ++cell) {
1648 for (
unsigned short i = 0; i < np; ++i) {
1652 }
else if (groupId == 2) {
1654 B2FATAL(
"CDCGeometryPar::setEDepToADCConversions(): groupId=2 not ready!");
1673 const double cth = fabs(costh) + dlt;
1674 const double iGen = edep / dx;
1675 const double tmp = cth - gam * iGen;
1676 const double disc = tmp * tmp + 4.*alf * cth * iGen;
1680 iMea = cth * iGen / tmp;
1681 }
else if (disc >= 0.) {
1682 iMea = (-tmp +
sqrt(disc)) / (2.*alf);
1690 double convF = mainF;
1692 convF = mainF * std::min(iMea / iGen, 1.);
1701 convF *= 1. + a * (costh - b);
1721 if (set == c_Misaligned) {
1725 }
else if (set == c_Base) {
1747 if (set == c_Misaligned) {
1750 }
else if (set == c_Base) {
1769 if (set == c_Misaligned) {
1773 }
else if (set == c_Base) {
1795 if (set == c_Misaligned) {
1798 }
else if (set == c_Base) {
1809 if (set == c_Misaligned) {
1811 }
else if (set == c_Aligned) {
1819 static double IRWL[c_maxNSenseLayers] = {0};
1831 static double ORWL[c_maxNSenseLayers] = {0};
1847 const unsigned nWires =
m_nWires[layerId];
1851 const double phiSize = 2 * M_PI / double(nWires);
1868 for (
unsigned i = 0; i < 1; ++i) {
1869 const double phiF = phiSize * (double(i) +
offset)
1876 const double beta = (position.
Z() - b.Z()) / u.Z();
1878 double dPhi = std::atan2(position.
Y(), position.
X())
1879 - std::atan2(p.Y(), p.X())
1881 while (dPhi < 0) dPhi += (2. * M_PI);
1882 j = int(dPhi / phiSize);
1883 while (j >= nWires) j -= nWires;
1892 std::ofstream ofs(of.c_str(), std::ios::out);
1894 B2ERROR(
"CDCGeometryPar::read !!! can not open file : "
1897 ofs <<
"<?xml version=\"1.0\" encoding=\"UTF-8\"?>"
1899 <<
"<Subdetector type=\"CDC\">"
1901 <<
" <Name>CDC BelleII </Name>"
1903 <<
" <Description>CDC geometry parameters</Description>"
1905 <<
" <Version>0</Version>"
1907 <<
" <GeoCreator>CDCBelleII</GeoCreator>"
1911 <<
" <Rotation desc=\"Rotation of the whole cdc detector (should be the same as beampipe)\" unit=\"mrad\">0.0</Rotation>"
1913 <<
" <OffsetZ desc=\"The offset of the whole cdc in z with respect to the IP (should be the same as beampipe)\" unit=\"mm\">0.0</OffsetZ>"
1915 <<
" <Material>CDCGas</Material>"
1919 ofs <<
" <SLayers>" << endl;
1922 ofs <<
" <SLayer id=\"" << i <<
"\">" << endl;
1923 ofs <<
" <Radius desc=\"Radius of wires in this layer\" unit=\"mm\">" <<
senseWireR(i) <<
"</Radius>" << endl;
1924 ofs <<
" <BackwardZ desc=\"z position of this wire layer at backward endplate\" unit=\"mm\">" <<
senseWireBZ(
1925 i) <<
"</BackwardZ>" << endl;
1926 ofs <<
" <ForwardZ desc=\"z position of this wire layer at forward endplate\" unit=\"mm\">" <<
senseWireFZ(
1927 i) <<
"</ForwardZ>" << endl;
1930 ofs <<
" <NHoles desc=\"the number of holes in this layer, 2*(cell number)\">" <<
nWiresInLayer(
1931 i) * 2 <<
"</NHoles>" << endl;
1932 ofs <<
" <NShift desc=\"the shifted hole number of each wire in this layer\">" <<
nShifts(i) <<
"</NShift>" << endl;
1933 ofs <<
" <Offset desc=\"wire offset in phi direction at endplate\">" <<
m_offSet[i] <<
"</Offset>" << endl;
1934 ofs <<
" </SLayer>" << endl;
1937 ofs <<
" </SLayers>" << endl;
1938 ofs <<
" <FLayers>" << endl;
1941 ofs <<
" <FLayer id=\"" << i <<
"\">" << endl;
1942 ofs <<
" <Radius desc=\"Radius of field wires in this layer\" unit=\"mm\">" <<
fieldWireR(i) <<
"</Radius>" << endl;
1943 ofs <<
" <BackwardZ desc=\"z position of this field wire layer at backward endplate\" unit=\"mm\">" <<
fieldWireBZ(
1944 i) <<
"</BackwardZ>" << endl;
1945 ofs <<
" <ForwardZ desc=\"z position of this field wire layer at forward endplate\" unit=\"mm\">" <<
fieldWireFZ(
1946 i) <<
"</ForwardZ>" << endl;
1947 ofs <<
" </FLayer>" << endl;
1950 ofs <<
" </FLayers>" << endl;
1952 ofs <<
" <InnerWall name=\"InnerWall\">" << endl;
1953 ofs <<
" <InnerR desc=\"Inner radius\" unit=\"mm\">" <<
innerRadiusInnerWall() <<
"</InnerR>" << endl;
1954 ofs <<
" <OuterR desc=\"Outer radius\" unit=\"mm\">" <<
outerRadiusInnerWall() <<
"</OuterR>" << endl;
1955 ofs <<
" <BackwardZ desc=\"z position at backward endplate\" unit=\"mm\">" <<
m_zWall[0][0] <<
"</BackwardZ>" << endl;
1956 ofs <<
" <ForwardZ desc=\"z position at forward endplate\" unit=\"mm\">" <<
m_zWall[0][1] <<
"</ForwardZ>" << endl;
1957 ofs <<
" </InnerWall>" << endl;
1959 ofs <<
" <OuterWall name=\"OuterWall\">" << endl;
1960 ofs <<
" <InnerR desc=\"Inner radius\" unit=\"mm\">" <<
innerRadiusOuterWall() <<
"</InnerR>" << endl;
1961 ofs <<
" <OuterR desc=\"Outer radius\" unit=\"mm\">" <<
outerRadiusOuterWall() <<
"</OuterR>" << endl;
1962 ofs <<
" <BackwardZ desc=\"z position at backward endplate\" unit=\"mm\">" <<
m_zWall[2][0] <<
"</BackwardZ>" << endl;
1963 ofs <<
" <ForwardZ desc=\"z position at forward endplate\" unit=\"mm\">" <<
m_zWall[2][1] <<
"</ForwardZ>" << endl;
1964 ofs <<
" </OuterWall>" << endl;
1966 ofs <<
" </Content>" << endl
1967 <<
"</Subdetector>" << endl;
1971 double& Yb_sag,
double& Yf_sag)
const
2002 if (set == c_Aligned) {
2016 }
else if (set == c_Misaligned) {
2030 }
else if (set == c_Base) {
2045 B2FATAL(
"CDCGeometryPar::getWireSagEffect: called with an invalid set: " <<
" " << set);
2048 const double dx = Xf - Xb;
2049 const double dy = Yf - Yb;
2050 const double dz = Zf - Zb;
2052 const double Zfp =
sqrt(dz * dz + dx * dx);
2053 const double Zp = (Z - Zb) * Zfp / dz;
2055 const double Y_sag = (Coef * (Zp - Zfp) + dy / Zfp) * Zp + Yb;
2056 const double dydz = (Coef * (2.*Zp - Zfp) * Zfp + dy) / dz;
2058 Yb_sag = Y_sag + dydz * (Zb - Z);
2059 Yf_sag = Y_sag + dydz * (Zf - Z);
2065 const unsigned L = layerID;
2066 const unsigned C = cellID;
2070 const double phiSize = 2 * M_PI / double(
m_nWires[L]);
2072 const double phiF = phiSize * (double(C) +
offset)
2085 for (
int i = 0; i < 3; ++i) {
2103 const unsigned L = layerID;
2104 const unsigned C = cellID;
2106 static bool first =
true;
2107 static ofstream ofs;
2110 ofs.open(
"alignment.dat");
2113 ofs << L <<
" " << C;
2115 ofs << setiosflags(ios::showpoint | ios::uppercase);
2117 for (
int i = 0; i < 3; ++i) ofs <<
" " << setw(15) << setprecision(8) <<
m_BWirPos[L][C][i];
2119 for (
int i = 0; i < 3; ++i) ofs <<
" " << setw(15) << setprecision(8) <<
m_FWirPos[L][C][i];
2120 ofs << setiosflags(ios::fixed);
2127 const double theta)
const
2129 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2137 double delta = time - minTime;
2143 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2146 unsigned short ial[2] = {0};
2147 unsigned short ilr[2] = {lro, lro};
2150 unsigned short ith[2] = {0};
2153 unsigned short jal(0), jlr(0), jth(0);
2157 double timep = delta < 0. ? minTime - delta : time;
2160 for (
unsigned k = 0; k < 4; ++k) {
2165 w = (1. - wal) * (1. - wth);
2166 }
else if (k == 1) {
2170 w = (1. - wal) * wth;
2171 }
else if (k == 2) {
2175 w = wal * (1. - wth);
2176 }
else if (k == 3) {
2183 double boundary =
m_XT[iCLayer][jlr][jal][jth][6];
2185 if (timep < boundary) {
2187 const double& c1 =
m_XT[iCLayer][jlr][jal][jth][1];
2188 const double& c2 =
m_XT[iCLayer][jlr][jal][jth][2];
2189 const double& c3 =
m_XT[iCLayer][jlr][jal][jth][3];
2190 const double& c4 =
m_XT[iCLayer][jlr][jal][jth][4];
2191 const double& c5 =
m_XT[iCLayer][jlr][jal][jth][5];
2192 dDdt += w * ROOT::Math::Chebyshev4(timep, c1 + 3.*c3 + 5.*c5, 4.*c2 + 8.*c4, 6.*c3 + 10.*c5, 8.*c4, 10.*c5);
2194 dDdt += w * (
m_XT[iCLayer][jlr][jal][jth][1] + timep
2195 * (2.*
m_XT[iCLayer][jlr][jal][jth][2] + timep
2196 * (3.*
m_XT[iCLayer][jlr][jal][jth][3] + timep
2197 * (4.*
m_XT[iCLayer][jlr][jal][jth][4] + timep
2198 * (5.*
m_XT[iCLayer][jlr][jal][jth][5])))));
2201 dDdt += w *
m_XT[iCLayer][jlr][jal][jth][7];
2215 const double theta)
const
2217 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2229 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2232 unsigned short ial[2] = {0};
2233 unsigned short ilr[2] = {lro, lro};
2236 unsigned short ith[2] = {0};
2239 unsigned short jal(0), jlr(0), jth(0);
2243 double timep = time;
2247 for (
unsigned k = 0; k < 4; ++k) {
2252 w = (1. - wal) * (1. - wth);
2253 }
else if (k == 1) {
2257 w = (1. - wal) * wth;
2258 }
else if (k == 2) {
2262 w = wal * (1. - wth);
2263 }
else if (k == 3) {
2270 double boundary =
m_XT[iCLayer][jlr][jal][jth][6];
2272 if (timep < boundary) {
2274 dist += w * ROOT::Math::Chebyshev5(timep,
m_XT[iCLayer][jlr][jal][jth][0],
m_XT[iCLayer][jlr][jal][jth][1],
2275 m_XT[iCLayer][jlr][jal][jth][2],
m_XT[iCLayer][jlr][jal][jth][3],
m_XT[iCLayer][jlr][jal][jth][4],
m_XT[iCLayer][jlr][jal][jth][5]);
2277 dist += w * (
m_XT[iCLayer][jlr][jal][jth][0] + timep
2278 * (
m_XT[iCLayer][jlr][jal][jth][1] + timep
2279 * (
m_XT[iCLayer][jlr][jal][jth][2] + timep
2280 * (
m_XT[iCLayer][jlr][jal][jth][3] + timep
2281 * (
m_XT[iCLayer][jlr][jal][jth][4] + timep
2282 * (
m_XT[iCLayer][jlr][jal][jth][5]))))));
2285 dist += w * (
m_XT[iCLayer][jlr][jal][jth][7] * (timep - boundary) +
m_XT[iCLayer][jlr][jal][jth][8]);
2298 const bool calculateMinTime,
2299 const double inputMinTime)
const
2301 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2308 double minTime = calculateMinTime ?
getMinDriftTime(iCLayer, lr, alpha, theta) : inputMinTime;
2309 double delta = time - minTime;
2317 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2320 unsigned short ial[2] = {0};
2321 unsigned short ilr[2] = {lro, lro};
2324 unsigned short ith[2] = {0};
2327 unsigned short jal(0), jlr(0), jth(0);
2331 double timep = delta < 0. ? minTime - delta : time;
2336 for (
unsigned k = 0; k < 4; ++k) {
2341 w = (1. - wal) * (1. - wth);
2342 }
else if (k == 1) {
2346 w = (1. - wal) * wth;
2347 }
else if (k == 2) {
2351 w = wal * (1. - wth);
2352 }
else if (k == 3) {
2372 double boundary =
m_XT[iCLayer][jlr][jal][jth][6];
2374 if (timep < boundary) {
2376 dist += w * ROOT::Math::Chebyshev5(timep,
m_XT[iCLayer][jlr][jal][jth][0],
m_XT[iCLayer][jlr][jal][jth][1],
2377 m_XT[iCLayer][jlr][jal][jth][2],
m_XT[iCLayer][jlr][jal][jth][3],
m_XT[iCLayer][jlr][jal][jth][4],
m_XT[iCLayer][jlr][jal][jth][5]);
2379 dist += w * (
m_XT[iCLayer][jlr][jal][jth][0] + timep
2380 * (
m_XT[iCLayer][jlr][jal][jth][1] + timep
2381 * (
m_XT[iCLayer][jlr][jal][jth][2] + timep
2382 * (
m_XT[iCLayer][jlr][jal][jth][3] + timep
2383 * (
m_XT[iCLayer][jlr][jal][jth][4] + timep
2384 * (
m_XT[iCLayer][jlr][jal][jth][5]))))));
2387 dist += w * (
m_XT[iCLayer][jlr][jal][jth][7] * (timep - boundary) +
m_XT[iCLayer][jlr][jal][jth][8]);
2394 if (delta < 0.) dist *= -1.;
2400 const double theta)
const
2402 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2406 double minTime = 0.;
2412 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2415 unsigned short ial[2] = {0};
2416 unsigned short ilr[2] = {lro, lro};
2419 unsigned short ith[2] = {0};
2422 unsigned short jal(0), jlr(0), jth(0);
2425 double c[6] = {0.}, a[6] = {0.};
2426 for (
unsigned k = 0; k < 4; ++k) {
2431 w = (1. - wal) * (1. - wth);
2432 }
else if (k == 1) {
2436 w = (1. - wal) * wth;
2437 }
else if (k == 2) {
2441 w = wal * (1. - wth);
2442 }
else if (k == 3) {
2449 for (
int i = 0; i < 5; ++i) {
2450 c[i] += w *
m_XT[iCLayer][jlr][jal][jth][i];
2455 a[0] = c[0] - c[2] + c[4];
2456 a[1] = c[1] - 3.*c[3] + 5.*c[5];
2457 a[2] = 2.*c[2] - 8.*c[4];
2458 a[3] = 4.*c[3] - 20.*c[5];
2462 for (
int i = 0; i < 5; ++i) a[i] = c[i];
2469 const double det = a[1] * a[1] - 4.*a[2] * a[0];
2472 minTime = (-a[1] +
sqrt(det)) / (2.*a[2]);
2476 minTime = -a[1] / (2.*a[2]);
2479 }
else if (a[1] != 0.) {
2480 minTime = -a[0] / a[1];
2482 B2WARNING(
"CDCGeometryPar::getMinDriftTime: minDriftTime not determined; assume zero.\n" <<
"layer(#0-55),lr,alpha(rad),theta= " <<
2483 iCLayer <<
" " << lr <<
" " << alpha <<
" " << theta);
2492 const double epsi4x = 5.e-6;
2494 const unsigned short maxIter = 8;
2495 const double maxDt = 20.;
2496 unsigned short nIter = 0;
2497 double minXsq = 1.e10;
2498 double minMinTime = minTime;
2501 for (nIter = 0; nIter <= maxIter; ++nIter) {
2504 double x = a[0] + t * (a[1] + t * (a[2] + t * (a[3] + t * (a[4] + t * a[5]))));
2510 double xp = a[1] + t * (2 * a[2] + t * (3 * a[3] + t * (4 * a[4] + t * 5 * a[5])));
2511 double xpp = 2 * a[2] + t * (6 * a[3] + t * (12 * a[4] + t * 20 * a[5]));
2512 double den = xp * xp + x * xpp;
2519 edm = fabs(x * xp) /
sqrt(den);
2520 if (edm < epsi4x)
break;
2527 dt = std::min(dt, maxDt);
2529 dt = std::max(dt, -maxDt);
2532 B2WARNING(
"CDCGeometryPar::getMinDriftTime: den = 0\n" <<
"layer(#0-55),lr,alpha(rad),theta= " <<
2535 " " << alpha <<
" " << theta);
2541 if (nIter == (maxIter + 1)) minTime = minMinTime;
2594 const double theta)
const
2596 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2602 const double eps = 2.5e-1;
2603 const double maxTrials = 100;
2611 double maxTime = 2000.;
2617 double t0 = minTime;
2619 const bool calMinTime =
false;
2624 double t1 = maxTime;
2626 while (((t1 - t0) > eps) && (i < maxTrials)) {
2627 time = 0.5 * (t0 + t1);
2628 double d1 =
getDriftLength(time, iCLayer, lr, alpha, theta, calMinTime, minTime) - dist;
2638 if (i >= maxTrials - 1 || time > maxTime) {
2639 B2WARNING(
"CDCGeometryPar::getDriftTime " << dist <<
" " << iCLayer <<
" " << alpha <<
" " << lr <<
" " << t0 <<
" " << t1 <<
" " <<
2650 const double theta)
const
2652 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2659 const double driftL = fabs(DriftL0);
2665 B2FATAL(
"linearInterpolationOfSgm = false is not allowed now !");
2669 unsigned short ial[2] = {0};
2670 unsigned short ilr[2] = {lro, lro};
2673 unsigned short ith[2] = {0};
2677 unsigned short jal(0), jlr(0), jth(0);
2679 for (
unsigned k = 0; k < 4; ++k) {
2684 w = (1. - wal) * (1. - wth);
2685 }
else if (k == 1) {
2689 w = (1. - wal) * wth;
2690 }
else if (k == 2) {
2694 w = wal * (1. - wth);
2695 }
else if (k == 3) {
2713 const double& P0 =
m_Sigma[iCLayer][jlr][jal][jth][0];
2714 const double& P1 =
m_Sigma[iCLayer][jlr][jal][jth][1];
2715 const double& P2 =
m_Sigma[iCLayer][jlr][jal][jth][2];
2716 const double& P3 =
m_Sigma[iCLayer][jlr][jal][jth][3];
2717 const double& P4 =
m_Sigma[iCLayer][jlr][jal][jth][4];
2718 const double& P5 =
m_Sigma[iCLayer][jlr][jal][jth][5];
2719 const double& P6 =
m_Sigma[iCLayer][jlr][jal][jth][6];
2721#if defined(CDC_DEBUG)
2722 cout <<
"driftL= " << driftL << endl;
2723 cout <<
"iCLayer= " << iCLayer <<
" " << jlr <<
" " << jal <<
" " << jth << endl;
2724 cout <<
"P0= " << P0 << endl;
2725 cout <<
"P1= " << P1 << endl;
2726 cout <<
"P2= " << P2 << endl;
2727 cout <<
"P3= " << P3 << endl;
2728 cout <<
"P4= " << P4 << endl;
2729 cout <<
"P5= " << P5 << endl;
2730 cout <<
"P6= " << P6 << endl;
2735 sigma += w *
sqrt(P0 / (driftL * driftL + P1) + P2 * driftL + P3 +
2736 P4 * exp(P5 * (driftL - P6) * (driftL - P6)));
2738 double forthTermAtP7 = P4 * exp(P5 * (P7 - P6) * (P7 - P6));
2739 const double& P8 =
m_Sigma[iCLayer][jlr][jal][jth][8];
2741 double sigmaAtP7 =
sqrt(P0 / (P7 * P7 + P1) + P2 * P7 + P3 + forthTermAtP7);
2742 sigma += w * (P8 * (driftL - P7) + sigmaAtP7);
2744 double onePls4AtP7 =
sqrt(P0 / (P7 * P7 + P1) + forthTermAtP7);
2745 const double onePls4 = P8 * (driftL - P7) + onePls4AtP7;
2746 sigma += w *
sqrt(P2 * driftL + P3 + onePls4 * onePls4);
2748 forthTermAtP7 =
sqrt(forthTermAtP7);
2749 const double forthTerm = P8 * (driftL - P7) + forthTermAtP7;
2750 sigma += w *
sqrt(P0 / (driftL * driftL + P1) + P2 * driftL + P3 +
2751 forthTerm * forthTerm);
2764 unsigned short lr = 0;
2765 double wCrossT = (posOnWire.
Cross(posOnTrack)).Z();
2769 }
else if (wCrossT > 0.) {
2772 if ((posOnTrack - posOnWire).Perp() != 0.) {
2773 double wCrossP = (posOnWire.
Cross(momentum)).Z();
2775 if (posOnTrack.
Perp() > posOnWire.
Perp()) {
2780 }
else if (wCrossP < 0.) {
2781 if (posOnTrack.
Perp() < posOnWire.
Perp()) {
2799 const double distanceCrossP = ((posOnWire - posOnTrack).Cross(momentum)).Z();
2800 unsigned short int lr = (distanceCrossP > 0.) ? 1 : 0;
2807 const double wx = posOnWire.
X();
2808 const double wy = posOnWire.
Y();
2809 const double px = momentum.X();
2810 const double py = momentum.Y();
2812 const double cross = wx * py - wy * px;
2813 const double dot = wx * px + wy * py;
2815 return atan2(cross,
dot);
2820 return atan2(momentum.Perp(), momentum.Z());
2826 unsigned short lro = (fabs(alpha) <= 0.5 * M_PI) ? lr : abs(lr - 1);
2834 double alphao = alpha;
2835 if (alpha > 0.5 * M_PI) {
2837 }
else if (alpha < -0.5 * M_PI) {
2847 double thetao = fabs(alpha) > 0.5 * M_PI ? M_PI - theta : theta;
2853 unsigned short lrs[2])
const
2862 lrs[0] = abs(lrs[0] - 1);
2869 lrs[1] = abs(lrs[1] - 1);
2876 points[0] = points[1] - 1;
2883 unsigned short lrs[2])
const
2892 lrs[0] = abs(lrs[0] - 1);
2899 lrs[1] = abs(lrs[1] - 1);
2906 points[0] = points[1] - 1;
2932 points[0] = points[1] - 1;
2939 unsigned short points[2])
const
2955 points[0] = points[1] - 1;
2987 const unsigned short nLayers[c_nSuperLayers] = {8, 6, 6, 6, 6, 6, 6, 6, 6};
2989 for (
unsigned short SLayer = 0; SLayer < c_nSuperLayers; ++SLayer) {
2990 unsigned short firstCLayer = 0;
2991 for (
unsigned short i = 0; i < SLayer; ++i) {
2992 firstCLayer += nLayers[i];
2997 for (
unsigned short Layer = 0; Layer < nLayers[SLayer]; ++Layer) {
2998 unsigned short CLayer = firstCLayer + Layer;
3000 if (CLayer == firstCLayer) {
3003 }
else if (CLayer == firstCLayer + 1) {
3009 if (Layer % 2 == 0) {
DataType Z() const
access variable Z (= .at(2) without boundary check)
void SetX(DataType x)
set X/1st-coordinate
B2Vector3< DataType > Cross(const B2Vector3< DataType > &p) const
Cross product.
DataType X() const
access variable X (= .at(0) without boundary check)
DataType Y() const
access variable Y (= .at(1) without boundary check)
void SetZ(DataType z)
set Z/3rd-coordinate
void SetY(DataType y)
set Y/2nd-coordinate
DataType Perp() const
The transverse component (R in cylindrical coordinate system).
B2Vector3< DataType > Unit() const
Unit vector parallel to this.
static const baseType layerDPhi
Layer rotation in global X-Y plane (gamma) dPhi = forward - backward endplate.
static const baseType layerDy
Layer shift in global Y dY = forward - backward endplate.
static const baseType wireBwdZ
Wire Z position w.r.t. nominal on backward endplate.
static const baseType layerDx
Layer shift in global X dX = forward - backward endplate.
static const baseType wireBwdY
Wire Y position w.r.t. nominal on backward endplate.
static const baseType wireFwdZ
Wire Z position w.r.t. nominal on forward endplate.
static const baseType wireFwdY
Wire Y position w.r.t. nominal on forward endplate.
static const baseType wireFwdX
Wire X position w.r.t. nominal on forward endplate.
static const baseType wireBwdX
Wire X position w.r.t. nominal on backward endplate.
static const baseType layerY
Layer shift in global Y at backward endplate.
static const baseType layerX
Layer shift in global X at backward endplate.
static const baseType layerPhi
Layer rotation in global X-Y plane (gamma) at backward endplate.
static const baseType wireTension
Wire tension w.r.t. nominal (=50. ?)
The Class for CDC geometry.
static const baseType wireBwdZ
Wire Z position w.r.t. nominal on backward endplate.
static const baseType wireBwdY
Wire Y position w.r.t. nominal on backward endplate.
static const baseType wireFwdZ
Wire Z position w.r.t. nominal on forward endplate.
static const baseType wireFwdY
Wire Y position w.r.t. nominal on forward endplate.
static const baseType wireFwdX
Wire X position w.r.t. nominal on forward endplate.
static const baseType wireBwdX
Wire X position w.r.t. nominal on backward endplate.
static const baseType wireTension
Wire tension w.r.t. nominal (=50. ?)
The Class for CDC Geometry Control Parameters.
bool getDebug() const
Get debug flag.
bool getSigmaInputType()
Get input type for sigma.
bool getMisalignmentInputType()
Get input type for wire misalignment.
std::string getT0File() const
Get input file name for t0.
bool getDisplacementInputType()
Get input type for wire displacement.
double getAddFudgeFactorForSigmaForMC() const
Get additional fudge factor for space resol for MC.
std::string getEDepToADCFile() const
Get input file name for edeptoadc.
std::string getDisplacementFile() const
Get input file name for wire displacement.
std::string getMisalignmentFile() const
Get input file name for wire misalignment.
std::string getAlignmentFile() const
Get input file name for wire alignment.
bool getAlignmentInputType()
Get input type for wire alignment.
double getMaterialDefinitionMode() const
Get material definition mode.
std::string getPropSpeedFile() const
Get input file name for prop-speed.
bool getT0InputType()
Get input type for t0.
bool getEDepToADCInputType()
Get input type for edeptoadc.
std::string getSigmaFile() const
Get input file name for sigma.
bool getAlignment() const
Get alignment switch.
bool getMisalignment() const
Get misalignment switch.
bool getDisplacement() const
Get displacement switch.
int getSenseWireZposMode() const
Get sense wire z position mode.
double getAddFudgeFactorForSigmaForData() const
Get additional fudge factor for space resol for data.
std::string getXtFile() const
Get input file name for xt-relation.
bool getTwInputType()
Get input type for time-walk.
std::string getFFactorFile() const
Get input file name for fudge factor.
bool getChMapInputType()
Get input type for channel map.
std::string getTwFile() const
Get input file name for time-walk.
bool getFFactorInputType()
Get input type for fuge factor.
bool getXtInputType()
Get input type for xt.
bool getBwInputType()
Get input type for bad wire.
bool getPropSpeedInputType()
Get input type for prop.
static CDCGeoControlPar & getInstance()
Static method to get a reference to the CDCGeoControlPar instance.
std::string getChMapFile() const
Get input file name for channel map.
The Class for CDC Geometry Parameters.
void outputDesignWirParam(unsigned layerID, unsigned cellID) const
Write the designed wire parameters to the alignment.dat (default).
std::map< WireID, unsigned short > m_wireToChannel
map relating wire-id and channel-id.
void setWirPosAlignParams()
Set wire alignment params.
int m_materialDefinitionMode
Control switch for gas and wire material definition.
unsigned short m_nAlphaPoints4Sgm
No.
void readSigma(const GearDir &gbxParams, int mode=0)
Read spatial resolution table.
void readTW(const GearDir &gbxParams, int mode=0)
Read time-walk parameter.
void setXtRel()
Set XT-relation table (from DB) (new).
float m_BWirPosMisalign[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
DBObjPtr< CDCBadWires > * m_badWireFromDB
bad-wires retrieved from DB.
double getTheta(const B2Vector3D &momentum) const
Returns track incident angle (theta in rad.).
void readT0(const GearDir &gbxParams, int mode=0)
Read t0 parameters (from a file).
void setT0()
Set t0 parameters (from DB)
float m_alphaPoints4Sgm[c_maxNAlphaPoints]
alpha sampling points for sigma (rad)
virtual ~CDCGeometryPar()
Destructor.
ushort m_maxNSuperLayers
Maximum number of Super Layers.
double m_fudgeFactorForSigma[3]
Fuge factor for space resol.
double outerRadiusInnerWall() const
Returns the outer radius of the inner wall.
int m_sigmaParamMode
Mode for sigma parameterization.
bool m_alignment
Switch for alignment.
float m_alphaPoints[c_maxNAlphaPoints]
alpha sampling points for xt (rad)
void getClosestAlphaPoints(const double alpha, double &wal, unsigned short points[2], unsigned short lrs[2]) const
Returns the two closest alpha points for the input track incident angle (alpha).
void readEDepToADC(const GearDir &gbxParams, int mode=0)
Read spatial edep-to-adc conv.
double m_globalPhiRotation
Global ratation in phi (rad.); only for sense wires now.
EWirePosition
Wire position set.
double innerRadiusOuterWall() const
Returns the inner radius of the outer wall.
unsigned cellId(unsigned layerId, const B2Vector3D &position) const
The method to get cell id based on given layer id and the position.
void newReadSigma(const GearDir &gbxParams, int mode=0)
Read spatial resolution table in new format.
void setEDepToADCConversions()
Set edep-to-ADC conversion params.
double m_nominalPropSpeed
Nominal propagation speed of the sense wire (27.25 cm/nsec).
int nShifts(int layerId) const
Returns number shift.
float m_thetaPoints[c_maxNThetaPoints]
theta sampling points for xt (rad)
void setDesignWirParam(unsigned layerID, unsigned cellID)
Set the desizend wire parameters.
void getWireSagEffect(EWirePosition set, unsigned layerID, unsigned cellID, double zw, double &ywb_sag, double &ywf_sag) const
Compute effects of the sense wire sag.
DBArray< CDCDisplacement > * m_displacementFromDB
displacement params.
bool m_XTetc
Switch for reading x-t etc.
void setShiftInSuperLayer()
Calculates and saves shifts in super-layers (to be used in searching hits in neighboring cells)
int m_nShifts[c_maxNSenseLayers]
The array to store shifted cell number in each sense wire layer.
bool m_wireSag
Switch for sense wire sag.
bool m_XTetc4Recon
Switch for selecting xt etc.
double getDriftLength0(double dt, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Return the drift dength to the sense wire; tentative ver.
unsigned short m_boardAndChannelToWire[c_nBoards][48]
array relating board-channel-id and wire-id.
unsigned short m_nThetaPoints4Sgm
No.
float m_WireSagCoef[c_maxNSenseLayers][c_maxNDriftCells]
Wire sag coefficient for each cell; ibid.
void generateXML(const std::string &of)
Generate an xml file used in gearbox.
void getClosestAlphaPoints4Sgm(const double alpha, double &wal, unsigned short points[2], unsigned short lrs[2]) const
Returns the two closest alpha points for sigma for the input track incident angle (alpha).
DBObjPtr< HardwareClockSettings > m_clockSettings
hardware clock settings
float m_eDepToADCParams[c_maxNSenseLayers][c_maxNDriftCells][7]
edep-to-ADC conv.
double m_minTrackLength
Minimum track length for G4 step.
double getAlpha(const B2Vector3D &posOnWire, const B2Vector3D &momentum) const
Returns track incident angle in rphi plane (alpha in rad.).
double fieldWireR(int layerId) const
Returns radius of field wire in each layer.
float m_FWirPosMisalign[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
const B2Vector3D wireForwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the forward position of the input sense wire.
bool isDeadWire(const WireID &wid, double &eff)
Inquire if the wire is dead.
unsigned m_nWires[c_maxNSenseLayers]
The array to store the wire number in each sense wire layre.
double getDriftV(double dt, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Get the realistic drift velocity.
double m_clockFreq4TDC
Clock frequency used for TDC (GHz).
int m_nFLayer
The number of field wire layer.
static CDCGeometryPar * m_B4CDCGeometryParDB
Pointer that saves the instance of this class.
DBObjPtr< CDCMisalignment > * m_misalignmentFromDB
misalignment params.
float m_XT[c_maxNSenseLayers][2][c_maxNAlphaPoints][c_maxNThetaPoints][c_nXTParams]
XT-relation coefficients for each layer, Left/Right, entrance angle and polar angle.
signed short m_shiftInSuperLayer[c_nSuperLayers][8]
shift in phi-direction wrt the 1st layer in each super layer
std::string m_version
The version of geometry parameters.
const B2Vector3D wireBackwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the backward position of the input sense wire.
bool m_linearInterpolationOfXT
Switch for linear interpolation of xt.
DBObjPtr< CDCSpaceResols > * m_sResolFromDB
sigma params.
void setDisplacement()
Set displacement of sense wire.
double getOutgoingAlpha(const double alpha) const
Converts incoming- to outgoing-alpha.
bool isHotWire(const WireID &wid)
Inquire if the wire is hot.
double m_meanT0
mean t0 over all wires; should be double.
bool m_debug
Switch for debug printing.
DBObjPtr< CDCTimeWalks > * m_timeWalkFromDB
time-walk coeffs.
ushort m_nSenseWires
Maximum number of Sense Wires.
unsigned short m_tdcOffset
Not used; to be removed later.
double getSigma(double dist, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Return the basic resolution of drift length (cm).
CDCGeometryPar(const CDCGeometry *=nullptr)
Singleton class.
unsigned short getOutgoingLR(const unsigned short lr, const double alpha) const
Converts incoming-lr to outgoing-lr.
ushort m_maxNSenseLayers
Maximum number of Sense Wire Layers.
void readWirePositionParams(EWirePosition set, const CDCGeometry *geom)
Read displacement or (mis)alignment params from text file.
float m_thetaPoints4Sgm[c_maxNThetaPoints]
theta sampling points for sigma (rad)
std::map< WireID, unsigned short > m_wireToBoard
map relating wire-id and board-id.
double getOutgoingTheta(const double alpha, const double theta) const
Converts incoming- to outgoing-theta.
double m_senseWireDiameter
The diameter of sense wires.
double offset(int layerID) const
Return wire offset in phi direction at endplate.
bool m_displacement
Switch for displacement.
int m_twParamMode
Mode for tw parameterization.
double senseWireBZ(int layerId) const
Returns backward z position of sense wire in each layer.
double m_dzSBackwardLayer[c_maxNSenseLayers]
Corrections for backward z position of sense wire layers.
double getMinDriftTime(unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI) const
Return the min.
DBObjPtr< CDCBadBoards > * m_badBoardsFromDB
bad-boards retrieved from DB.
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
double fieldWireBZ(int layerId) const
Returns backward z position of field wire in each layer.
void getClosestThetaPoints(const double alpha, const double theta, double &wth, unsigned short points[2]) const
Returns the two closest theta points for the input track incident angle (theta).
double senseWireFZ(int layerId) const
Returns forward z position of sense wire in each layer.
void readPropSpeed(const GearDir &gbxParams, int mode=0)
Read the propagation speed along the sense wire.
double innerRadiusInnerWall() const
Returns the inner radius of the inner wall.
DBObjPtr< CDCTimeZeros > * m_t0FromDB
t0s retrieved from DB.
int m_xtFileFormat
Format of xt input file.
int m_sigmaFileFormat
Format of sigma input file.
ushort m_firstSuperLayerOffset
Offset of the first super layer (for reduced CDC studies)
void setChMap()
Set channel map (from DB)
DBObjPtr< CDCFudgeFactorsForSigma > * m_fFactorFromDB
fudge factors retrieved from DB.
int m_nSLayer
The number of sense wire layer.
DBArray< CDCChannelMap > * m_chMapFromDB
channel map retrieved from DB.
void calcMeanT0(double minT0=3800, double maxT0=5800, int maxIt=10, double nStdv=3, double epsi=0.1)
Calculate mean t0 in ns (over all good wires)
void readFromDB(const CDCGeometry &)
Gets geometry parameters from database.
DBObjPtr< CDCXtRelations > * m_xtRelFromDB
xt params.
double m_nominalDriftV
Nominal drift velocity (4.0x10^-3 cm/nsec).
float m_FWirPos[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
float m_WireSagCoefAlign[c_maxNSenseLayers][c_maxNDriftCells]
Wire sag coefficient incl.
double m_zSForwardLayer[c_maxNSenseLayers]
The array to store forward z position of sense wire layers.
bool m_linearInterpolationOfSgm
Switch for linear interpolation of sigma.
DBObjPtr< CDCPropSpeeds > * m_propSpeedFromDB
prop.
void setBadWire()
Set bad-wires (from DB)
int m_xtParamMode
Mode for xt parameterization.
float m_FWirPosAlign[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
double fieldWireFZ(int layerId) const
Returns forward z position of field wire in each layer.
double m_dzSForwardLayer[c_maxNSenseLayers]
Corrections for forward z position of sense wire layers.
double getEDepToADCConvFactor(unsigned short layer, unsigned short cell, double edep, double dx, double costh)
Return edep-to-ADC conversion factor.
double m_nominalDriftVInv
Inverse of the nominal drift velocity.
const double * innerRadiusWireLayer() const
Returns an array of inner radius of wire layers.
double getDriftTime(double dist, unsigned short layer, unsigned short lr, double alpha, double theta) const
Return the drift time to the sense wire.
void setPropSpeed()
Set prop.
signed short getShiftInSuperLayer(unsigned short iSuperLayer, unsigned short iLayer) const
Returns shift in the super-layer.
void setBadBoard()
Set bad-boards (from DB)
double getWireSagCoef(EWirePosition set, uint layerId, int cellId) const
Returns coefficient for the sense wire sag.
void readChMap()
Read channel map between wire-id and electronics-id.
double m_zFForwardLayer[c_maxNFieldLayers]
The array to store forward z position of field wire layers.
double m_rSLayer[c_maxNSenseLayers]
The array to store radius of sense wire layers.
double m_fieldWireDiameter
The diameter of field wires.
double m_zFBackwardLayer[c_maxNFieldLayers]
The array to store backward z position of field wire layers.
ushort m_nFieldWires
Maximum number of Field Wires.
float m_WireSagCoefMisalign[c_maxNSenseLayers][c_maxNDriftCells]
Wire sag coefficient incl.
double getDriftLength(double dt, unsigned short layer, unsigned short lr, double alpha=0., double theta=0.5 *M_PI, bool calculateMinTime=true, double minTime=0.) const
Return the drift dength to the sense wire.
double m_cellSize[c_maxNSenseLayers]
The array to store cell size in each sense wire layer.
float m_BWirPosAlign[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
double outerRadiusOuterWall() const
Returns the outer radius of the outer wall.
float m_t0[c_maxNSenseLayers][c_maxNDriftCells]
t0 for each sense-wire (in nsec).
void setWirPosMisalignParams()
Set wire misalignment params.
double m_offSet[c_maxNSenseLayers]
The array to store z offset of sense wire layers.
void setTW()
Set time-walk parameters.
void setSResol()
Set spatial resolution (from DB).
double m_zSBackwardLayer[c_maxNSenseLayers]
The array to store backward z position of sense wire layers.
unsigned short getOldLeftRight(const B2Vector3D &posOnWire, const B2Vector3D &posOnTrack, const B2Vector3D &momentum) const
Returns old left/right.
float m_propSpeedInv[c_maxNSenseLayers]
Inverse of propagation speed of the sense wire.
ushort m_maxNCellsPerLayer
Maximum number wires within a layer.
double m_zWall[4][2]
The array to store z position of inner wall and outer wall.
bool isBadWire(const WireID &wid)
Inquire if the wire is totally-dead.
DBObjPtr< CDCAlignment > * m_alignmentFromDB
alignment params.
ushort m_maxNFieldLayers
Maximum number of Field Wire Layers.
bool m_misalignment
Switch for misalignment.
double m_senseWireDensity
The density of sense wires.
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
ushort m_firstLayerOffset
Offset of the first layer (for reduced CDC studies)
void Print() const
Print some debug information.
float m_timeWalkCoef[c_nBoards][2]
coefficients for time walk.
void readFFactor(const GearDir &gbxParams, int mode=0)
Read fudge factors.
unsigned short m_nThetaPoints
No.
double m_maxSpaceResol
max space resolution allowed (cm).
double m_thresholdEnergyDeposit
Energy thresh.
float m_BWirPos[c_maxNSenseLayers][c_maxNDriftCells][3]
Wire position incl.
float m_Sigma[c_maxNSenseLayers][2][c_maxNAlphaPoints][c_maxNThetaPoints][c_nSigmaParams]
position resolution for each layer.
void readXT(const GearDir &gbxParams, int mode=0)
Read XT-relation table.
double m_rFLayer[c_maxNFieldLayers]
The array to store radius of field wire layers.
void newReadXT(const GearDir &gbxParams, int mode=0)
Read XT-relation table in new format.
double m_nominalSpaceResol
Nominal spatial resolution (0.0130 cm).
unsigned nWireLayers() const
Returns a number of wire layers.
unsigned short m_nAlphaPoints
No.
unsigned short getNewLeftRightRaw(const B2Vector3D &posOnWire, const B2Vector3D &posOnTrack, const B2Vector3D &momentum) const
Returns new left/right_raw.
int m_senseWireZposMode
Mode for sense wire z position corr.
double m_rWall[4]
The array to store radius of inner wall and outer wall.
double m_tdcBinWidth
TDC bin width (nsec/bin).
bool m_modLeftRightFlag
Switch for modified left/right flag.
void getClosestThetaPoints4Sgm(const double alpha, const double theta, double &wth, unsigned short points[2]) const
Returns the two closest theta points for sigma for the input track incident angle (theta).
DBObjPtr< CDCEDepToADCConversions > * m_eDepToADCConversionsFromDB
Pointer to edep-to-ADC conv.
const double * outerRadiusWireLayer() const
Returns an array of outer radius of wire layers.
double senseWireR(int layerId) const
Returns radius of sense wire in each layer.
double m_senseWireTension
The tension of sense wires.
void setFFactor()
Set fudge factors (from DB).
bool getModLeftRightFlag() const
Get modified left/right flag.
double getMinTrackLength() const
Get minimum track length.
double getThresholdEnergyDeposit() const
Get threshold for Energy Deposit;.
static CDCSimControlPar & getInstance()
Static method to get a reference to the CDCSimControlPar instance.
bool getWireSag() const
Get wiresag flag.
Class for accessing arrays of objects in the database.
Class for accessing objects in the database.
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
GearDir is the basic class used for accessing the parameter store.
virtual std::string getString(const std::string &path="") const noexcept(false) override
Get the parameter path as a string.
Optional DBObjPtr: This class behaves the same as the DBObjPtr except that it will not raise errors w...
Class to identify a wire inside the CDC.
unsigned short getICLayer() const
Getter for continuous layer numbering.
unsigned short getIWire() const
Getter for wire within the layer.
unsigned short getEWire() const
Getter for encoded wire number.
B2Vector3< double > B2Vector3D
typedef for common usage with double
double sqrt(double a)
sqrt for double
T dot(GeneralVector< T > a, GeneralVector< T > b)
dot product of two general vectors
void openFileB(boost::iostreams::filtering_istream &ifs, const std::string &fileName0)
Open a file using boost (to be able to read a gzipped file)
void openFileA(std::ifstream &ifs, const std::string &fileName0)
Open a file.
Abstract base class for different kinds of events.