9#include <framework/gearbox/GearDir.h>
10#include <framework/logging/Logger.h>
11#include <framework/utilities/FileSystem.h>
13#include <cdc/geometry/CDCGeometryPar.h>
14#include <cdc/geometry/CDCGeoControlPar.h>
15#include <cdc/simulation/CDCSimControlPar.h>
16#include <cdc/utilities/OpenFile.h>
21#include <boost/iostreams/filtering_stream.hpp>
23#include <Math/ChebyshevPol.h>
44 if ((*m_t0FromDB).isValid()) {
51 if ((*m_badWireFromDB).isValid()) {
57 if ((*m_badBoardsFromDB).isValid()) {
63 if ((*m_propSpeedFromDB).isValid()) {
70 if ((*m_timeWalkFromDB).isValid()) {
77 if ((*m_xtRelFromDB).isValid()) {
84 if ((*m_sResolFromDB).isValid()) {
91 if ((*m_fFactorFromDB).isValid()) {
98 if ((*m_chMapFromDB).isValid()) {
105 if ((*m_displacementFromDB).isValid()) {
112 if ((*m_alignmentFromDB).isValid()) {
120 if ((*m_misalignmentFromDB).isValid()) {
129 if ((*m_eDepToADCConversionsFromDB).isValid()) {
142 B2WARNING(
"CDCGeometryPar: Strange that readFromDB is not called! Please make sure that CDC is included in Geometry.");
181 for (
unsigned i = 0; i < 4; ++i) {
183 for (
unsigned j = 0; j < 2; ++j)
186 for (
unsigned i = 0; i < c_maxNSenseLayers; ++i) {
198 for (
unsigned i = 0; i < c_maxNFieldLayers; ++i) {
204 for (
unsigned L = 0; L < c_maxNSenseLayers; ++L) {
205 for (
unsigned C = 0; C < c_maxNDriftCells; ++C) {
206 for (
unsigned i = 0; i < 3; ++i) {
214 for (
unsigned i = 0; i < 7; ++i) {
224 for (
unsigned L = 0; L < c_maxNSenseLayers; ++L) {
225 for (
unsigned i = 0; i < 2; ++i) {
226 for (
unsigned alpha = 0; alpha < c_maxNAlphaPoints; ++alpha) {
227 for (
unsigned theta = 0; theta < c_maxNThetaPoints; ++theta) {
228 for (
unsigned xtparam = 0; xtparam < c_nXTParams; ++xtparam) {
229 m_XT[L][i][alpha][theta][xtparam] = 0.;
232 for (
unsigned sigmaparam = 0; sigmaparam < c_nSigmaParams; ++sigmaparam) {
233 m_Sigma[L][i][alpha][theta][sigmaparam] = 0.;
240 for (
unsigned board = 0; board < c_nBoards; ++board) {
241 for (
unsigned i = 0; i < 2; ++i) {
244 for (
unsigned channel = 0; channel < 48; ++channel) {
249 for (
unsigned superLayer = 0; superLayer < c_nSuperLayers; ++superLayer) {
250 for (
unsigned layer = 0; layer < 8; ++layer) {
263 m_rWall[0] = geom.getInnerWall(2).getRmin();
264 m_zWall[0][0] = geom.getInnerWall(0).getZbwd();
265 m_zWall[0][1] = geom.getInnerWall(0).getZfwd();
267 m_rWall[1] = geom.getInnerWall(0).getRmax();
268 m_zWall[1][0] = geom.getInnerWall(0).getZbwd();
269 m_zWall[1][1] = geom.getInnerWall(0).getZbwd();
272 m_rWall[2] = geom.getOuterWall(0).getRmin();
273 m_zWall[2][0] = geom.getOuterWall(0).getZbwd();
274 m_zWall[2][1] = geom.getOuterWall(0).getZfwd();
276 m_rWall[3] = geom.getOuterWall(1).getRmax();
277 m_zWall[3][0] = geom.getOuterWall(0).getZbwd();
278 m_zWall[3][1] = geom.getOuterWall(0).getZfwd();
287 B2DEBUG(100,
"CDCGeometryPar: Define a mixture of gases and wires in the tracking volume.");
290 B2FATAL(
"CDCGeometryPar: Materialdefinition=2 is disabled for now.");
292 B2FATAL(
"CDCGeometryPar: Materialdefinition mode you specify is invalid.");
303 GearDir content =
GearDir(
"/Detector/DetectorComponent[@name=\"CDC\"]/Content/");
311 for (
const auto& sense : geom.getSenseLayers()) {
312 uint layerId = sense.getId();
317 m_nWires[layerId] = sense.getNWires();
319 m_offSet[layerId] = sense.getOffset();
332 B2FATAL(
"CDCGeometryPar: invalid wire z definition mode specified");
336 const int nWires =
m_nWires[layerId];
337 for (
int iCell = 0; iCell < nWires; ++iCell) {
344 for (
const auto& field : geom.getFieldLayers()) {
345 uint layerId = field.getId();
377 B2FATAL(
"HardwareClockSettings payloads are not valid.");
378 const double officialClockFreq4TDC = 2 *
m_clockSettings->getAcceleratorRF();
380 B2WARNING(
"ClockFreq4TDC changed from cdclocal " << scientific << setprecision(6) <<
m_clockFreq4TDC <<
" to official " <<
381 officialClockFreq4TDC <<
" (GHz) (difference larger than 0.01%)");
384 B2DEBUG(100,
"CDCGeometryPar: Clock freq. for TDC= " <<
m_clockFreq4TDC <<
" (GHz).");
386 B2DEBUG(100,
"CDCGeometryPar: TDC bin width= " <<
m_tdcBinWidth <<
" (ns).");
397 B2DEBUG(100,
"CDCGeometryPar: Load displacement params. (=1); not load (=0):" <<
m_displacement);
400 B2DEBUG(100,
"CDCGeometryPar: Read displacement from DB");
409 B2DEBUG(100,
"CDCGeometryPar: Load alignment params. (=1); not load (=0):" <<
413 B2DEBUG(100,
"CDCGeometryPar: Read alignment from DB");
422 B2DEBUG(100,
"CDCGeometryPar: Load misalignment params. (=1); not load (=0):" <<
426 B2DEBUG(100,
"CDCGeometryPar: Read misalignment from DB");
439 B2FATAL(
"ModifiedLeftRightFlag = true is disabled for now; need to update a G4-related code in framework...");
448 B2DEBUG(100,
"CDCGeometryPar: Read xt from DB");
455 B2DEBUG(100,
"CDCGeometryPar: Read sigma from DB");
462 B2DEBUG(100,
"CDCGeometryPar: Read fudge factors from DB");
469 B2DEBUG(100,
"CDCGeometryPar: Read prop-speed from DB");
476 B2DEBUG(100,
"CDCGeometryPar: Read t0 from DB");
483 B2DEBUG(100,
"CDCGeometryPar: Read badwire from DB");
487 B2FATAL(
"Text file input mode for bdwires is disabled now!");
491 B2DEBUG(100,
"CDCGeometryPar: Read ch-map from DB");
498 B2DEBUG(100,
"CDCGeometryPar: Read time-walk from DB");
503 B2DEBUG(100,
"CDCGeometryPar: Time-walk param. mode= " <<
m_twParamMode);
506 B2DEBUG(29,
"CDCGeometryPar: Read EDepToADC from DB");
507 if ((*m_eDepToADCConversionsFromDB).isValid()) {
579 std::string fileName0;
584 }
else if (set == c_Misaligned) {
586 }
else if (set == c_Aligned) {
592 }
else if (set == c_Misaligned) {
594 }
else if (set == c_Aligned) {
601 boost::iostreams::filtering_istream ifs;
606 double back[np], fwrd[np], tension;
611 for (
int i = 0; i < np; ++i) {
614 for (
int i = 0; i < np; ++i) {
620 if (ifs.eof())
break;
628 for (
int i = 0; i < np; ++i) {
632 }
else if (set == c_Misaligned) {
635 }
else if (set == c_Aligned) {
647 }
else if (set == c_Misaligned) {
652 }
else if (set == c_Aligned) {
671 if (nRead !=
m_nSenseWires) B2FATAL(
"CDCGeometryPar::readWirePositionParams: #lines read-in (=" << nRead <<
672 ") is inconsistent with total #sense wires (=" <<
m_nSenseWires <<
") !");
675 boost::iostreams::close(ifs);
683 for (
unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
690 auto layerID =
WireID(iL, 511);
701 for (
unsigned iC = 0; iC <
m_nWires[iL]; ++iC) {
713 m_BWirPosAlign[iL][iC][0] = d_layerXbwd + cos(d_layerPhiBwd) * wireXbwd + sin(d_layerPhiBwd) * wireYbwd;
714 m_BWirPosAlign[iL][iC][1] = d_layerYbwd - sin(d_layerPhiBwd) * wireXbwd + cos(d_layerPhiBwd) * wireYbwd;
717 m_FWirPosAlign[iL][iC][0] = d_layerXfwd + cos(d_layerPhiFwd) * wireXfwd + sin(d_layerPhiFwd) * wireYfwd;
718 m_FWirPosAlign[iL][iC][1] = d_layerYfwd - sin(d_layerPhiFwd) * wireXfwd + cos(d_layerPhiFwd) * wireYfwd;
724 double back[np], fwrd[np];
726 for (
unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
732 for (
unsigned iC = 0; iC <
m_nWires[iL]; ++iC) {
743 for (
int i = 0; i < np; ++i) {
767 double back[np], fwrd[np];
769 for (
unsigned iL = 0; iL < c_maxNSenseLayers; ++iL) {
775 for (
unsigned iC = 0; iC <
m_nWires[iL]; ++iC) {
786 for (
int i = 0; i < np; ++i) {
821 fileName0 = gbxParams.
getString(
"xt4ReconFileName");
824 boost::iostreams::filtering_istream ifs;
848 unsigned short nAlphaBins = 0;
849 if (ifs >> nAlphaBins) {
850 if (nAlphaBins == 0 || nAlphaBins > c_maxNAlphaPoints) B2FATAL(
"Fail to read alpha bins !");
852 B2FATAL(
"Fail to read alpha bins !");
855 double alpha0, alpha1, alpha2;
856 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
857 ifs >> alpha0 >> alpha1 >> alpha2;
862 unsigned short nThetaBins = 0;
863 if (ifs >> nThetaBins) {
864 if (nThetaBins == 0 || nThetaBins > c_maxNThetaPoints) B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
866 B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
869 double theta0, theta1, theta2;
871 for (
unsigned short i = 0; i < nThetaBins; ++i) {
872 ifs >> theta0 >> theta1 >> theta2;
877 unsigned short iCL, iLR;
878 const unsigned short npx = c_nXTParams - 1;
880 double theta, alpha, dummy1;
885 if (np <= 0 || np > npx) B2FATAL(
"CDCGeometryPar: no. of xt-params. outside limits !");
887 const double epsi = 0.1;
895 ifs >> theta >> alpha >> dummy1 >> iLR;
896 for (
int i = 0; i < np; ++i) {
901 for (
unsigned short i = 0; i < nThetaBins; ++i) {
907 if (itheta < 0) B2FATAL(
"CDCGeometryPar: thetas in xt.dat are inconsistent !");
910 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
916 if (ialpha < 0) B2FATAL(
"CDCGeometryPar: alphas in xt.dat are inconsistent !");
918 for (
int i = 0; i < np; ++i) {
919 m_XT[iCL][iLR][ialpha][itheta][i] = xtc[i];
922 double boundT = xtc[6];
924 m_XT[iCL][iLR][ialpha][itheta][np] = ROOT::Math::Chebyshev5(boundT, xtc[0], xtc[1], xtc[2], xtc[3], xtc[4], xtc[5]);
926 m_XT[iCL][iLR][ialpha][itheta][np] =
937 boost::iostreams::close(ifs);
940 const double degrad = M_PI / 180.;
941 for (
unsigned i = 0; i < nAlphaBins; ++i) {
944 for (
unsigned i = 0; i < nThetaBins; ++i) {
967 fileName0 = gbxParams.
getString(
"sigma4ReconFileName");
975 unsigned short nAlphaBins = 0;
976 if (ifs >> nAlphaBins) {
977 if (nAlphaBins == 0 || nAlphaBins > c_maxNAlphaPoints) B2FATAL(
"Fail to read alpha bins !");
979 B2FATAL(
"Fail to read alpha bins !");
983 double alpha0, alpha1, alpha2;
984 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
985 ifs >> alpha0 >> alpha1 >> alpha2;
991 unsigned short nThetaBins = 0;
992 if (ifs >> nThetaBins) {
993 if (nThetaBins == 0 || nThetaBins > c_maxNThetaPoints) B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
995 B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
999 double theta0, theta1, theta2;
1001 for (
unsigned short i = 0; i < nThetaBins; ++i) {
1002 ifs >> theta0 >> theta1 >> theta2;
1007 unsigned short np = 0;
1008 unsigned short iCL, iLR;
1009 double sigma[c_nSigmaParams];
1010 double theta, alpha;
1016 if (np > c_nSigmaParams) B2FATAL(
"CDCGeometryPar: no. of sigma-params. outside limits !");
1020 const double epsi = 0.1;
1022 while (ifs >> iCL) {
1028 ifs >> theta >> alpha >> iLR;
1030 for (
int i = 0; i < np; ++i) {
1035 for (
unsigned short i = 0; i < nThetaBins; ++i) {
1041 if (itheta < 0) B2FATAL(
"CDCGeometryPar: thetas in sigma.dat are inconsistent !");
1044 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
1050 if (ialpha < 0) B2FATAL(
"CDCGeometryPar: alphas in sigma.dat are inconsistent !");
1052 for (
int i = 0; i < np; ++i) {
1053 m_Sigma[iCL][iLR][ialpha][itheta][i] = sigma[i];
1060 const double degrad = M_PI / 180.;
1061 for (
unsigned i = 0; i < nAlphaBins; ++i) {
1064 for (
unsigned i = 0; i < nThetaBins; ++i) {
1077 fileName0 = gbxParams.
getString(
"fudgeFactorFileName");
1079 B2WARNING(
"readFFactor is not ready! " << fileName0);
1089 fileName0 = gbxParams.
getString(
"propSpeed4ReconFileName");
1102 if (ifs.eof())
break;
1108 if (
m_debug) B2DEBUG(150, iL <<
" " << speed);
1111 if (nRead != c_maxNSenseLayers) B2FATAL(
"CDCGeometryPar::readPropSpeed: #lines read-in (=" << nRead <<
1112 ") is inconsistent with total #layers (=" << c_maxNSenseLayers <<
") !");
1157 fileName0 = gbxParams.
getString(
"t04ReconFileName");
1169 ifs >> iL >> iC >> t0;
1175 if (ifs.eof())
break;
1182 B2DEBUG(150, iL <<
" " << iC <<
" " << t0);
1186 if (nRead !=
m_nSenseWires) B2FATAL(
"CDCGeometryPar::readT0: #lines read-in (=" << nRead <<
1187 ") is inconsistent with total #sense wires (=" <<
m_nSenseWires <<
") !");
1238 fileName0 = gbxParams.
getString(
"tw4ReconFileName");
1245 unsigned short nPars(0);
1248 B2FATAL(
"CDCGeometryPar::readTW: invalid mode specified!");
1251 B2FATAL(
"CDCGeometryPar::readTW: invalid #params specified!");
1254 unsigned iBoard = 0;
1257 while (ifs >> iBoard) {
1258 for (
unsigned short i = 0; i < nPars; ++i) {
1264 if (nRead != c_nBoards) B2FATAL(
"CDCGeometryPar::readTW: #lines read-in (=" << nRead <<
") is inconsistent with #boards (=" <<
1282 unsigned short iSL, iL, iW, iB, iC;
1287 ifs >> iSL >> iL >> iW >> iB >> iC;
1288 if (ifs.eof())
break;
1293 m_wireToBoard.insert(pair<WireID, unsigned short>(wID, iB));
1296 if (nRead !=
m_nSenseWires) B2FATAL(
"CDCGeometryPar::readChMap: #lines read-in (=" << nRead <<
1297 ") is inconsistent with #sense-wires (="
1310 fileName0 = gbxParams.
getString(
"fudgeFactorFileName");
1315 std::string fileName1 =
"/data/cdc/" + fileName0;
1318 if (fileName ==
"") {
1322 if (fileName ==
"") {
1323 B2FATAL(
"CDC::openFile: " << fileName0 <<
" not exist!");
1326 B2DEBUG(29,
"CDC::openFile: open " << fileName);
1327 ifs.open(fileName.c_str());
1328 if (!ifs) B2FATAL(
"CDC::openFile: cannot open " << fileName <<
" !");
1331 unsigned short paramMode(0), nParams(0);
1332 ifs >> paramMode >> nParams;
1333 if (paramMode > 1) B2FATAL(
"Param mode > 1!");
1334 if (nParams > 7) B2FATAL(
"No. of params. > 7!");
1335 unsigned short groupId(0);
1337 B2DEBUG(29, paramMode <<
" " << nParams <<
" " << groupId);
1338 if (groupId > 0) B2FATAL(
"GgroupId > 0!");
1341 unsigned short cLMin[c_nSuperLayers], cLMax[c_nSuperLayers];
1344 for (
unsigned int sl = 1; sl < c_nSuperLayers; ++sl) {
1345 cLMin[sl] = cLMax[0] + 6 * sl - 5;
1346 cLMax[sl] = cLMax[0] + 6 * sl;
1349 unsigned short id = 0;
1351 unsigned short nRead = 0;
1353 for (
unsigned short i = 0; i < nParams; ++i) {
1355 for (
unsigned short cL = cLMin[
id]; cL <= cLMax[id]; ++cL) {
1356 for (
unsigned short cell = 0; cell <
m_nWires[cL]; ++cell) {
1363 if (nRead > c_nSuperLayers) B2FATAL(
"No. of read in lines > " << c_nSuperLayers <<
" !");
1373 for (
unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1374 for (
unsigned short iW = 0; iW < c_maxNDriftCells; ++iW) {
1379 for (
auto const& ent : (*m_t0FromDB)->getT0s()) {
1382 const unsigned short iW = wid.
getIWire();
1393 double oldMeanT0 = 0;
1394 unsigned short it1 = 0;
1395 for (
unsigned short it = 0; it < maxIt; ++it) {
1397 double effiSum = 0.;
1400 for (
unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1401 for (
unsigned short iW = 0; iW <
m_nWires[iCL]; ++iW) {
1402 if (
m_t0[iCL][iW] < minT0 ||
m_t0[iCL][iW] > maxT0)
continue;
1417 B2DEBUG(29, it <<
" " << effiSum <<
" " <<
m_meanT0 <<
" " << stdvT0);
1418 if (fabs(
m_meanT0 - oldMeanT0) < epsi)
break;
1423 B2FATAL(
"Wire efficiency sum <= 0!");
1426 if (it1 == maxIt - 1) B2WARNING(
"Max. iterations(=" << maxIt <<
") needed to calculate the mean t0. Strange.");
1446 for (
unsigned short iCL = 0; iCL < (*m_propSpeedFromDB)->getEntries(); ++iCL) {
1458 for (
unsigned short iBd = 0; iBd < (*m_timeWalkFromDB)->getEntries(); ++iBd) {
1459 int np = ((*m_timeWalkFromDB)->getTimeWalkParams(iBd)).size();
1460 for (
int i = 0; i < np; ++i) {
1461 m_timeWalkCoef[iBd][i] = ((*m_timeWalkFromDB)->getTimeWalkParams(iBd))[i];
1487 for (
unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1494 for (
unsigned short iLR = 0; iLR < 2; ++iLR) {
1497 const std::vector<float> params = (*m_xtRelFromDB)->getXtParams(iCL, iLR, iA, iT);
1498 unsigned short np = params.size();
1500 for (
unsigned short i = 0; i < np; ++i) {
1501 m_XT[iCL][iLR][iA][iT][i] = params[i];
1504 double boundT =
m_XT[iCL][iLR][iA][iT][6];
1506 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],
1507 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]);
1509 m_XT[iCL][iLR][iA][iT][np] =
1510 m_XT[iCL][iLR][iA][iT][0] + boundT
1511 * (
m_XT[iCL][iLR][iA][iT][1] + boundT
1512 * (
m_XT[iCL][iLR][iA][iT][2] + boundT
1513 * (
m_XT[iCL][iLR][iA][iT][3] + boundT
1514 * (
m_XT[iCL][iLR][iA][iT][4] + boundT
1515 * (
m_XT[iCL][iLR][iA][iT][5])))));
1550 for (
unsigned short iCL = 0; iCL < c_maxNSenseLayers; ++iCL) {
1551 for (
unsigned short iLR = 0; iLR < 2; ++iLR) {
1554 const std::vector<float> params = (*m_sResolFromDB)->getSigmaParams(iCL, iLR, iA, iT);
1555 unsigned short np = params.size();
1557 for (
unsigned short i = 0; i < np; ++i) {
1571 unsigned short groupId = (*m_fFactorFromDB)->getGroupID();
1572 unsigned short nEnt = (*m_fFactorFromDB)->getEntries();
1573 B2DEBUG(29,
"setFFactor called: groupId,nEnt= " << groupId <<
" " << nEnt);
1577 B2FATAL(
"CDCGeometryPar:: Invalid group-id " << groupId <<
" specified!");
1580 for (
unsigned short id = 0;
id < nEnt; ++id) {
1581 unsigned short np = ((*m_fFactorFromDB)->getFactors(
id)).size();
1582 if (np != 3) B2FATAL(
"CDCGeometryPar:: No. of fudge factors != 3!");
1583 for (
unsigned short i = 0; i < np; ++i) {
1600 const unsigned short isl = cm.getISuperLayer();
1602 const uint il = cm.getILayer();
1603 const int iw = cm.getIWire();
1604 const int iBd = cm.getBoardID();
1605 const WireID wID(isl, il, iw);
1606 m_wireToBoard.insert(pair<WireID, unsigned short>(wID, iBd));
1607 const int iCh = cm.getBoardChannel();
1616 unsigned short groupId = (*m_eDepToADCConversionsFromDB)->getGroupID();
1617 unsigned short nEnt = (*m_eDepToADCConversionsFromDB)->getEntries();
1619 if (nEnt > c_nSuperLayers) B2FATAL(
"CDCGeometryPar:: group-id " << groupId <<
" and #entries " << nEnt <<
" are inconsistent!");
1620 }
else if (groupId == 1) {
1621 if (nEnt > c_maxNSenseLayers) B2FATAL(
"CDCGeometryPar:: group-id " << groupId <<
" and #entries " << nEnt <<
" are inconsistent!");
1623 B2FATAL(
"CDCGeometryPar:: Invalid group-id " << groupId <<
" specified !");
1626 unsigned short cLMin[c_nSuperLayers], cLMax[c_nSuperLayers];
1629 for (
unsigned int sl = 1; sl < c_nSuperLayers; ++sl) {
1630 cLMin[sl] = cLMax[0] + 6 * sl - 5;
1631 cLMax[sl] = cLMax[0] + 6 * sl;
1634 for (
unsigned short id = 0;
id < nEnt; ++id) {
1635 unsigned short np = ((*m_eDepToADCConversionsFromDB)->getParams(
id)).size();
1636 if (np > 7) B2FATAL(
"CDCGeometryPar:: No. of edep-to-ADC conversion params. > 7");
1638 for (
unsigned short cL = cLMin[
id]; cL <= cLMax[id]; ++cL) {
1639 for (
unsigned short cell = 0; cell <
m_nWires[cL]; ++cell) {
1640 for (
unsigned short i = 0; i < np; ++i) {
1645 }
else if (groupId == 1) {
1646 for (
unsigned short cell = 0; cell <
m_nWires[id]; ++cell) {
1647 for (
unsigned short i = 0; i < np; ++i) {
1651 }
else if (groupId == 2) {
1653 B2FATAL(
"CDCGeometryPar::setEDepToADCConversions(): groupId=2 not ready!");
1672 const double cth = fabs(costh) + dlt;
1673 const double iGen = edep / dx;
1674 const double tmp = cth - gam * iGen;
1675 const double disc = tmp * tmp + 4.*alf * cth * iGen;
1679 iMea = cth * iGen / tmp;
1680 }
else if (disc >= 0.) {
1681 iMea = (-tmp +
sqrt(disc)) / (2.*alf);
1689 double convF = mainF;
1691 convF = mainF * std::min(iMea / iGen, 1.);
1700 convF *= 1. + a * (costh - b);
1720 if (set == c_Misaligned) {
1724 }
else if (set == c_Base) {
1746 if (set == c_Misaligned) {
1749 }
else if (set == c_Base) {
1768 if (set == c_Misaligned) {
1772 }
else if (set == c_Base) {
1794 if (set == c_Misaligned) {
1797 }
else if (set == c_Base) {
1808 if (set == c_Misaligned) {
1810 }
else if (set == c_Aligned) {
1818 static double IRWL[c_maxNSenseLayers] = {0};
1830 static double ORWL[c_maxNSenseLayers] = {0};
1846 const unsigned nWires =
m_nWires[layerId];
1850 const double phiSize = 2 * M_PI / double(nWires);
1867 for (
unsigned i = 0; i < 1; ++i) {
1868 const double phiF = phiSize * (double(i) +
offset)
1875 const double beta = (position.
Z() - b.Z()) / u.Z();
1877 double dPhi = std::atan2(position.
Y(), position.
X())
1878 - std::atan2(p.Y(), p.X())
1880 while (dPhi < 0) dPhi += (2. * M_PI);
1881 j = int(dPhi / phiSize);
1882 while (j >= nWires) j -= nWires;
1891 std::ofstream ofs(of.c_str(), std::ios::out);
1893 B2ERROR(
"CDCGeometryPar::read !!! can not open file : "
1896 ofs <<
"<?xml version=\"1.0\" encoding=\"UTF-8\"?>"
1898 <<
"<Subdetector type=\"CDC\">"
1900 <<
" <Name>CDC BelleII </Name>"
1902 <<
" <Description>CDC geometry parameters</Description>"
1904 <<
" <Version>0</Version>"
1906 <<
" <GeoCreator>CDCBelleII</GeoCreator>"
1910 <<
" <Rotation desc=\"Rotation of the whole cdc detector (should be the same as beampipe)\" unit=\"mrad\">0.0</Rotation>"
1912 <<
" <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>"
1914 <<
" <Material>CDCGas</Material>"
1918 ofs <<
" <SLayers>" << endl;
1921 ofs <<
" <SLayer id=\"" << i <<
"\">" << endl;
1922 ofs <<
" <Radius desc=\"Radius of wires in this layer\" unit=\"mm\">" <<
senseWireR(i) <<
"</Radius>" << endl;
1923 ofs <<
" <BackwardZ desc=\"z position of this wire layer at backward endplate\" unit=\"mm\">" <<
senseWireBZ(
1924 i) <<
"</BackwardZ>" << endl;
1925 ofs <<
" <ForwardZ desc=\"z position of this wire layer at forward endplate\" unit=\"mm\">" <<
senseWireFZ(
1926 i) <<
"</ForwardZ>" << endl;
1929 ofs <<
" <NHoles desc=\"the number of holes in this layer, 2*(cell number)\">" <<
nWiresInLayer(
1930 i) * 2 <<
"</NHoles>" << endl;
1931 ofs <<
" <NShift desc=\"the shifted hole number of each wire in this layer\">" <<
nShifts(i) <<
"</NShift>" << endl;
1932 ofs <<
" <Offset desc=\"wire offset in phi direction at endplate\">" <<
m_offSet[i] <<
"</Offset>" << endl;
1933 ofs <<
" </SLayer>" << endl;
1936 ofs <<
" </SLayers>" << endl;
1937 ofs <<
" <FLayers>" << endl;
1940 ofs <<
" <FLayer id=\"" << i <<
"\">" << endl;
1941 ofs <<
" <Radius desc=\"Radius of field wires in this layer\" unit=\"mm\">" <<
fieldWireR(i) <<
"</Radius>" << endl;
1942 ofs <<
" <BackwardZ desc=\"z position of this field wire layer at backward endplate\" unit=\"mm\">" <<
fieldWireBZ(
1943 i) <<
"</BackwardZ>" << endl;
1944 ofs <<
" <ForwardZ desc=\"z position of this field wire layer at forward endplate\" unit=\"mm\">" <<
fieldWireFZ(
1945 i) <<
"</ForwardZ>" << endl;
1946 ofs <<
" </FLayer>" << endl;
1949 ofs <<
" </FLayers>" << endl;
1951 ofs <<
" <InnerWall name=\"InnerWall\">" << endl;
1952 ofs <<
" <InnerR desc=\"Inner radius\" unit=\"mm\">" <<
innerRadiusInnerWall() <<
"</InnerR>" << endl;
1953 ofs <<
" <OuterR desc=\"Outer radius\" unit=\"mm\">" <<
outerRadiusInnerWall() <<
"</OuterR>" << endl;
1954 ofs <<
" <BackwardZ desc=\"z position at backward endplate\" unit=\"mm\">" <<
m_zWall[0][0] <<
"</BackwardZ>" << endl;
1955 ofs <<
" <ForwardZ desc=\"z position at forward endplate\" unit=\"mm\">" <<
m_zWall[0][1] <<
"</ForwardZ>" << endl;
1956 ofs <<
" </InnerWall>" << endl;
1958 ofs <<
" <OuterWall name=\"OuterWall\">" << endl;
1959 ofs <<
" <InnerR desc=\"Inner radius\" unit=\"mm\">" <<
innerRadiusOuterWall() <<
"</InnerR>" << endl;
1960 ofs <<
" <OuterR desc=\"Outer radius\" unit=\"mm\">" <<
outerRadiusOuterWall() <<
"</OuterR>" << endl;
1961 ofs <<
" <BackwardZ desc=\"z position at backward endplate\" unit=\"mm\">" <<
m_zWall[2][0] <<
"</BackwardZ>" << endl;
1962 ofs <<
" <ForwardZ desc=\"z position at forward endplate\" unit=\"mm\">" <<
m_zWall[2][1] <<
"</ForwardZ>" << endl;
1963 ofs <<
" </OuterWall>" << endl;
1965 ofs <<
" </Content>" << endl
1966 <<
"</Subdetector>" << endl;
1970 double& Yb_sag,
double& Yf_sag)
const
2001 if (set == c_Aligned) {
2015 }
else if (set == c_Misaligned) {
2029 }
else if (set == c_Base) {
2044 B2FATAL(
"CDCGeometryPar::getWireSagEffect: called with an invalid set: " <<
" " << set);
2047 const double dx = Xf - Xb;
2048 const double dy = Yf - Yb;
2049 const double dz = Zf - Zb;
2051 const double Zfp =
sqrt(dz * dz + dx * dx);
2052 const double Zp = (Z - Zb) * Zfp / dz;
2054 const double Y_sag = (Coef * (Zp - Zfp) + dy / Zfp) * Zp + Yb;
2055 const double dydz = (Coef * (2.*Zp - Zfp) * Zfp + dy) / dz;
2057 Yb_sag = Y_sag + dydz * (Zb - Z);
2058 Yf_sag = Y_sag + dydz * (Zf - Z);
2064 const unsigned L = layerID;
2065 const unsigned C = cellID;
2069 const double phiSize = 2 * M_PI / double(
m_nWires[L]);
2071 const double phiF = phiSize * (double(C) +
offset)
2084 for (
int i = 0; i < 3; ++i) {
2102 const unsigned L = layerID;
2103 const unsigned C = cellID;
2105 static bool first =
true;
2106 static ofstream ofs;
2109 ofs.open(
"alignment.dat");
2112 ofs << L <<
" " << C;
2114 ofs << setiosflags(ios::showpoint | ios::uppercase);
2116 for (
int i = 0; i < 3; ++i) ofs <<
" " << setw(15) << setprecision(8) <<
m_BWirPos[L][C][i];
2118 for (
int i = 0; i < 3; ++i) ofs <<
" " << setw(15) << setprecision(8) <<
m_FWirPos[L][C][i];
2119 ofs << setiosflags(ios::fixed);
2126 const double theta)
const
2128 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2136 double delta = time - minTime;
2142 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2145 unsigned short ial[2] = {0};
2146 unsigned short ilr[2] = {lro, lro};
2149 unsigned short ith[2] = {0};
2152 unsigned short jal(0), jlr(0), jth(0);
2156 double timep = delta < 0. ? minTime - delta : time;
2159 for (
unsigned k = 0; k < 4; ++k) {
2164 w = (1. - wal) * (1. - wth);
2165 }
else if (k == 1) {
2169 w = (1. - wal) * wth;
2170 }
else if (k == 2) {
2174 w = wal * (1. - wth);
2175 }
else if (k == 3) {
2182 double boundary =
m_XT[iCLayer][jlr][jal][jth][6];
2184 if (timep < boundary) {
2186 const double& c1 =
m_XT[iCLayer][jlr][jal][jth][1];
2187 const double& c2 =
m_XT[iCLayer][jlr][jal][jth][2];
2188 const double& c3 =
m_XT[iCLayer][jlr][jal][jth][3];
2189 const double& c4 =
m_XT[iCLayer][jlr][jal][jth][4];
2190 const double& c5 =
m_XT[iCLayer][jlr][jal][jth][5];
2191 dDdt += w * ROOT::Math::Chebyshev4(timep, c1 + 3.*c3 + 5.*c5, 4.*c2 + 8.*c4, 6.*c3 + 10.*c5, 8.*c4, 10.*c5);
2193 dDdt += w * (
m_XT[iCLayer][jlr][jal][jth][1] + timep
2194 * (2.*
m_XT[iCLayer][jlr][jal][jth][2] + timep
2195 * (3.*
m_XT[iCLayer][jlr][jal][jth][3] + timep
2196 * (4.*
m_XT[iCLayer][jlr][jal][jth][4] + timep
2197 * (5.*
m_XT[iCLayer][jlr][jal][jth][5])))));
2200 dDdt += w *
m_XT[iCLayer][jlr][jal][jth][7];
2214 const double theta)
const
2216 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2228 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2231 unsigned short ial[2] = {0};
2232 unsigned short ilr[2] = {lro, lro};
2235 unsigned short ith[2] = {0};
2238 unsigned short jal(0), jlr(0), jth(0);
2242 double timep = time;
2246 for (
unsigned k = 0; k < 4; ++k) {
2251 w = (1. - wal) * (1. - wth);
2252 }
else if (k == 1) {
2256 w = (1. - wal) * wth;
2257 }
else if (k == 2) {
2261 w = wal * (1. - wth);
2262 }
else if (k == 3) {
2269 double boundary =
m_XT[iCLayer][jlr][jal][jth][6];
2271 if (timep < boundary) {
2273 dist += w * ROOT::Math::Chebyshev5(timep,
m_XT[iCLayer][jlr][jal][jth][0],
m_XT[iCLayer][jlr][jal][jth][1],
2274 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]);
2276 dist += w * (
m_XT[iCLayer][jlr][jal][jth][0] + timep
2277 * (
m_XT[iCLayer][jlr][jal][jth][1] + timep
2278 * (
m_XT[iCLayer][jlr][jal][jth][2] + timep
2279 * (
m_XT[iCLayer][jlr][jal][jth][3] + timep
2280 * (
m_XT[iCLayer][jlr][jal][jth][4] + timep
2281 * (
m_XT[iCLayer][jlr][jal][jth][5]))))));
2284 dist += w * (
m_XT[iCLayer][jlr][jal][jth][7] * (timep - boundary) +
m_XT[iCLayer][jlr][jal][jth][8]);
2297 const bool calculateMinTime,
2298 const double inputMinTime)
const
2300 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2307 double minTime = calculateMinTime ?
getMinDriftTime(iCLayer, lr, alpha, theta) : inputMinTime;
2308 double delta = time - minTime;
2316 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2319 unsigned short ial[2] = {0};
2320 unsigned short ilr[2] = {lro, lro};
2323 unsigned short ith[2] = {0};
2326 unsigned short jal(0), jlr(0), jth(0);
2330 double timep = delta < 0. ? minTime - delta : time;
2335 for (
unsigned k = 0; k < 4; ++k) {
2340 w = (1. - wal) * (1. - wth);
2341 }
else if (k == 1) {
2345 w = (1. - wal) * wth;
2346 }
else if (k == 2) {
2350 w = wal * (1. - wth);
2351 }
else if (k == 3) {
2371 double boundary =
m_XT[iCLayer][jlr][jal][jth][6];
2373 if (timep < boundary) {
2375 dist += w * ROOT::Math::Chebyshev5(timep,
m_XT[iCLayer][jlr][jal][jth][0],
m_XT[iCLayer][jlr][jal][jth][1],
2376 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]);
2378 dist += w * (
m_XT[iCLayer][jlr][jal][jth][0] + timep
2379 * (
m_XT[iCLayer][jlr][jal][jth][1] + timep
2380 * (
m_XT[iCLayer][jlr][jal][jth][2] + timep
2381 * (
m_XT[iCLayer][jlr][jal][jth][3] + timep
2382 * (
m_XT[iCLayer][jlr][jal][jth][4] + timep
2383 * (
m_XT[iCLayer][jlr][jal][jth][5]))))));
2386 dist += w * (
m_XT[iCLayer][jlr][jal][jth][7] * (timep - boundary) +
m_XT[iCLayer][jlr][jal][jth][8]);
2393 if (delta < 0.) dist *= -1.;
2399 const double theta)
const
2401 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2405 double minTime = 0.;
2411 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2414 unsigned short ial[2] = {0};
2415 unsigned short ilr[2] = {lro, lro};
2418 unsigned short ith[2] = {0};
2421 unsigned short jal(0), jlr(0), jth(0);
2424 double c[6] = {0.}, a[6] = {0.};
2425 for (
unsigned k = 0; k < 4; ++k) {
2430 w = (1. - wal) * (1. - wth);
2431 }
else if (k == 1) {
2435 w = (1. - wal) * wth;
2436 }
else if (k == 2) {
2440 w = wal * (1. - wth);
2441 }
else if (k == 3) {
2448 for (
int i = 0; i < 5; ++i) {
2449 c[i] += w *
m_XT[iCLayer][jlr][jal][jth][i];
2454 a[0] = c[0] - c[2] + c[4];
2455 a[1] = c[1] - 3.*c[3] + 5.*c[5];
2456 a[2] = 2.*c[2] - 8.*c[4];
2457 a[3] = 4.*c[3] - 20.*c[5];
2461 for (
int i = 0; i < 5; ++i) a[i] = c[i];
2468 const double det = a[1] * a[1] - 4.*a[2] * a[0];
2471 minTime = (-a[1] +
sqrt(det)) / (2.*a[2]);
2475 minTime = -a[1] / (2.*a[2]);
2478 }
else if (a[1] != 0.) {
2479 minTime = -a[0] / a[1];
2481 B2WARNING(
"CDCGeometryPar::getMinDriftTime: minDriftTime not determined; assume zero.\n" <<
"layer(#0-55),lr,alpha(rad),theta= " <<
2482 iCLayer <<
" " << lr <<
" " << alpha <<
" " << theta);
2491 const double epsi4x = 5.e-6;
2493 const unsigned short maxIter = 8;
2494 const double maxDt = 20.;
2495 unsigned short nIter = 0;
2496 double minXsq = 1.e10;
2497 double minMinTime = minTime;
2500 for (nIter = 0; nIter <= maxIter; ++nIter) {
2503 double x = a[0] + t * (a[1] + t * (a[2] + t * (a[3] + t * (a[4] + t * a[5]))));
2509 double xp = a[1] + t * (2 * a[2] + t * (3 * a[3] + t * (4 * a[4] + t * 5 * a[5])));
2510 double xpp = 2 * a[2] + t * (6 * a[3] + t * (12 * a[4] + t * 20 * a[5]));
2511 double den = xp * xp + x * xpp;
2518 edm = fabs(x * xp) /
sqrt(den);
2519 if (edm < epsi4x)
break;
2526 dt = std::min(dt, maxDt);
2528 dt = std::max(dt, -maxDt);
2531 B2WARNING(
"CDCGeometryPar::getMinDriftTime: den = 0\n" <<
"layer(#0-55),lr,alpha(rad),theta= " <<
2534 " " << alpha <<
" " << theta);
2540 if (nIter == (maxIter + 1)) minTime = minMinTime;
2593 const double theta)
const
2595 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2601 const double eps = 2.5e-1;
2602 const double maxTrials = 100;
2610 double maxTime = 2000.;
2616 double t0 = minTime;
2618 const bool calMinTime =
false;
2623 double t1 = maxTime;
2625 while (((t1 - t0) > eps) && (i < maxTrials)) {
2626 time = 0.5 * (t0 + t1);
2627 double d1 =
getDriftLength(time, iCLayer, lr, alpha, theta, calMinTime, minTime) - dist;
2637 if (i >= maxTrials - 1 || time > maxTime) {
2638 B2WARNING(
"CDCGeometryPar::getDriftTime " << dist <<
" " << iCLayer <<
" " << alpha <<
" " << lr <<
" " << t0 <<
" " << t1 <<
" " <<
2649 const double theta)
const
2651 if (iCLayer < m_firstLayerOffset || iCLayer >= c_maxNSenseLayers) {
2658 const double driftL = fabs(DriftL0);
2664 B2FATAL(
"linearInterpolationOfSgm = false is not allowed now !");
2668 unsigned short ial[2] = {0};
2669 unsigned short ilr[2] = {lro, lro};
2672 unsigned short ith[2] = {0};
2676 unsigned short jal(0), jlr(0), jth(0);
2678 for (
unsigned k = 0; k < 4; ++k) {
2683 w = (1. - wal) * (1. - wth);
2684 }
else if (k == 1) {
2688 w = (1. - wal) * wth;
2689 }
else if (k == 2) {
2693 w = wal * (1. - wth);
2694 }
else if (k == 3) {
2712 const double& P0 =
m_Sigma[iCLayer][jlr][jal][jth][0];
2713 const double& P1 =
m_Sigma[iCLayer][jlr][jal][jth][1];
2714 const double& P2 =
m_Sigma[iCLayer][jlr][jal][jth][2];
2715 const double& P3 =
m_Sigma[iCLayer][jlr][jal][jth][3];
2716 const double& P4 =
m_Sigma[iCLayer][jlr][jal][jth][4];
2717 const double& P5 =
m_Sigma[iCLayer][jlr][jal][jth][5];
2718 const double& P6 =
m_Sigma[iCLayer][jlr][jal][jth][6];
2720#if defined(CDC_DEBUG)
2721 cout <<
"driftL= " << driftL << endl;
2722 cout <<
"iCLayer= " << iCLayer <<
" " << jlr <<
" " << jal <<
" " << jth << endl;
2723 cout <<
"P0= " << P0 << endl;
2724 cout <<
"P1= " << P1 << endl;
2725 cout <<
"P2= " << P2 << endl;
2726 cout <<
"P3= " << P3 << endl;
2727 cout <<
"P4= " << P4 << endl;
2728 cout <<
"P5= " << P5 << endl;
2729 cout <<
"P6= " << P6 << endl;
2734 sigma += w *
sqrt(P0 / (driftL * driftL + P1) + P2 * driftL + P3 +
2735 P4 * exp(P5 * (driftL - P6) * (driftL - P6)));
2737 double forthTermAtP7 = P4 * exp(P5 * (P7 - P6) * (P7 - P6));
2738 const double& P8 =
m_Sigma[iCLayer][jlr][jal][jth][8];
2740 double sigmaAtP7 =
sqrt(P0 / (P7 * P7 + P1) + P2 * P7 + P3 + forthTermAtP7);
2741 sigma += w * (P8 * (driftL - P7) + sigmaAtP7);
2743 double onePls4AtP7 =
sqrt(P0 / (P7 * P7 + P1) + forthTermAtP7);
2744 const double onePls4 = P8 * (driftL - P7) + onePls4AtP7;
2745 sigma += w *
sqrt(P2 * driftL + P3 + onePls4 * onePls4);
2747 forthTermAtP7 =
sqrt(forthTermAtP7);
2748 const double forthTerm = P8 * (driftL - P7) + forthTermAtP7;
2749 sigma += w *
sqrt(P0 / (driftL * driftL + P1) + P2 * driftL + P3 +
2750 forthTerm * forthTerm);
2763 unsigned short lr = 0;
2764 double wCrossT = (posOnWire.
Cross(posOnTrack)).Z();
2768 }
else if (wCrossT > 0.) {
2771 if ((posOnTrack - posOnWire).Perp() != 0.) {
2772 double wCrossP = (posOnWire.
Cross(momentum)).Z();
2774 if (posOnTrack.
Perp() > posOnWire.
Perp()) {
2779 }
else if (wCrossP < 0.) {
2780 if (posOnTrack.
Perp() < posOnWire.
Perp()) {
2798 const double distanceCrossP = ((posOnWire - posOnTrack).Cross(momentum)).Z();
2799 unsigned short int lr = (distanceCrossP > 0.) ? 1 : 0;
2806 const double wx = posOnWire.
X();
2807 const double wy = posOnWire.
Y();
2808 const double px = momentum.X();
2809 const double py = momentum.Y();
2811 const double cross = wx * py - wy * px;
2812 const double dot = wx * px + wy * py;
2814 return atan2(cross,
dot);
2819 return atan2(momentum.Perp(), momentum.Z());
2825 unsigned short lro = (fabs(alpha) <= 0.5 * M_PI) ? lr : abs(lr - 1);
2833 double alphao = alpha;
2834 if (alpha > 0.5 * M_PI) {
2836 }
else if (alpha < -0.5 * M_PI) {
2846 double thetao = fabs(alpha) > 0.5 * M_PI ? M_PI - theta : theta;
2852 unsigned short lrs[2])
const
2861 lrs[0] = abs(lrs[0] - 1);
2868 lrs[1] = abs(lrs[1] - 1);
2875 points[0] = points[1] - 1;
2882 unsigned short lrs[2])
const
2891 lrs[0] = abs(lrs[0] - 1);
2898 lrs[1] = abs(lrs[1] - 1);
2905 points[0] = points[1] - 1;
2931 points[0] = points[1] - 1;
2938 unsigned short points[2])
const
2954 points[0] = points[1] - 1;
2986 const unsigned short nLayers[c_nSuperLayers] = {8, 6, 6, 6, 6, 6, 6, 6, 6};
2988 for (
unsigned short SLayer = 0; SLayer < c_nSuperLayers; ++SLayer) {
2989 unsigned short firstCLayer = 0;
2990 for (
unsigned short i = 0; i < SLayer; ++i) {
2991 firstCLayer += nLayers[i];
2996 for (
unsigned short Layer = 0; Layer < nLayers[SLayer]; ++Layer) {
2997 unsigned short CLayer = firstCLayer + Layer;
2999 if (CLayer == firstCLayer) {
3002 }
else if (CLayer == firstCLayer + 1) {
3008 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.