11 #include <framework/gearbox/Gearbox.h>
12 #include <framework/gearbox/GearDir.h>
13 #include <framework/logging/Logger.h>
14 #include <framework/utilities/FileSystem.h>
16 #include <cdc/geometry/CDCGeometryPar.h>
17 #include <cdc/geometry/CDCGeoControlPar.h>
18 #include <cdc/simulation/CDCSimControlPar.h>
19 #include <cdc/utilities/OpenFile.h>
24 #include <boost/format.hpp>
28 #include <boost/iostreams/filtering_stream.hpp>
32 #include <Math/ChebyshevPol.h>
35 using namespace boost;
43 if (!m_B4CDCGeometryParDB) m_B4CDCGeometryParDB =
new CDCGeometryPar(geom);
44 return *m_B4CDCGeometryParDB;
54 if ((*m_t0FromDB).isValid()) {
55 (*m_t0FromDB).
addCallback(
this, &CDCGeometryPar::setT0);
61 if ((*m_badWireFromDB).isValid()) {
62 (*m_badWireFromDB).
addCallback(
this, &CDCGeometryPar::setBadWire);
68 if ((*m_propSpeedFromDB).isValid()) {
69 (*m_propSpeedFromDB).
addCallback(
this, &CDCGeometryPar::setPropSpeed);
75 if ((*m_timeWalkFromDB).isValid()) {
76 (*m_timeWalkFromDB).
addCallback(
this, &CDCGeometryPar::setTW);
82 if ((*m_xtRelFromDB).isValid()) {
83 (*m_xtRelFromDB).
addCallback(
this, &CDCGeometryPar::setXtRel);
89 if ((*m_sResolFromDB).isValid()) {
90 (*m_sResolFromDB).
addCallback(
this, &CDCGeometryPar::setSResol);
96 if ((*m_fFactorFromDB).isValid()) {
97 (*m_fFactorFromDB).
addCallback(
this, &CDCGeometryPar::setFFactor);
103 if ((*m_chMapFromDB).isValid()) {
104 (*m_chMapFromDB).
addCallback(
this, &CDCGeometryPar::setChMap);
110 if ((*m_displacementFromDB).isValid()) {
111 (*m_displacementFromDB).
addCallback(
this, &CDCGeometryPar::setDisplacement);
117 if ((*m_alignmentFromDB).isValid()) {
118 (*m_alignmentFromDB).
addCallback(
this, &CDCGeometryPar::setWirPosAlignParams);
125 if ((*m_misalignmentFromDB).isValid()) {
126 (*m_misalignmentFromDB).
addCallback(
this, &CDCGeometryPar::setWirPosMisalignParams);
133 if ((*m_eDepToADCConversionsFromDB).isValid()) {
134 (*m_eDepToADCConversionsFromDB).
addCallback(
this, &CDCGeometryPar::setEDepToADCConversions);
146 B2WARNING(
"CDCGeometryPar: Strange that readFromDB is not called! Please make sure that CDC is included in Geometry.");
150 CDCGeometryPar::~CDCGeometryPar()
168 void CDCGeometryPar::clear()
170 m_version =
"unknown";
173 m_senseWireDiameter = 0.0;
174 m_senseWireTension = 0.0;
175 m_senseWireDensity = 0.0;
176 m_fieldWireDiameter = 0.0;
179 m_clockFreq4TDC = 0.0;
181 m_nominalDriftV = 0.0;
182 m_nominalPropSpeed = 0.0;
183 m_nominalSpaceResol = 0.0;
185 for (
unsigned i = 0; i < 4; ++i) {
187 for (
unsigned j = 0; j < 2; ++j)
190 for (
unsigned i = 0; i < MAX_N_SLAYERS; ++i) {
192 m_zSForwardLayer[i] = 0;
193 m_zSBackwardLayer[i] = 0;
199 for (
unsigned i = 0; i < MAX_N_FLAYERS; ++i) {
201 m_zFForwardLayer[i] = 0;
202 m_zFBackwardLayer[i] = 0;
205 for (
unsigned L = 0; L < MAX_N_SLAYERS; ++L) {
206 for (
unsigned C = 0; C < MAX_N_SCELLS; ++C) {
207 for (
unsigned i = 0; i < 3; ++i) {
208 m_FWirPos [L][C][i] = 0.;
209 m_BWirPos [L][C][i] = 0.;
210 m_FWirPosMisalign[L][C][i] = 0.;
211 m_BWirPosMisalign[L][C][i] = 0.;
212 m_FWirPosAlign [L][C][i] = 0.;
213 m_BWirPosAlign [L][C][i] = 0.;
215 m_WireSagCoef [L][C] = 0.;
216 m_WireSagCoefMisalign[L][C] = 0.;
217 m_WireSagCoefAlign [L][C] = 0.;
225 m_globalPhiRotation = geom.getGlobalPhiRotation();
229 m_rWall[0] = geom.getInnerWall(2).getRmin();
230 m_zWall[0][0] = geom.getInnerWall(0).getZbwd();
231 m_zWall[0][1] = geom.getInnerWall(0).getZfwd();
233 m_rWall[1] = geom.getInnerWall(0).getRmax();
234 m_zWall[1][0] = geom.getInnerWall(0).getZbwd();
235 m_zWall[1][1] = geom.getInnerWall(0).getZbwd();
238 m_rWall[2] = geom.getOuterWall(0).getRmin();
239 m_zWall[2][0] = geom.getOuterWall(0).getZbwd();
240 m_zWall[2][1] = geom.getOuterWall(0).getZfwd();
242 m_rWall[3] = geom.getOuterWall(1).getRmax();
243 m_zWall[3][0] = geom.getOuterWall(0).getZbwd();
244 m_zWall[3][1] = geom.getOuterWall(0).getZfwd();
247 m_debug = CDCGeoControlPar::getInstance().getDebug();
248 m_nSLayer = geom.getNSenseLayers();
250 m_materialDefinitionMode = CDCGeoControlPar::getInstance().getMaterialDefinitionMode();
252 if (m_materialDefinitionMode == 0) {
253 B2DEBUG(100,
"CDCGeometryPar: Define a mixture of gases and wires in the tracking volume.");
254 }
else if (m_materialDefinitionMode == 2) {
256 B2FATAL(
"CDCGeometryPar: Materialdefinition=2 is disabled for now.");
258 B2FATAL(
"CDCGeometryPar: Materialdefinition mode you specify is invalid.");
262 m_senseWireZposMode = CDCGeoControlPar::getInstance().getSenseWireZposMode();
264 B2DEBUG(100,
"CDCGeometryPar: Sense wire z mode:" << m_senseWireZposMode);
269 GearDir content =
GearDir(
"/Detector/DetectorComponent[@name=\"CDC\"]/Content/");
277 for (
const auto& sense : geom.getSenseLayers()) {
278 int layerId = sense.getId();
279 m_rSLayer[layerId] = sense.getR();
280 m_zSBackwardLayer[layerId] = sense.getZbwd();
281 m_zSForwardLayer[layerId] = sense.getZfwd();
282 m_nWires[layerId] = sense.getNWires();
283 m_nShifts[layerId] = sense.getNShifts();
284 m_offSet[layerId] = sense.getOffset();
285 m_cellSize[layerId] = 2 * M_PI * m_rSLayer[layerId] / (double) m_nWires[layerId];
286 m_dzSBackwardLayer[layerId] = sense.getDZbwd();
287 m_dzSForwardLayer[layerId] = sense.getDZfwd();
290 if (m_senseWireZposMode == 0) {
291 }
else if (m_senseWireZposMode == 1) {
294 m_zSBackwardLayer[layerId] += m_dzSBackwardLayer[layerId];
295 m_zSForwardLayer [layerId] -= m_dzSForwardLayer [layerId];
297 B2FATAL(
"CDCGeometryPar: invalid wire z definition mode specified");
301 const int nWires = m_nWires[layerId];
302 for (
int iCell = 0; iCell < nWires; ++iCell) {
303 setDesignWirParam(layerId, iCell);
309 for (
const auto& field : geom.getFieldLayers()) {
310 int layerId = field.getId();
311 m_rFLayer[layerId] = field.getR();
312 m_zFBackwardLayer[layerId] = field.getZbwd();
313 m_zFForwardLayer[layerId] = field.getZfwd();
317 m_senseWireDiameter = geom.getSenseDiameter();
320 m_senseWireTension = geom.getSenseTension();
323 m_senseWireDensity = 19.3;
326 m_fieldWireDiameter = geom.getFieldDiameter();
329 m_clockFreq4TDC = 1.017774;
330 double tmp = geom.getClockFrequency();
332 if (tmp != m_clockFreq4TDC) {
333 B2WARNING(
"CDCGeometryPar: The default clock freq. for TDC (" << m_clockFreq4TDC <<
" GHz) is replaced with " << tmp <<
" (GHz).");
334 m_clockFreq4TDC = tmp;
336 B2DEBUG(100,
"CDCGeometryPar: Clock freq. for TDC= " << m_clockFreq4TDC <<
" (GHz).");
337 m_tdcBinWidth = 1. / m_clockFreq4TDC;
338 B2DEBUG(100,
"CDCGeometryPar: TDC bin width= " << m_tdcBinWidth <<
" (ns).");
340 m_nominalDriftV = 4.e-3;
341 m_nominalDriftVInv = 1. / m_nominalDriftV;
342 m_nominalPropSpeed = 27.25;
344 m_nominalSpaceResol = geom.getNominalSpaceResolution();
350 m_displacement = CDCGeoControlPar::getInstance().getDisplacement();
351 B2DEBUG(100,
"CDCGeometryPar: Load displacement params. (=1); not load (=0):" << m_displacement);
352 if (m_displacement) {
354 B2DEBUG(100,
"CDCGeometryPar: Read displacement from DB");
357 readWirePositionParams(c_Base, &geom);
362 m_alignment = CDCGeoControlPar::getInstance().getAlignment();
363 B2DEBUG(100,
"CDCGeometryPar: Load alignment params. (=1); not load (=0):" <<
367 B2DEBUG(100,
"CDCGeometryPar: Read alignment from DB");
368 setWirPosAlignParams();
370 readWirePositionParams(c_Aligned, &geom);
375 m_misalignment = CDCGeoControlPar::getInstance().getMisalignment();
376 B2DEBUG(100,
"CDCGeometryPar: Load misalignment params. (=1); not load (=0):" <<
378 if (m_misalignment) {
380 B2DEBUG(100,
"CDCGeometryPar: Read misalignment from DB");
381 setWirPosMisalignParams();
383 readWirePositionParams(c_Misaligned, &geom);
388 m_thresholdEnergyDeposit = CDCSimControlPar::getInstance().getThresholdEnergyDeposit();
389 m_minTrackLength = CDCSimControlPar::getInstance().getMinTrackLength();
390 m_wireSag = CDCSimControlPar::getInstance().getWireSag();
391 m_modLeftRightFlag = CDCSimControlPar::getInstance().getModLeftRightFlag();
392 if (m_modLeftRightFlag) {
393 B2FATAL(
"ModifiedLeftRightFlag = true is disabled for now; need to update a G4-related code in framework...");
397 m_sigmaFileFormat = 1;
402 B2DEBUG(100,
"CDCGeometryPar: Read xt from DB");
409 B2DEBUG(100,
"CDCGeometryPar: Read sigma from DB");
412 readSigma(gbxParams);
416 B2DEBUG(100,
"CDCGeometryPar: Read fudge factors from DB");
419 readFFactor(gbxParams);
423 B2DEBUG(100,
"CDCGeometryPar: Read prop-speed from DB");
426 readPropSpeed(gbxParams);
430 B2DEBUG(100,
"CDCGeometryPar: Read t0 from DB");
437 B2DEBUG(100,
"CDCGeometryPar: Read badwire from DB");
441 B2FATAL(
"Text file input mode for bdwires is disabled now!");
445 B2DEBUG(100,
"CDCGeometryPar: Read ch-map from DB");
452 B2DEBUG(100,
"CDCGeometryPar: Read time-walk from DB");
457 B2DEBUG(100,
"CDCGeometryPar: Time-walk param. mode= " << m_twParamMode);
460 B2DEBUG(29,
"CDCGeometryPar: Read EDepToADC from DB");
461 if ((*m_eDepToADCConversionsFromDB).isValid()) {
462 setEDepToADCConversions();
465 readEDepToADC(gbxParams);
473 readXT(gbxParams, 1);
474 readSigma(gbxParams, 1);
475 readPropSpeed(gbxParams, 1);
476 readT0(gbxParams, 1);
477 readTW(gbxParams, 1);
481 setShiftInSuperLayer();
533 std::string fileName0;
538 }
else if (set == c_Misaligned) {
540 }
else if (set == c_Aligned) {
546 }
else if (set == c_Misaligned) {
548 }
else if (set == c_Aligned) {
555 boost::iostreams::filtering_istream ifs;
560 double back[np], fwrd[np], tension;
565 for (
int i = 0; i < np; ++i) {
568 for (
int i = 0; i < np; ++i) {
574 if (ifs.eof())
break;
578 for (
int i = 0; i < np; ++i) {
580 m_BWirPos[iL][iC][i] += back[i];
581 m_FWirPos[iL][iC][i] += fwrd[i];
582 }
else if (set == c_Misaligned) {
583 m_BWirPosMisalign[iL][iC][i] = m_BWirPos[iL][iC][i] + back[i];
584 m_FWirPosMisalign[iL][iC][i] = m_FWirPos[iL][iC][i] + fwrd[i];
585 }
else if (set == c_Aligned) {
586 m_BWirPosAlign[iL][iC][i] = m_BWirPos[iL][iC][i] + back[i];
587 m_FWirPosAlign[iL][iC][i] = m_FWirPos[iL][iC][i] + fwrd[i];
594 m_WireSagCoef[iL][iC] = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.*(m_senseWireTension + tension));
596 }
else if (set == c_Misaligned) {
597 double baseTension = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.* m_WireSagCoef[iL][iC]);
598 m_WireSagCoefMisalign[iL][iC] = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.*
599 (baseTension + tension));
601 }
else if (set == c_Aligned) {
602 double baseTension = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.* m_WireSagCoef[iL][iC]);
603 m_WireSagCoefAlign[iL][iC] = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.*(baseTension + tension));
619 if (nRead != nSenseWires) B2FATAL(
"CDCGeometryPar::readWirePositionParams: #lines read-in (=" << nRead <<
620 ") is inconsistent with total #sense wires (=" << nSenseWires <<
") !");
623 boost::iostreams::close(ifs);
628 void CDCGeometryPar::setWirPosAlignParams()
631 for (
unsigned iL = 0; iL < MAX_N_SLAYERS; ++iL) {
633 auto layerID =
WireID(iL, 511);
636 double d_layerXbwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerX);
637 double d_layerYbwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerY);
638 double d_layerPhiBwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerPhi);
640 double d_layerXfwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerDx) + d_layerXbwd;
641 double d_layerYfwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerDy) + d_layerYbwd;
642 double d_layerPhiFwd = (*m_alignmentFromDB)->get(layerID, CDCAlignment::layerDPhi) + d_layerPhiBwd;
644 for (
unsigned iC = 0; iC < m_nWires[iL]; ++iC) {
646 double wireXbwd = m_BWirPos[iL][iC][0];
647 double wireYbwd = m_BWirPos[iL][iC][1];
648 double wireZbwd = m_BWirPos[iL][iC][2];
650 double wireXfwd = m_FWirPos[iL][iC][0];
651 double wireYfwd = m_FWirPos[iL][iC][1];
652 double wireZfwd = m_FWirPos[iL][iC][2];
656 m_BWirPosAlign[iL][iC][0] = d_layerXbwd + cos(d_layerPhiBwd) * wireXbwd + sin(d_layerPhiBwd) * wireYbwd;
657 m_BWirPosAlign[iL][iC][1] = d_layerYbwd - sin(d_layerPhiBwd) * wireXbwd + cos(d_layerPhiBwd) * wireYbwd;
658 m_BWirPosAlign[iL][iC][2] = wireZbwd;
660 m_FWirPosAlign[iL][iC][0] = d_layerXfwd + cos(d_layerPhiFwd) * wireXfwd + sin(d_layerPhiFwd) * wireYfwd;
661 m_FWirPosAlign[iL][iC][1] = d_layerYfwd - sin(d_layerPhiFwd) * wireXfwd + cos(d_layerPhiFwd) * wireYfwd;
662 m_FWirPosAlign[iL][iC][2] = wireZfwd;
667 double back[np], fwrd[np];
669 for (
unsigned iL = 0; iL < MAX_N_SLAYERS; ++iL) {
670 for (
unsigned iC = 0; iC < m_nWires[iL]; ++iC) {
673 back[0] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireBwdX);
674 back[1] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireBwdY);
675 back[2] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireBwdZ);
677 fwrd[0] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireFwdX);
678 fwrd[1] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireFwdY);
679 fwrd[2] = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireFwdZ);
681 for (
int i = 0; i < np; ++i) {
684 m_BWirPosAlign[iL][iC][i] += back[i];
685 m_FWirPosAlign[iL][iC][i] += fwrd[i];
689 double baseTension = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.* m_WireSagCoef[iL][iC]);
690 double tension = (*m_alignmentFromDB)->get(wire, CDCAlignment::wireTension);
692 m_WireSagCoefAlign[iL][iC] = M_PI * m_senseWireDensity *
693 m_senseWireDiameter * m_senseWireDiameter / (8.*(baseTension + tension));
702 void CDCGeometryPar::setWirPosMisalignParams()
705 double back[np], fwrd[np];
707 for (
unsigned iL = 0; iL < MAX_N_SLAYERS; ++iL) {
708 for (
unsigned iC = 0; iC < m_nWires[iL]; ++iC) {
711 back[0] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireBwdX);
712 back[1] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireBwdY);
713 back[2] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireBwdZ);
715 fwrd[0] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireFwdX);
716 fwrd[1] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireFwdY);
717 fwrd[2] = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireFwdZ);
719 for (
int i = 0; i < np; ++i) {
720 m_BWirPosMisalign[iL][iC][i] = m_BWirPos[iL][iC][i] + back[i];
721 m_FWirPosMisalign[iL][iC][i] = m_FWirPos[iL][iC][i] + fwrd[i];
725 double baseTension = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.* m_WireSagCoef[iL][iC]);
726 double tension = (*m_misalignmentFromDB)->get(wire, CDCMisalignment::wireTension);
728 m_WireSagCoefMisalign[iL][iC] = M_PI * m_senseWireDensity *
729 m_senseWireDiameter * m_senseWireDiameter / (8.*(baseTension + tension));
737 void CDCGeometryPar::readXT(
const GearDir gbxParams,
const int mode)
739 if (m_xtFileFormat == 0) {
742 newReadXT(gbxParams, mode);
748 void CDCGeometryPar::newReadXT(
const GearDir gbxParams,
const int mode)
750 m_linearInterpolationOfXT =
true;
752 std::string fileName0 = CDCGeoControlPar::getInstance().getXtFile();
754 fileName0 = gbxParams.
getString(
"xt4ReconFileName");
757 boost::iostreams::filtering_istream ifs;
781 unsigned short nAlphaBins = 0;
782 if (ifs >> nAlphaBins) {
783 if (nAlphaBins == 0 || nAlphaBins > maxNAlphaPoints) B2FATAL(
"Fail to read alpha bins !");
785 B2FATAL(
"Fail to read alpha bins !");
787 m_nAlphaPoints = nAlphaBins;
788 double alpha0, alpha1, alpha2;
789 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
790 ifs >> alpha0 >> alpha1 >> alpha2;
791 m_alphaPoints[i] = alpha2;
795 unsigned short nThetaBins = 0;
796 if (ifs >> nThetaBins) {
797 if (nThetaBins == 0 || nThetaBins > maxNThetaPoints) B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
799 B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
801 m_nThetaPoints = nThetaBins;
802 double theta0, theta1, theta2;
804 for (
unsigned short i = 0; i < nThetaBins; ++i) {
805 ifs >> theta0 >> theta1 >> theta2;
806 m_thetaPoints[i] = theta2;
810 unsigned short iCL, iLR;
811 const unsigned short npx = nXTParams - 1;
813 double theta, alpha, dummy1;
816 ifs >> m_xtParamMode >> np;
817 if (m_xtParamMode < 0 || m_xtParamMode > 3) B2FATAL(
"CDCGeometryPar: invalid xt-parameterization mode read !");
819 if (np <= 0 || np > npx) B2FATAL(
"CDCGeometryPar: no. of xt-params. outside limits !");
821 const double epsi = 0.1;
824 ifs >> theta >> alpha >> dummy1 >> iLR;
825 for (
int i = 0; i < np; ++i) {
831 for (
unsigned short i = 0; i < nThetaBins; ++i) {
832 if (fabs(theta - m_thetaPoints[i]) < epsi) {
837 if (itheta < 0) B2FATAL(
"CDCGeometryPar: thetas in xt.dat are inconsistent !");
840 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
841 if (fabs(alpha - m_alphaPoints[i]) < epsi) {
846 if (ialpha < 0) B2FATAL(
"CDCGeometryPar: alphas in xt.dat are inconsistent !");
848 for (
int i = 0; i < np; ++i) {
849 m_XT[iCL][iLR][ialpha][itheta][i] = xtc[i];
852 double boundT = xtc[6];
853 if (m_xtParamMode == 1) {
854 m_XT[iCL][iLR][ialpha][itheta][np] = ROOT::Math::Chebyshev5(boundT, xtc[0], xtc[1], xtc[2], xtc[3], xtc[4], xtc[5]);
856 m_XT[iCL][iLR][ialpha][itheta][np] =
867 boost::iostreams::close(ifs);
870 const double degrad = M_PI / 180.;
871 for (
unsigned i = 0; i < nAlphaBins; ++i) {
872 m_alphaPoints[i] *= degrad;
874 for (
unsigned i = 0; i < nThetaBins; ++i) {
875 m_thetaPoints[i] *= degrad;
882 void CDCGeometryPar::readSigma(
const GearDir gbxParams,
const int mode)
884 if (m_sigmaFileFormat == 0) {
887 newReadSigma(gbxParams, mode);
891 void CDCGeometryPar::newReadSigma(
const GearDir gbxParams,
const int mode)
893 m_linearInterpolationOfSgm =
true;
895 std::string fileName0 = CDCGeoControlPar::getInstance().getSigmaFile();
897 fileName0 = gbxParams.
getString(
"sigma4ReconFileName");
905 unsigned short nAlphaBins = 0;
906 if (ifs >> nAlphaBins) {
907 if (nAlphaBins == 0 || nAlphaBins > maxNAlphaPoints) B2FATAL(
"Fail to read alpha bins !");
909 B2FATAL(
"Fail to read alpha bins !");
911 m_nAlphaPoints4Sgm = nAlphaBins;
913 double alpha0, alpha1, alpha2;
914 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
915 ifs >> alpha0 >> alpha1 >> alpha2;
916 m_alphaPoints4Sgm[i] = alpha2;
921 unsigned short nThetaBins = 0;
922 if (ifs >> nThetaBins) {
923 if (nThetaBins == 0 || nThetaBins > maxNThetaPoints) B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
925 B2FATAL(
"CDCGeometryPar: fail to read theta bins !");
927 m_nThetaPoints4Sgm = nThetaBins;
929 double theta0, theta1, theta2;
931 for (
unsigned short i = 0; i < nThetaBins; ++i) {
932 ifs >> theta0 >> theta1 >> theta2;
933 m_thetaPoints4Sgm[i] = theta2;
937 unsigned short np = 0;
938 unsigned short iCL, iLR;
939 double sigma[nSigmaParams];
943 ifs >> m_sigmaParamMode >> np;
945 if (m_sigmaParamMode < 0 || m_sigmaParamMode > 4) B2FATAL(
"CDCGeometryPar: invalid sigma-parameterization mode read !");
947 if (np > nSigmaParams) B2FATAL(
"CDCGeometryPar: no. of sigma-params. outside limits !");
949 const double epsi = 0.1;
952 ifs >> theta >> alpha >> iLR;
954 for (
int i = 0; i < np; ++i) {
960 for (
unsigned short i = 0; i < nThetaBins; ++i) {
961 if (fabs(theta - m_thetaPoints4Sgm[i]) < epsi) {
966 if (itheta < 0) B2FATAL(
"CDCGeometryPar: thetas in sigma.dat are inconsistent !");
969 for (
unsigned short i = 0; i < nAlphaBins; ++i) {
970 if (fabs(alpha - m_alphaPoints4Sgm[i]) < epsi) {
975 if (ialpha < 0) B2FATAL(
"CDCGeometryPar: alphas in sigma.dat are inconsistent !");
977 for (
int i = 0; i < np; ++i) {
978 m_Sigma[iCL][iLR][ialpha][itheta][i] = sigma[i];
985 const double degrad = M_PI / 180.;
986 for (
unsigned i = 0; i < nAlphaBins; ++i) {
987 m_alphaPoints4Sgm[i] *= degrad;
989 for (
unsigned i = 0; i < nThetaBins; ++i) {
990 m_thetaPoints4Sgm[i] *= degrad;
998 void CDCGeometryPar::readFFactor(
const GearDir gbxParams,
const int mode)
1000 std::string fileName0 = CDCGeoControlPar::getInstance().getFFactorFile();
1002 fileName0 = gbxParams.
getString(
"fudgeFactorFileName");
1004 B2WARNING(
"readFFactor is not ready! " << fileName0);
1010 void CDCGeometryPar::readPropSpeed(
const GearDir gbxParams,
const int mode)
1012 std::string fileName0 = CDCGeoControlPar::getInstance().getPropSpeedFile();
1014 fileName0 = gbxParams.
getString(
"propSpeed4ReconFileName");
1027 if (ifs.eof())
break;
1031 m_propSpeedInv[iL] = 1. / speed;
1033 if (m_debug) B2DEBUG(150, iL <<
" " << speed);
1036 if (nRead != MAX_N_SLAYERS) B2FATAL(
"CDCGeometryPar::readPropSpeed: #lines read-in (=" << nRead <<
1037 ") is inconsistent with total #layers (=" << MAX_N_SLAYERS <<
") !");
1078 void CDCGeometryPar::readT0(
const GearDir gbxParams,
int mode)
1080 std::string fileName0 = CDCGeoControlPar::getInstance().getT0File();
1082 fileName0 = gbxParams.
getString(
"t04ReconFileName");
1094 ifs >> iL >> iC >> t0;
1096 if (ifs.eof())
break;
1103 B2DEBUG(150, iL <<
" " << iC <<
" " << t0);
1107 if (nRead != nSenseWires) B2FATAL(
"CDCGeometryPar::readT0: #lines read-in (=" << nRead <<
1108 ") is inconsistent with total #sense wires (=" << nSenseWires <<
") !");
1155 void CDCGeometryPar::readTW(
const GearDir gbxParams,
const int mode)
1157 std::string fileName0 = CDCGeoControlPar::getInstance().getTwFile();
1159 fileName0 = gbxParams.
getString(
"tw4ReconFileName");
1166 unsigned short nPars(0);
1167 ifs >> m_twParamMode >> nPars;
1168 if (m_twParamMode > 1) {
1169 B2FATAL(
"CDCGeometryPar::readTW: invalid mode specified!");
1172 B2FATAL(
"CDCGeometryPar::readTW: invalid #params specified!");
1175 unsigned iBoard = 0;
1178 while (ifs >> iBoard) {
1179 for (
unsigned short i = 0; i < nPars; ++i) {
1180 ifs >> m_timeWalkCoef[iBoard][i];
1185 if (nRead != nBoards) B2FATAL(
"CDCGeometryPar::readTW: #lines read-in (=" << nRead <<
") is inconsistent with #boards (=" << nBoards
1194 void CDCGeometryPar::readChMap()
1196 std::string fileName0 = CDCGeoControlPar::getInstance().getChMapFile();
1202 unsigned short iSL, iL, iW, iB, iC;
1207 ifs >> iSL >> iL >> iW >> iB >> iC;
1208 if (ifs.eof())
break;
1210 if (iSL >= nSuperLayers)
continue;
1213 m_wireToBoard.insert(pair<WireID, unsigned short>(wID, iB));
1216 if (nRead != nSenseWires) B2FATAL(
"CDCGeometryPar::readChMap: #lines read-in (=" << nRead <<
1217 ") is inconsistent with #sense-wires (="
1218 << nSenseWires <<
") !");
1225 void CDCGeometryPar::readEDepToADC(
const GearDir gbxParams,
const int mode)
1228 std::string fileName0 = CDCGeoControlPar::getInstance().getEDepToADCFile();
1230 fileName0 = gbxParams.
getString(
"fudgeFactorFileName");
1235 std::string fileName1 =
"/data/cdc/" + fileName0;
1236 std::string fileName = FileSystem::findFile(fileName1,
true);
1238 if (fileName ==
"") {
1239 fileName = FileSystem::findFile(fileName0,
true);
1242 if (fileName ==
"") {
1243 B2FATAL(
"CDC::openFile: " << fileName0 <<
" not exist!");
1246 B2DEBUG(29,
"CDC::openFile: open " << fileName);
1247 ifs.open(fileName.c_str());
1248 if (!ifs) B2FATAL(
"CDC::openFile: cannot open " << fileName <<
" !");
1251 unsigned short paramMode(0), nParams(0);
1252 ifs >> paramMode >> nParams;
1253 if (paramMode > 0) B2FATAL(
"Param mode > 0!");
1254 if (nParams > 6) B2FATAL(
"No. of params. > 6!");
1255 unsigned short groupId(0);
1257 B2DEBUG(29, paramMode <<
" " << nParams <<
" " << groupId);
1258 if (groupId > 0) B2FATAL(
"GgroupId > 0!");
1261 unsigned short cLMin[nSuperLayers], cLMax[nSuperLayers];
1264 for (
unsigned int sl = 1; sl < nSuperLayers; ++sl) {
1265 cLMin[sl] = cLMax[0] + 6 * sl - 5;
1266 cLMax[sl] = cLMax[0] + 6 * sl;
1268 int nCell[56] = {160, 160, 160, 160, 160, 160, 160, 160,
1269 160, 160, 160, 160, 160, 160,
1270 192, 192, 192, 192, 192, 192,
1271 224, 224, 224, 224, 224, 224,
1272 256, 256, 256, 256, 256, 256,
1273 288, 288, 288, 288, 288, 288,
1274 320, 320, 320, 320, 320, 320,
1275 352, 352, 352, 352, 352, 352,
1276 384, 384, 384, 384, 384, 384
1279 unsigned short id = 0;
1281 unsigned short nRead = 0;
1283 for (
unsigned short i = 0; i < nParams; ++i) {
1285 for (
unsigned short cL = cLMin[
id]; cL <= cLMax[id]; ++cL) {
1286 for (
unsigned short cell = 0; cell < nCell[cL]; ++cell) {
1287 m_eDepToADCParams[cL][cell][i] = coef;
1293 if (nRead > nSuperLayers) B2FATAL(
"No. of read in lines > " << nSuperLayers <<
" !");
1301 void CDCGeometryPar::setT0()
1303 for (
unsigned short iCL = 0; iCL < MAX_N_SLAYERS; ++iCL) {
1304 for (
unsigned short iW = 0; iW < MAX_N_SCELLS; ++iW) {
1309 for (
auto const& ent : (*m_t0FromDB)->getT0s()) {
1312 const unsigned short iW = wid.
getIWire();
1313 m_t0[iCL][iW] = ent.second;
1321 void CDCGeometryPar::calcMeanT0()
1323 B2DEBUG(29,
"calcMeanT0 start");
1324 double effiSum = 0.;
1326 for (
unsigned short iCL = 0; iCL < MAX_N_SLAYERS; ++iCL) {
1327 for (
unsigned short iW = 0; iW < MAX_N_SCELLS; ++iW) {
1328 if (m_t0[iCL][iW] <= 0.)
continue;
1330 if (isHotWire(wid))
continue;
1331 if (isBadWire(wid))
continue;
1334 isDeadWire(wid, effi);
1336 m_meanT0 += effi * m_t0[iCL][iW];
1340 m_meanT0 /= effiSum;
1342 B2FATAL(
"Wire efficiency sum <= 0!");
1344 B2DEBUG(29,
"calcMeanT0 end");
1349 void CDCGeometryPar::setBadWire()
1357 void CDCGeometryPar::setPropSpeed()
1359 for (
unsigned short iCL = 0; iCL < (*m_propSpeedFromDB)->getEntries(); ++iCL) {
1360 m_propSpeedInv[iCL] = 1. / (*m_propSpeedFromDB)->getSpeed(iCL);
1366 void CDCGeometryPar::setTW()
1369 m_twParamMode = (*m_timeWalkFromDB)->getTwParamMode();
1371 for (
unsigned short iBd = 0; iBd < (*m_timeWalkFromDB)->getEntries(); ++iBd) {
1372 int np = ((*m_timeWalkFromDB)->getTimeWalkParams(iBd)).size();
1373 for (
int i = 0; i < np; ++i) {
1374 m_timeWalkCoef[iBd][i] = ((*m_timeWalkFromDB)->getTimeWalkParams(iBd))[i];
1381 void CDCGeometryPar::setXtRel()
1383 m_linearInterpolationOfXT =
true;
1386 m_nAlphaPoints = (*m_xtRelFromDB)->getNoOfAlphaBins();
1387 for (
unsigned short i = 0; i < m_nAlphaPoints; ++i) {
1388 m_alphaPoints[i] = (*m_xtRelFromDB)->getAlphaPoint(i);
1392 m_nThetaPoints = (*m_xtRelFromDB)->getNoOfThetaBins();
1393 for (
unsigned short i = 0; i < m_nThetaPoints; ++i) {
1394 m_thetaPoints[i] = (*m_xtRelFromDB)->getThetaPoint(i);
1398 m_xtParamMode = (*m_xtRelFromDB)->getXtParamMode();
1400 for (
unsigned short iCL = 0; iCL < MAX_N_SLAYERS; ++iCL) {
1401 for (
unsigned short iLR = 0; iLR < 2; ++iLR) {
1402 for (
unsigned short iA = 0; iA < m_nAlphaPoints; ++iA) {
1403 for (
unsigned short iT = 0; iT < m_nThetaPoints; ++iT) {
1404 const std::vector<float> params = (*m_xtRelFromDB)->getXtParams(iCL, iLR, iA, iT);
1405 unsigned short np = params.size();
1407 for (
unsigned short i = 0; i < np; ++i) {
1408 m_XT[iCL][iLR][iA][iT][i] = params[i];
1411 double boundT = m_XT[iCL][iLR][iA][iT][6];
1412 if (m_xtParamMode == 1) {
1413 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],
1414 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]);
1416 m_XT[iCL][iLR][iA][iT][np] =
1417 m_XT[iCL][iLR][iA][iT][0] + boundT
1418 * (m_XT[iCL][iLR][iA][iT][1] + boundT
1419 * (m_XT[iCL][iLR][iA][iT][2] + boundT
1420 * (m_XT[iCL][iLR][iA][iT][3] + boundT
1421 * (m_XT[iCL][iLR][iA][iT][4] + boundT
1422 * (m_XT[iCL][iLR][iA][iT][5])))));
1433 void CDCGeometryPar::setSResol()
1435 m_linearInterpolationOfSgm =
true;
1438 m_nAlphaPoints4Sgm = (*m_sResolFromDB)->getNoOfAlphaBins();
1439 for (
unsigned short i = 0; i < m_nAlphaPoints4Sgm; ++i) {
1440 m_alphaPoints4Sgm[i] = (*m_sResolFromDB)->getAlphaPoint(i);
1444 m_nThetaPoints4Sgm = (*m_sResolFromDB)->getNoOfThetaBins();
1445 for (
unsigned short i = 0; i < m_nThetaPoints4Sgm; ++i) {
1446 m_thetaPoints4Sgm[i] = (*m_sResolFromDB)->getThetaPoint(i);
1453 m_sigmaParamMode = (*m_sResolFromDB)->getSigmaParamMode();
1455 for (
unsigned short iCL = 0; iCL < MAX_N_SLAYERS; ++iCL) {
1456 for (
unsigned short iLR = 0; iLR < 2; ++iLR) {
1457 for (
unsigned short iA = 0; iA < m_nAlphaPoints4Sgm; ++iA) {
1458 for (
unsigned short iT = 0; iT < m_nThetaPoints4Sgm; ++iT) {
1459 const std::vector<float> params = (*m_sResolFromDB)->getSigmaParams(iCL, iLR, iA, iT);
1460 unsigned short np = params.size();
1462 for (
unsigned short i = 0; i < np; ++i) {
1463 m_Sigma[iCL][iLR][iA][iT][i] = params[i];
1474 void CDCGeometryPar::setFFactor()
1476 unsigned short groupId = (*m_fFactorFromDB)->getGroupID();
1477 unsigned short nEnt = (*m_fFactorFromDB)->getEntries();
1478 B2DEBUG(29,
"setFFactor called: groupId,nEnt= " << groupId <<
" " << nEnt);
1482 B2FATAL(
"CDCGeometryPar:: Invalid group-id " << groupId <<
" specified!");
1485 for (
unsigned short id = 0;
id < nEnt; ++id) {
1486 unsigned short np = ((*m_fFactorFromDB)->getFactors(
id)).size();
1487 if (np != 3) B2FATAL(
"CDCGeometryPar:: No. of fudge factors != 3!");
1489 for (
unsigned short i = 0; i < np; ++i) {
1490 m_fudgeFactorForSigma[i] = ((*m_fFactorFromDB)->getFactors(
id))[i];
1491 B2DEBUG(29, i <<
" " << m_fudgeFactorForSigma[i]);
1499 B2DEBUG(29,
"fudge factors= " << m_fudgeFactorForSigma[0] <<
" " << m_fudgeFactorForSigma[1] <<
" " << m_fudgeFactorForSigma[2]);
1504 void CDCGeometryPar::setChMap()
1506 for (
const auto& cm : (*m_chMapFromDB)) {
1507 const unsigned short isl = cm.getISuperLayer();
1508 if (isl >= nSuperLayers)
continue;
1509 const int il = cm.getILayer();
1510 const int iw = cm.getIWire();
1511 const int iBd = cm.getBoardID();
1512 const WireID wID(isl, il, iw);
1513 m_wireToBoard.insert(pair<WireID, unsigned short>(wID, iBd));
1514 const int iCh = cm.getBoardChannel();
1515 m_wireToChannel.insert(pair<WireID, unsigned short>(wID, iCh));
1516 m_boardAndChannelToWire[iBd][iCh] = wID.
getEWire();
1521 void CDCGeometryPar::setEDepToADCConversions()
1523 unsigned short groupId = (*m_eDepToADCConversionsFromDB)->getGroupID();
1524 unsigned short nEnt = (*m_eDepToADCConversionsFromDB)->getEntries();
1526 if (nEnt > nSuperLayers) B2FATAL(
"CDCGeometryPar:: group-id " << groupId <<
" and #entries " << nEnt <<
" are inconsistent!");
1527 }
else if (groupId == 1) {
1528 if (nEnt > MAX_N_SLAYERS) B2FATAL(
"CDCGeometryPar:: group-id " << groupId <<
" and #entries " << nEnt <<
" are inconsistent!");
1530 B2FATAL(
"CDCGeometryPar:: Invalid group-id " << groupId <<
" specified !");
1533 unsigned short cLMin[nSuperLayers], cLMax[nSuperLayers];
1536 for (
unsigned int sl = 1; sl < nSuperLayers; ++sl) {
1537 cLMin[sl] = cLMax[0] + 6 * sl - 5;
1538 cLMax[sl] = cLMax[0] + 6 * sl;
1540 int nCell[56] = {160, 160, 160, 160, 160, 160, 160, 160,
1541 160, 160, 160, 160, 160, 160,
1542 192, 192, 192, 192, 192, 192,
1543 224, 224, 224, 224, 224, 224,
1544 256, 256, 256, 256, 256, 256,
1545 288, 288, 288, 288, 288, 288,
1546 320, 320, 320, 320, 320, 320,
1547 352, 352, 352, 352, 352, 352,
1548 384, 384, 384, 384, 384, 384
1551 for (
unsigned short id = 0;
id < nEnt; ++id) {
1552 unsigned short np = ((*m_eDepToADCConversionsFromDB)->getParams(
id)).size();
1553 if (np > 6) B2FATAL(
"CDCGeometryPar:: No. of edep-to-ADC conversion params. > 6");
1555 for (
unsigned short cL = cLMin[
id]; cL <= cLMax[id]; ++cL) {
1556 for (
unsigned short cell = 0; cell < nCell[cL]; ++cell) {
1557 for (
unsigned short i = 0; i < np; ++i) {
1558 m_eDepToADCParams[cL][cell][i] = ((*m_eDepToADCConversionsFromDB)->getParams(
id))[i];
1562 }
else if (groupId == 1) {
1563 for (
unsigned short cell = 0; cell < nCell[id]; ++cell) {
1564 for (
unsigned short i = 0; i < np; ++i) {
1565 m_eDepToADCParams[id][cell][i] = ((*m_eDepToADCConversionsFromDB)->getParams(
id))[i];
1568 }
else if (groupId == 2) {
1570 B2FATAL(
"CDCGeometryPar::setEDepToADCConversions(): groupId=2 not ready!");
1576 double CDCGeometryPar::getEDepToADCConvFactor(
unsigned short iCL,
unsigned short iW,
double edep,
double dx,
double costh)
1583 const double mainF = m_eDepToADCParams[iCL][iW][0];
1584 const double& alf = m_eDepToADCParams[iCL][iW][1];
1585 const double& gam = m_eDepToADCParams[iCL][iW][2];
1586 const double& dlt = m_eDepToADCParams[iCL][iW][3];
1587 const double& a = m_eDepToADCParams[iCL][iW][4];
1588 const double& b = m_eDepToADCParams[iCL][iW][5];
1589 const double cth = fabs(costh) + dlt;
1590 const double iGen = edep / dx;
1591 const double tmp = cth - gam * iGen;
1592 const double disc = tmp * tmp + 4.*alf * cth * iGen;
1596 iMea = cth * iGen / tmp;
1597 }
else if (disc >= 0.) {
1598 iMea = (-tmp + sqrt(disc)) / (2.*alf);
1606 double convF = mainF;
1608 convF = mainF * std::min(iMea / iGen, 1.);
1617 convF *= 1. + a * (costh - b);
1622 void CDCGeometryPar::Print()
const
1625 const TVector3 CDCGeometryPar::wireForwardPosition(
int layerID,
int cellID,
EWirePosition set)
const
1628 TVector3 wPos(m_FWirPosAlign[layerID][cellID][0],
1629 m_FWirPosAlign[layerID][cellID][1],
1630 m_FWirPosAlign[layerID][cellID][2]);
1632 if (set == c_Misaligned) {
1633 wPos.SetX(m_FWirPosMisalign[layerID][cellID][0]);
1634 wPos.SetY(m_FWirPosMisalign[layerID][cellID][1]);
1635 wPos.SetZ(m_FWirPosMisalign[layerID][cellID][2]);
1636 }
else if (set == c_Base) {
1637 wPos.SetX(m_FWirPos [layerID][cellID][0]);
1638 wPos.SetY(m_FWirPos [layerID][cellID][1]);
1639 wPos.SetZ(m_FWirPos [layerID][cellID][2]);
1644 const TVector3 CDCGeometryPar::wireForwardPosition(
int layerID,
int cellID,
double z,
EWirePosition set)
const
1648 getWireSagEffect(set, layerID, cellID, z, yb_sag, yf_sag);
1650 TVector3 wPos(m_FWirPosAlign[layerID][cellID][0], yf_sag,
1651 m_FWirPosAlign[layerID][cellID][2]);
1652 if (set == c_Misaligned) {
1653 wPos.SetX(m_FWirPosMisalign[layerID][cellID][0]);
1654 wPos.SetZ(m_FWirPosMisalign[layerID][cellID][2]);
1655 }
else if (set == c_Base) {
1656 wPos.SetX(m_FWirPos [layerID][cellID][0]);
1657 wPos.SetZ(m_FWirPos [layerID][cellID][2]);
1662 const TVector3 CDCGeometryPar::wireBackwardPosition(
int layerID,
int cellID,
EWirePosition set)
const
1664 TVector3 wPos(m_BWirPosAlign[layerID][cellID][0],
1665 m_BWirPosAlign[layerID][cellID][1],
1666 m_BWirPosAlign[layerID][cellID][2]);
1668 if (set == c_Misaligned) {
1669 wPos.SetX(m_BWirPosMisalign[layerID][cellID][0]);
1670 wPos.SetY(m_BWirPosMisalign[layerID][cellID][1]);
1671 wPos.SetZ(m_BWirPosMisalign[layerID][cellID][2]);
1672 }
else if (set == c_Base) {
1673 wPos.SetX(m_BWirPos [layerID][cellID][0]);
1674 wPos.SetY(m_BWirPos [layerID][cellID][1]);
1675 wPos.SetZ(m_BWirPos [layerID][cellID][2]);
1680 const TVector3 CDCGeometryPar::wireBackwardPosition(
int layerID,
int cellID,
double z,
EWirePosition set)
const
1684 getWireSagEffect(set, layerID, cellID, z, yb_sag, yf_sag);
1686 TVector3 wPos(m_BWirPosAlign[layerID][cellID][0], yb_sag,
1687 m_BWirPosAlign[layerID][cellID][2]);
1688 if (set == c_Misaligned) {
1689 wPos.SetX(m_BWirPosMisalign[layerID][cellID][0]);
1690 wPos.SetZ(m_BWirPosMisalign[layerID][cellID][2]);
1691 }
else if (set == c_Base) {
1692 wPos.SetX(m_BWirPos [layerID][cellID][0]);
1693 wPos.SetZ(m_BWirPos [layerID][cellID][2]);
1698 double CDCGeometryPar::getWireSagCoef(
EWirePosition set,
int layerID,
int cellID)
const
1700 double coef = m_WireSagCoef[layerID][cellID];
1701 if (set == c_Misaligned) {
1702 coef = m_WireSagCoefMisalign[layerID][cellID];
1703 }
else if (set == c_Aligned) {
1704 coef = m_WireSagCoefAlign [layerID][cellID];
1709 const double* CDCGeometryPar::innerRadiusWireLayer()
const
1711 static double IRWL[MAX_N_SLAYERS] = {0};
1713 IRWL[0] = outerRadiusInnerWall();
1714 for (
unsigned i = 1; i < nWireLayers(); i++)
1716 IRWL[i] = m_rFLayer[i - 1];
1721 const double* CDCGeometryPar::outerRadiusWireLayer()
const
1723 static double ORWL[MAX_N_SLAYERS] = {0};
1725 ORWL[nWireLayers() - 1] = innerRadiusOuterWall();
1726 for (
unsigned i = 0; i < nWireLayers() - 1; i++)
1728 ORWL[i] = m_rFLayer[i];
1733 unsigned CDCGeometryPar::cellId(
unsigned layerId,
const TVector3& position)
const
1735 const unsigned nWires = m_nWires[layerId];
1737 double offset = m_offSet[layerId];
1739 const double phiSize = 2 * M_PI / double(nWires);
1756 for (
unsigned i = 0; i < 1; ++i) {
1757 const double phiF = phiSize * (double(i) + offset)
1758 + phiSize * 0.5 *
double(m_nShifts[layerId]) + m_globalPhiRotation;
1759 const double phiB = phiSize * (double(i) + offset) + m_globalPhiRotation;
1760 const TVector3 f(m_rSLayer[layerId] * cos(phiF), m_rSLayer[layerId] * sin(phiF), m_zSForwardLayer[layerId]);
1761 const TVector3 b(m_rSLayer[layerId] * cos(phiB), m_rSLayer[layerId] * sin(phiB), m_zSBackwardLayer[layerId]);
1762 const TVector3 v = f - b;
1763 const TVector3 u = v.Unit();
1764 const double beta = (position.z() - b.z()) / u.z();
1765 const TVector3 p = b + beta * u;
1766 double dPhi = std::atan2(position.y(), position.x())
1767 - std::atan2(p.y(), p.x())
1769 while (dPhi < 0) dPhi += (2. * M_PI);
1770 j = int(dPhi / phiSize);
1771 while (j >= nWires) j -= nWires;
1777 void CDCGeometryPar::generateXML(
const string& of)
1780 std::ofstream ofs(of.c_str(), std::ios::out);
1782 B2ERROR(
"CDCGeometryPar::read !!! can not open file : "
1785 ofs <<
"<?xml version=\"1.0\" encoding=\"UTF-8\"?>"
1787 <<
"<Subdetector type=\"CDC\">"
1789 <<
" <Name>CDC BelleII </Name>"
1791 <<
" <Description>CDC geometry parameters</Description>"
1793 <<
" <Version>0</Version>"
1795 <<
" <GeoCreator>CDCBelleII</GeoCreator>"
1799 <<
" <Rotation desc=\"Rotation of the whole cdc detector (should be the same as beampipe)\" unit=\"mrad\">0.0</Rotation>"
1801 <<
" <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>"
1803 <<
" <Material>CDCGas</Material>"
1807 ofs <<
" <SLayers>" << endl;
1809 for (
int i = 0; i < m_nSLayer; i++) {
1810 ofs <<
" <SLayer id=\"" << i <<
"\">" << endl;
1811 ofs <<
" <Radius desc=\"Radius of wires in this layer\" unit=\"mm\">" << senseWireR(i) <<
"</Radius>" << endl;
1812 ofs <<
" <BackwardZ desc=\"z position of this wire layer at backward endplate\" unit=\"mm\">" << senseWireBZ(
1813 i) <<
"</BackwardZ>" << endl;
1814 ofs <<
" <ForwardZ desc=\"z position of this wire layer at forward endplate\" unit=\"mm\">" << senseWireFZ(
1815 i) <<
"</ForwardZ>" << endl;
1818 ofs <<
" <NHoles desc=\"the number of holes in this layer, 2*(cell number)\">" << nWiresInLayer(
1819 i) * 2 <<
"</NHoles>" << endl;
1820 ofs <<
" <NShift desc=\"the shifted hole number of each wire in this layer\">" << nShifts(i) <<
"</NShift>" << endl;
1821 ofs <<
" <Offset desc=\"wire offset in phi direction at endplate\">" << m_offSet[i] <<
"</Offset>" << endl;
1822 ofs <<
" </SLayer>" << endl;
1825 ofs <<
" </SLayers>" << endl;
1826 ofs <<
" <FLayers>" << endl;
1828 for (
int i = 0; i < m_nFLayer; i++) {
1829 ofs <<
" <FLayer id=\"" << i <<
"\">" << endl;
1830 ofs <<
" <Radius desc=\"Radius of field wires in this layer\" unit=\"mm\">" << fieldWireR(i) <<
"</Radius>" << endl;
1831 ofs <<
" <BackwardZ desc=\"z position of this field wire layer at backward endplate\" unit=\"mm\">" << fieldWireBZ(
1832 i) <<
"</BackwardZ>" << endl;
1833 ofs <<
" <ForwardZ desc=\"z position of this field wire layer at forward endplate\" unit=\"mm\">" << fieldWireFZ(
1834 i) <<
"</ForwardZ>" << endl;
1835 ofs <<
" </FLayer>" << endl;
1838 ofs <<
" </FLayers>" << endl;
1840 ofs <<
" <InnerWall name=\"InnerWall\">" << endl;
1841 ofs <<
" <InnerR desc=\"Inner radius\" unit=\"mm\">" << innerRadiusInnerWall() <<
"</InnerR>" << endl;
1842 ofs <<
" <OuterR desc=\"Outer radius\" unit=\"mm\">" << outerRadiusInnerWall() <<
"</OuterR>" << endl;
1843 ofs <<
" <BackwardZ desc=\"z position at backward endplate\" unit=\"mm\">" << m_zWall[0][0] <<
"</BackwardZ>" << endl;
1844 ofs <<
" <ForwardZ desc=\"z position at forward endplate\" unit=\"mm\">" << m_zWall[0][1] <<
"</ForwardZ>" << endl;
1845 ofs <<
" </InnerWall>" << endl;
1847 ofs <<
" <OuterWall name=\"OuterWall\">" << endl;
1848 ofs <<
" <InnerR desc=\"Inner radius\" unit=\"mm\">" << innerRadiusOuterWall() <<
"</InnerR>" << endl;
1849 ofs <<
" <OuterR desc=\"Outer radius\" unit=\"mm\">" << outerRadiusOuterWall() <<
"</OuterR>" << endl;
1850 ofs <<
" <BackwardZ desc=\"z position at backward endplate\" unit=\"mm\">" << m_zWall[2][0] <<
"</BackwardZ>" << endl;
1851 ofs <<
" <ForwardZ desc=\"z position at forward endplate\" unit=\"mm\">" << m_zWall[2][1] <<
"</ForwardZ>" << endl;
1852 ofs <<
" </OuterWall>" << endl;
1854 ofs <<
" </Content>" << endl
1855 <<
"</Subdetector>" << endl;
1858 void CDCGeometryPar::getWireSagEffect(
const EWirePosition set,
const unsigned layerID,
const unsigned cellID,
const double Z,
1859 double& Yb_sag,
double& Yf_sag)
const
1883 if (set == c_Aligned) {
1884 Coef = m_WireSagCoefAlign[layerID][cellID];
1885 Yb = m_BWirPosAlign[layerID][cellID][1];
1886 Yf = m_FWirPosAlign[layerID][cellID][1];
1892 Xb = m_BWirPosAlign[layerID][cellID][0];
1893 Xf = m_FWirPosAlign[layerID][cellID][0];
1894 Zb = m_BWirPosAlign[layerID][cellID][2];
1895 Zf = m_FWirPosAlign[layerID][cellID][2];
1897 }
else if (set == c_Misaligned) {
1898 Coef = m_WireSagCoefMisalign[layerID][cellID];
1899 Yb = m_BWirPosMisalign[layerID][cellID][1];
1900 Yf = m_FWirPosMisalign[layerID][cellID][1];
1906 Xb = m_BWirPosMisalign[layerID][cellID][0];
1907 Xf = m_FWirPosMisalign[layerID][cellID][0];
1908 Zb = m_BWirPosMisalign[layerID][cellID][2];
1909 Zf = m_FWirPosMisalign[layerID][cellID][2];
1911 }
else if (set == c_Base) {
1912 Coef = m_WireSagCoef[layerID][cellID];
1913 Yb = m_BWirPos[layerID][cellID][1];
1914 Yf = m_FWirPos[layerID][cellID][1];
1920 Xb = m_BWirPos[layerID][cellID][0];
1921 Xf = m_FWirPos[layerID][cellID][0];
1922 Zb = m_BWirPos[layerID][cellID][2];
1923 Zf = m_FWirPos[layerID][cellID][2];
1926 B2FATAL(
"CDCGeometryPar::getWireSagEffect: called with an invalid set: " <<
" " << set);
1929 const double dx = Xf - Xb;
1930 const double dy = Yf - Yb;
1931 const double dz = Zf - Zb;
1933 const double Zfp = sqrt(dz * dz + dx * dx);
1934 const double Zp = (Z - Zb) * Zfp / dz;
1936 const double Y_sag = (Coef * (Zp - Zfp) + dy / Zfp) * Zp + Yb;
1937 const double dydz = (Coef * (2.*Zp - Zfp) * Zfp + dy) / dz;
1939 Yb_sag = Y_sag + dydz * (Zb - Z);
1940 Yf_sag = Y_sag + dydz * (Zf - Z);
1944 void CDCGeometryPar::setDesignWirParam(
const unsigned layerID,
const unsigned cellID)
1946 const unsigned L = layerID;
1947 const unsigned C = cellID;
1949 const double offset = m_offSet[L];
1951 const double phiSize = 2 * M_PI / double(m_nWires[L]);
1953 const double phiF = phiSize * (double(C) + offset)
1954 + phiSize * 0.5 *
double(m_nShifts[L]) + m_globalPhiRotation;
1956 m_FWirPos[L][C][0] = m_rSLayer[L] * cos(phiF);
1957 m_FWirPos[L][C][1] = m_rSLayer[L] * sin(phiF);
1958 m_FWirPos[L][C][2] = m_zSForwardLayer[L];
1960 const double phiB = phiSize * (double(C) + offset) + m_globalPhiRotation;
1962 m_BWirPos[L][C][0] = m_rSLayer[L] * cos(phiB);
1963 m_BWirPos[L][C][1] = m_rSLayer[L] * sin(phiB);
1964 m_BWirPos[L][C][2] = m_zSBackwardLayer[L];
1966 for (
int i = 0; i < 3; ++i) {
1967 m_FWirPosMisalign[L][C][i] = m_FWirPos[L][C][i];
1968 m_BWirPosMisalign[L][C][i] = m_BWirPos[L][C][i];
1969 m_FWirPosAlign [L][C][i] = m_FWirPos[L][C][i];
1970 m_BWirPosAlign [L][C][i] = m_BWirPos[L][C][i];
1973 m_WireSagCoef[L][C] = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8. * m_senseWireTension);
1975 m_WireSagCoefMisalign[L][C] = m_WireSagCoef[L][C];
1976 m_WireSagCoefAlign [L][C] = m_WireSagCoef [L][C];
1980 void CDCGeometryPar::outputDesignWirParam(
const unsigned layerID,
const unsigned cellID)
const
1983 const unsigned L = layerID;
1984 const unsigned C = cellID;
1986 static bool first =
true;
1987 static ofstream ofs;
1990 ofs.open(
"alignment.dat");
1993 ofs << L <<
" " << C;
1995 ofs << setiosflags(ios::showpoint | ios::uppercase);
1997 for (
int i = 0; i < 3; ++i) ofs <<
" " << setw(15) << setprecision(8) << m_BWirPos[L][C][i];
1999 for (
int i = 0; i < 3; ++i) ofs <<
" " << setw(15) << setprecision(8) << m_FWirPos[L][C][i];
2000 ofs << setiosflags(ios::fixed);
2001 ofs <<
" " << setw(4) << setprecision(1) << m_senseWireTension;
2006 double CDCGeometryPar::getDriftV(
const double time,
const unsigned short iCLayer,
const unsigned short lr,
const double alpha,
2007 const double theta)
const
2012 double minTime = getMinDriftTime(iCLayer, lr, alpha, theta);
2013 double delta = time - minTime;
2016 unsigned short lro = getOutgoingLR(lr, alpha);
2018 if (!m_linearInterpolationOfXT) {
2019 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2022 unsigned short ial[2] = {0};
2023 unsigned short ilr[2] = {lro, lro};
2024 getClosestAlphaPoints(alpha, wal, ial, ilr);
2026 unsigned short ith[2] = {0};
2027 getClosestThetaPoints(alpha, theta, wth, ith);
2029 unsigned short jal(0), jlr(0), jth(0);
2033 double timep = delta < 0. ? minTime - delta : time;
2036 for (
unsigned k = 0; k < 4; ++k) {
2041 w = (1. - wal) * (1. - wth);
2042 }
else if (k == 1) {
2046 w = (1. - wal) * wth;
2047 }
else if (k == 2) {
2051 w = wal * (1. - wth);
2052 }
else if (k == 3) {
2059 double boundary = m_XT[iCLayer][jlr][jal][jth][6];
2061 if (timep < boundary) {
2062 if (m_xtParamMode == 1) {
2063 const double& c1 = m_XT[iCLayer][jlr][jal][jth][1];
2064 const double& c2 = m_XT[iCLayer][jlr][jal][jth][2];
2065 const double& c3 = m_XT[iCLayer][jlr][jal][jth][3];
2066 const double& c4 = m_XT[iCLayer][jlr][jal][jth][4];
2067 const double& c5 = m_XT[iCLayer][jlr][jal][jth][5];
2068 dDdt += w * ROOT::Math::Chebyshev4(timep, c1 + 3.*c3 + 5.*c5, 4.*c2 + 8.*c4, 6.*c3 + 10.*c5, 8.*c4, 10.*c5);
2070 dDdt += w * (m_XT[iCLayer][jlr][jal][jth][1] + timep
2071 * (2.*m_XT[iCLayer][jlr][jal][jth][2] + timep
2072 * (3.*m_XT[iCLayer][jlr][jal][jth][3] + timep
2073 * (4.*m_XT[iCLayer][jlr][jal][jth][4] + timep
2074 * (5.*m_XT[iCLayer][jlr][jal][jth][5])))));
2077 dDdt += w * m_XT[iCLayer][jlr][jal][jth][7];
2090 double CDCGeometryPar::getDriftLength0(
const double time,
const unsigned short iCLayer,
const unsigned short lr,
const double alpha,
2091 const double theta)
const
2096 unsigned short lro = getOutgoingLR(lr, alpha);
2100 if (!m_linearInterpolationOfXT) {
2101 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2104 unsigned short ial[2] = {0};
2105 unsigned short ilr[2] = {lro, lro};
2106 getClosestAlphaPoints(alpha, wal, ial, ilr);
2108 unsigned short ith[2] = {0};
2109 getClosestThetaPoints(alpha, theta, wth, ith);
2111 unsigned short jal(0), jlr(0), jth(0);
2115 double timep = time;
2119 for (
unsigned k = 0; k < 4; ++k) {
2124 w = (1. - wal) * (1. - wth);
2125 }
else if (k == 1) {
2129 w = (1. - wal) * wth;
2130 }
else if (k == 2) {
2134 w = wal * (1. - wth);
2135 }
else if (k == 3) {
2142 double boundary = m_XT[iCLayer][jlr][jal][jth][6];
2144 if (timep < boundary) {
2145 if (m_xtParamMode == 1) {
2146 dist += w * ROOT::Math::Chebyshev5(timep, m_XT[iCLayer][jlr][jal][jth][0], m_XT[iCLayer][jlr][jal][jth][1],
2147 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]);
2149 dist += w * (m_XT[iCLayer][jlr][jal][jth][0] + timep
2150 * (m_XT[iCLayer][jlr][jal][jth][1] + timep
2151 * (m_XT[iCLayer][jlr][jal][jth][2] + timep
2152 * (m_XT[iCLayer][jlr][jal][jth][3] + timep
2153 * (m_XT[iCLayer][jlr][jal][jth][4] + timep
2154 * (m_XT[iCLayer][jlr][jal][jth][5]))))));
2157 dist += w * (m_XT[iCLayer][jlr][jal][jth][7] * (timep - boundary) + m_XT[iCLayer][jlr][jal][jth][8]);
2168 double CDCGeometryPar::getDriftLength(
const double time,
const unsigned short iCLayer,
const unsigned short lr,
const double alpha,
2170 const bool calculateMinTime,
2171 const double inputMinTime)
const
2176 double minTime = calculateMinTime ? getMinDriftTime(iCLayer, lr, alpha, theta) : inputMinTime;
2177 double delta = time - minTime;
2180 unsigned short lro = getOutgoingLR(lr, alpha);
2184 if (!m_linearInterpolationOfXT) {
2185 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2188 unsigned short ial[2] = {0};
2189 unsigned short ilr[2] = {lro, lro};
2190 getClosestAlphaPoints(alpha, wal, ial, ilr);
2192 unsigned short ith[2] = {0};
2193 getClosestThetaPoints(alpha, theta, wth, ith);
2195 unsigned short jal(0), jlr(0), jth(0);
2199 double timep = delta < 0. ? minTime - delta : time;
2204 for (
unsigned k = 0; k < 4; ++k) {
2209 w = (1. - wal) * (1. - wth);
2210 }
else if (k == 1) {
2214 w = (1. - wal) * wth;
2215 }
else if (k == 2) {
2219 w = wal * (1. - wth);
2220 }
else if (k == 3) {
2240 double boundary = m_XT[iCLayer][jlr][jal][jth][6];
2242 if (timep < boundary) {
2243 if (m_xtParamMode == 1) {
2244 dist += w * ROOT::Math::Chebyshev5(timep, m_XT[iCLayer][jlr][jal][jth][0], m_XT[iCLayer][jlr][jal][jth][1],
2245 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]);
2247 dist += w * (m_XT[iCLayer][jlr][jal][jth][0] + timep
2248 * (m_XT[iCLayer][jlr][jal][jth][1] + timep
2249 * (m_XT[iCLayer][jlr][jal][jth][2] + timep
2250 * (m_XT[iCLayer][jlr][jal][jth][3] + timep
2251 * (m_XT[iCLayer][jlr][jal][jth][4] + timep
2252 * (m_XT[iCLayer][jlr][jal][jth][5]))))));
2255 dist += w * (m_XT[iCLayer][jlr][jal][jth][7] * (timep - boundary) + m_XT[iCLayer][jlr][jal][jth][8]);
2262 if (delta < 0.) dist *= -1.;
2267 double CDCGeometryPar::getMinDriftTime(
const unsigned short iCLayer,
const unsigned short lr,
const double alpha,
2268 const double theta)
const
2270 double minTime = 0.;
2273 unsigned short lro = getOutgoingLR(lr, alpha);
2275 if (!m_linearInterpolationOfXT) {
2276 B2FATAL(
"linearInterpolationOfXT = false is not allowed now !");
2279 unsigned short ial[2] = {0};
2280 unsigned short ilr[2] = {lro, lro};
2281 getClosestAlphaPoints(alpha, wal, ial, ilr);
2283 unsigned short ith[2] = {0};
2284 getClosestThetaPoints(alpha, theta, wth, ith);
2286 unsigned short jal(0), jlr(0), jth(0);
2289 double c[6] = {0.}, a[6] = {0.};
2290 for (
unsigned k = 0; k < 4; ++k) {
2295 w = (1. - wal) * (1. - wth);
2296 }
else if (k == 1) {
2300 w = (1. - wal) * wth;
2301 }
else if (k == 2) {
2305 w = wal * (1. - wth);
2306 }
else if (k == 3) {
2313 for (
int i = 0; i < 5; ++i) {
2314 c[i] += w * m_XT[iCLayer][jlr][jal][jth][i];
2318 if (m_xtParamMode == 1) {
2319 a[0] = c[0] - c[2] + c[4];
2320 a[1] = c[1] - 3.*c[3] + 5.*c[5];
2321 a[2] = 2.*c[2] - 8.*c[4];
2322 a[3] = 4.*c[3] - 20.*c[5];
2326 for (
int i = 0; i < 5; ++i) a[i] = c[i];
2333 const double det = a[1] * a[1] - 4.*a[2] * a[0];
2336 minTime = (-a[1] + sqrt(det)) / (2.*a[2]);
2340 minTime = -a[1] / (2.*a[2]);
2343 }
else if (a[1] != 0.) {
2344 minTime = -a[0] / a[1];
2346 B2WARNING(
"CDCGeometryPar::getMinDriftTime: minDriftTime not determined; assume zero.\n" <<
"layer(#0-55),lr,alpha(rad),theta= " <<
2347 iCLayer <<
" " << lr <<
" " << alpha <<
" " << theta);
2356 const double epsi4x = 5.e-6;
2358 const unsigned short maxIter = 8;
2359 const double maxDt = 20.;
2360 unsigned short nIter = 0;
2361 double minXsq = 1.e10;
2362 double minMinTime = minTime;
2365 for (nIter = 0; nIter <= maxIter; ++nIter) {
2368 double x = a[0] + t * (a[1] + t * (a[2] + t * (a[3] + t * (a[4] + t * a[5]))));
2374 double xp = a[1] + t * (2 * a[2] + t * (3 * a[3] + t * (4 * a[4] + t * 5 * a[5])));
2375 double xpp = 2 * a[2] + t * (6 * a[3] + t * (12 * a[4] + t * 20 * a[5]));
2376 double den = xp * xp + x * xpp;
2383 edm = fabs(x * xp) / sqrt(den);
2384 if (edm < epsi4x)
break;
2391 dt = std::min(dt, maxDt);
2393 dt = std::max(dt, -maxDt);
2396 B2WARNING(
"CDCGeometryPar::getMinDriftTime: den = 0\n" <<
"layer(#0-55),lr,alpha(rad),theta= " <<
2399 " " << alpha <<
" " << theta);
2405 if (nIter == (maxIter + 1)) minTime = minMinTime;
2457 double CDCGeometryPar::getDriftTime(
const double dist,
const unsigned short iCLayer,
const unsigned short lr,
const double alpha,
2458 const double theta)
const
2462 const double eps = 2.5e-1;
2463 const double maxTrials = 100;
2471 double maxTime = 2000.;
2476 double minTime = getMinDriftTime(iCLayer, lr, alpha, theta);
2477 double t0 = minTime;
2479 const bool calMinTime =
false;
2484 double t1 = maxTime;
2485 double time = dist * m_nominalDriftVInv;
2486 while (((t1 - t0) > eps) && (i < maxTrials)) {
2487 time = 0.5 * (t0 + t1);
2488 double d1 = getDriftLength(time, iCLayer, lr, alpha, theta, calMinTime, minTime) - dist;
2498 if (i >= maxTrials - 1 || time > maxTime) {
2499 B2WARNING(
"CDCGeometryPar::getDriftTime " << dist <<
" " << iCLayer <<
" " << alpha <<
" " << lr <<
" " << t0 <<
" " << t1 <<
" " <<
2509 double CDCGeometryPar::getSigma(
const double DriftL0,
const unsigned short iCLayer,
const unsigned short lr,
const double alpha,
2510 const double theta)
const
2515 const double driftL = fabs(DriftL0);
2518 unsigned short lro = getOutgoingLR(lr, alpha);
2520 if (!m_linearInterpolationOfSgm) {
2521 B2FATAL(
"linearInterpolationOfSgm = false is not allowed now !");
2523 if (m_linearInterpolationOfSgm) {
2525 unsigned short ial[2] = {0};
2526 unsigned short ilr[2] = {lro, lro};
2527 getClosestAlphaPoints4Sgm(alpha, wal, ial, ilr);
2529 unsigned short ith[2] = {0};
2530 getClosestThetaPoints4Sgm(alpha, theta, wth, ith);
2533 unsigned short jal(0), jlr(0), jth(0);
2535 for (
unsigned k = 0; k < 4; ++k) {
2540 w = (1. - wal) * (1. - wth);
2541 }
else if (k == 1) {
2545 w = (1. - wal) * wth;
2546 }
else if (k == 2) {
2550 w = wal * (1. - wth);
2551 }
else if (k == 3) {
2569 const double& P0 = m_Sigma[iCLayer][jlr][jal][jth][0];
2570 const double& P1 = m_Sigma[iCLayer][jlr][jal][jth][1];
2571 const double& P2 = m_Sigma[iCLayer][jlr][jal][jth][2];
2572 const double& P3 = m_Sigma[iCLayer][jlr][jal][jth][3];
2573 const double& P4 = m_Sigma[iCLayer][jlr][jal][jth][4];
2574 const double& P5 = m_Sigma[iCLayer][jlr][jal][jth][5];
2575 const double& P6 = m_Sigma[iCLayer][jlr][jal][jth][6];
2577 #if defined(CDC_DEBUG)
2578 cout <<
"driftL= " << driftL << endl;
2579 cout <<
"iCLayer= " << iCLayer <<
" " << jlr <<
" " << jal <<
" " << jth << endl;
2580 cout <<
"P0= " << P0 << endl;
2581 cout <<
"P1= " << P1 << endl;
2582 cout <<
"P2= " << P2 << endl;
2583 cout <<
"P3= " << P3 << endl;
2584 cout <<
"P4= " << P4 << endl;
2585 cout <<
"P5= " << P5 << endl;
2586 cout <<
"P6= " << P6 << endl;
2588 const double P7 = m_sigmaParamMode == 0 ? DBL_MAX : m_Sigma[iCLayer][jlr][jal][jth][7];
2591 sigma += w * sqrt(P0 / (driftL * driftL + P1) + P2 * driftL + P3 +
2592 P4 * exp(P5 * (driftL - P6) * (driftL - P6)));
2594 double forthTermAtP7 = P4 * exp(P5 * (P7 - P6) * (P7 - P6));
2595 const double& P8 = m_Sigma[iCLayer][jlr][jal][jth][8];
2596 if (m_sigmaParamMode == 1) {
2597 double sigmaAtP7 = sqrt(P0 / (P7 * P7 + P1) + P2 * P7 + P3 + forthTermAtP7);
2598 sigma += w * (P8 * (driftL - P7) + sigmaAtP7);
2599 }
else if (m_sigmaParamMode == 2) {
2600 double onePls4AtP7 = sqrt(P0 / (P7 * P7 + P1) + forthTermAtP7);
2601 const double onePls4 = P8 * (driftL - P7) + onePls4AtP7;
2602 sigma += w * sqrt(P2 * driftL + P3 + onePls4 * onePls4);
2603 }
else if (m_sigmaParamMode == 3) {
2604 forthTermAtP7 = sqrt(forthTermAtP7);
2605 const double forthTerm = P8 * (driftL - P7) + forthTermAtP7;
2606 sigma += w * sqrt(P0 / (driftL * driftL + P1) + P2 * driftL + P3 +
2607 forthTerm * forthTerm);
2613 sigma = std::min(sigma, m_maxSpaceResol);
2617 unsigned short CDCGeometryPar::getOldLeftRight(
const TVector3& posOnWire,
const TVector3& posOnTrack,
2618 const TVector3& momentum)
const
2620 unsigned short lr = 0;
2621 double wCrossT = (posOnWire.Cross(posOnTrack)).z();
2625 }
else if (wCrossT > 0.) {
2628 if ((posOnTrack - posOnWire).Perp() != 0.) {
2629 double wCrossP = (posOnWire.Cross(momentum)).z();
2631 if (posOnTrack.Perp() > posOnWire.Perp()) {
2636 }
else if (wCrossP < 0.) {
2637 if (posOnTrack.Perp() < posOnWire.Perp()) {
2652 unsigned short CDCGeometryPar::getNewLeftRightRaw(
const TVector3& posOnWire,
const TVector3& posOnTrack,
2653 const TVector3& momentum)
const
2655 const double distanceCrossP = ((posOnWire - posOnTrack).Cross(momentum)).z();
2656 unsigned short int lr = (distanceCrossP > 0.) ? 1 : 0;
2661 double CDCGeometryPar::getAlpha(
const TVector3& posOnWire,
const TVector3& momentum)
const
2663 const double wx = posOnWire.x();
2664 const double wy = posOnWire.y();
2665 const double px = momentum.x();
2666 const double py = momentum.y();
2668 const double cross = wx * py - wy * px;
2669 const double dot = wx * px + wy * py;
2671 return atan2(cross, dot);
2674 double CDCGeometryPar::getTheta(
const TVector3& momentum)
const
2676 return atan2(momentum.Perp(), momentum.z());
2680 unsigned short CDCGeometryPar::getOutgoingLR(
const unsigned short lr,
const double alpha)
const
2682 unsigned short lro = (fabs(alpha) <= 0.5 * M_PI) ? lr : abs(lr - 1);
2687 double CDCGeometryPar::getOutgoingAlpha(
const double alpha)
const
2690 double alphao = alpha;
2691 if (alpha > 0.5 * M_PI) {
2693 }
else if (alpha < -0.5 * M_PI) {
2700 double CDCGeometryPar::getOutgoingTheta(
const double alpha,
const double theta)
const
2703 double thetao = fabs(alpha) > 0.5 * M_PI ? M_PI - theta : theta;
2708 void CDCGeometryPar::getClosestAlphaPoints(
const double alpha,
double& weight,
unsigned short points[2],
2709 unsigned short lrs[2])
const
2711 double alphao = getOutgoingAlpha(alpha);
2714 if (alphao < m_alphaPoints[0]) {
2715 points[0] = m_nAlphaPoints - 1;
2717 if (m_nAlphaPoints > 1) {
2718 lrs[0] = abs(lrs[0] - 1);
2719 weight = (alphao - (m_alphaPoints[points[0]] - M_PI)) / (m_alphaPoints[points[1]] - (m_alphaPoints[points[0]] - M_PI));
2721 }
else if (m_alphaPoints[m_nAlphaPoints - 1] <= alphao) {
2722 points[0] = m_nAlphaPoints - 1;
2724 if (m_nAlphaPoints > 1) {
2725 lrs[1] = abs(lrs[1] - 1);
2726 weight = (alphao - m_alphaPoints[points[0]]) / (m_alphaPoints[points[1]] + M_PI - m_alphaPoints[points[0]]);
2729 for (
unsigned short i = 0; i <= m_nAlphaPoints - 2; ++i) {
2730 if (m_alphaPoints[i] <= alphao && alphao < m_alphaPoints[i + 1]) {
2733 weight = (alphao - m_alphaPoints[points[0]]) / (m_alphaPoints[points[1]] - m_alphaPoints[points[0]]);
2742 void CDCGeometryPar::getClosestAlphaPoints4Sgm(
const double alpha,
double& weight,
unsigned short points[2],
2743 unsigned short lrs[2])
const
2745 double alphao = getOutgoingAlpha(alpha);
2748 if (alphao < m_alphaPoints4Sgm[0]) {
2749 points[0] = m_nAlphaPoints4Sgm - 1;
2751 if (m_nAlphaPoints4Sgm > 1) {
2752 lrs[0] = abs(lrs[0] - 1);
2753 weight = (alphao - (m_alphaPoints4Sgm[points[0]] - M_PI)) / (m_alphaPoints4Sgm[points[1]] - (m_alphaPoints4Sgm[points[0]] - M_PI));
2755 }
else if (m_alphaPoints4Sgm[m_nAlphaPoints4Sgm - 1] <= alphao) {
2756 points[0] = m_nAlphaPoints4Sgm - 1;
2758 if (m_nAlphaPoints4Sgm > 1) {
2759 lrs[1] = abs(lrs[1] - 1);
2760 weight = (alphao - m_alphaPoints4Sgm[points[0]]) / (m_alphaPoints4Sgm[points[1]] + M_PI - m_alphaPoints4Sgm[points[0]]);
2763 for (
unsigned short i = 0; i <= m_nAlphaPoints4Sgm - 2; ++i) {
2764 if (m_alphaPoints4Sgm[i] <= alphao && alphao < m_alphaPoints4Sgm[i + 1]) {
2767 weight = (alphao - m_alphaPoints4Sgm[points[0]]) / (m_alphaPoints4Sgm[points[1]] - m_alphaPoints4Sgm[points[0]]);
2775 void CDCGeometryPar::getClosestThetaPoints(
const double alpha,
const double theta,
double& weight,
unsigned short points[2])
const
2777 const double thetao = getOutgoingTheta(alpha, theta);
2779 if (thetao < m_thetaPoints[0]) {
2785 }
else if (m_thetaPoints[m_nThetaPoints - 1] <= thetao) {
2788 points[0] = m_nThetaPoints - 1;
2789 points[1] = m_nThetaPoints - 1;
2792 for (
unsigned short i = 0; i <= m_nThetaPoints - 2; ++i) {
2793 if (m_thetaPoints[i] <= thetao && thetao < m_thetaPoints[i + 1]) {
2796 weight = (thetao - m_thetaPoints[points[0]]) / (m_thetaPoints[points[1]] - m_thetaPoints[points[0]]);
2805 void CDCGeometryPar::getClosestThetaPoints4Sgm(
const double alpha,
const double theta,
double& weight,
2806 unsigned short points[2])
const
2808 const double thetao = getOutgoingTheta(alpha, theta);
2810 if (thetao < m_thetaPoints4Sgm[0]) {
2814 }
else if (m_thetaPoints4Sgm[m_nThetaPoints4Sgm - 1] <= thetao) {
2815 points[0] = m_nThetaPoints4Sgm - 1;
2816 points[1] = m_nThetaPoints4Sgm - 1;
2819 for (
unsigned short i = 0; i <= m_nThetaPoints4Sgm - 2; ++i) {
2820 if (m_thetaPoints4Sgm[i] <= thetao && thetao < m_thetaPoints4Sgm[i + 1]) {
2823 weight = (thetao - m_thetaPoints4Sgm[points[0]]) / (m_thetaPoints4Sgm[points[1]] - m_thetaPoints4Sgm[points[0]]);
2831 void CDCGeometryPar::setDisplacement()
2834 for (
const auto& disp : (*m_displacementFromDB)) {
2841 m_FWirPos[iLayer][iWire][0] += disp.getXFwd();
2842 m_FWirPos[iLayer][iWire][1] += disp.getYFwd();
2843 m_FWirPos[iLayer][iWire][2] += disp.getZFwd();
2844 m_BWirPos[iLayer][iWire][0] += disp.getXBwd();
2845 m_BWirPos[iLayer][iWire][1] += disp.getYBwd();
2846 m_BWirPos[iLayer][iWire][2] += disp.getZBwd();
2847 m_WireSagCoef[iLayer][iWire] = M_PI * m_senseWireDensity * m_senseWireDiameter * m_senseWireDiameter / (8.*
2848 (m_senseWireTension + disp.getTension()));
2854 void CDCGeometryPar::setShiftInSuperLayer()
2856 const unsigned short nLayers[nSuperLayers] = {8, 6, 6, 6, 6, 6, 6, 6, 6};
2858 for (
unsigned short SLayer = 0; SLayer < nSuperLayers; ++SLayer) {
2859 unsigned short firstCLayer = 0;
2860 for (
unsigned short i = 0; i < SLayer; ++i) {
2861 firstCLayer += nLayers[i];
2865 TVector3 firstBPos = wireBackwardPosition(firstCLayer, 0);
2866 for (
unsigned short Layer = 0; Layer < nLayers[SLayer]; ++Layer) {
2867 unsigned short CLayer = firstCLayer + Layer;
2869 if (CLayer == firstCLayer) {
2870 m_shiftInSuperLayer[SLayer][Layer] = 0;
2872 }
else if (CLayer == firstCLayer + 1) {
2873 TVector3 BPos = wireBackwardPosition(CLayer, 0);
2874 m_shiftInSuperLayer[SLayer][Layer] = (BPos.Cross(firstBPos)).Z() > 0. ? -1 : 1;
2878 if (Layer % 2 == 0) {
2879 m_shiftInSuperLayer[SLayer][Layer] = 0;
2881 m_shiftInSuperLayer[SLayer][Layer] = m_shiftInSuperLayer[SLayer][1];
2889 signed short CDCGeometryPar::getShiftInSuperLayer(
unsigned short iSuperLayer,
unsigned short iLayer)
const
2891 return m_shiftInSuperLayer[iSuperLayer][iLayer];