Belle II Software development
Modules | Namespaces | Classes | Typedefs | Functions
cdc
Packages

Modules

 cdc data objects
 
 cdc modules
 

Namespaces

namespace  Belle2::CDC
 

Classes

class  CDCDedx1DCellAlgorithm
 A calibration algorithm for CDC dE/dx electron: 1D enta cleanup correction. More...
 
class  CDCDedx2DCellAlgorithm
 A calibration algorithm for CDC dE/dx electron 2D enta vs doca correction. More...
 
class  CDCDedxBadWireAlgorithm
 A calibration algorithm for CDC dE/dx to find the bad wires. More...
 
class  CDCDedxCosEdgeAlgorithm
 A calibration algorithm for CDC dE/dx electron cos(theta) dependence. More...
 
class  CDCDedxCosineAlgorithm
 A calibration algorithm for CDC dE/dx electron cos(theta) dependence. More...
 
class  CDCDedxHadBGAlgorithm
 A calibration algorithm for beta gamma curve and resolution fitting and save payloads. More...
 
class  CDCDedxHadSatAlgorithm
 A calibration algorithm for hadron saturation. More...
 
class  CDCDedxInjectTimeAlgorithm
 A calibration algorithm for CDC dE/dx injection time (HER/LER) More...
 
class  CDCDedxMomentumAlgorithm
 A calibration algorithm for CDC dE/dx electron cos(theta) dependence. More...
 
class  CDCDedxRunGainAlgorithm
 A calibration algorithm for CDC dE/dx run gains. More...
 
class  CDCDedxWireGainAlgorithm
 A calibration algorithm for CDC dE/dx wire gains. More...
 
class  HadronBgPrep
 Class to prepare sample for fitting in beta gamma bins. More...
 
class  HadronCalibration
 Class to perform the fitting in beta gamma bins. More...
 
class  HadronPrep
 Class to prepare sample for hadron saturation calibration. More...
 
class  HadronSaturation
 Class to perform the hadron saturation calibration. More...
 
class  CDCDatabaseImporter
 CDC database importer. More...
 
class  SliceFit
 Class to do the slice fit. More...
 
class  CDCADCDeltaPedestals
 Database object for ADC pedestals. More...
 
class  CDCAlignment
 CDC alignment constants. More...
 
class  CDCBadWires
 Database object for bad wires. More...
 
class  CDCChannelMap
 Database object of CDC channel map. More...
 
class  CDCCorrToThresholds
 Database object for correcting a simple threshold model in MC. More...
 
struct  asicChannel
 record to be used to store ASIC info More...
 
struct  adcChannelPair
 pair ADC, channel More...
 
struct  adcAsicTuple
 tuple to store ADC,Channel -> 8 asicChannels More...
 
struct  adc_search
 functions to search in the sorted list of tuples More...
 
class  CDCCrossTalkLibrary
 Database object for ASIC crosstalk library. More...
 
class  CDCDedx1DCell
 dE/dx wire gain calibration constants More...
 
class  CDCDedx2DCell
 dE/dx wire gain calibration constants More...
 
class  CDCDedxADCNonLinearity
 dE/dx electronic ADC non-linearity correction for highly ionising particles (used in offline hadron saturation calibration) parameters are for X vs Y relation and sep for inner and outer layer vector array 0,1 for inner and 2,3 for outer layers More...
 
class  CDCDedxBadWires
 dE/dx wire gain calibration constants More...
 
class  CDCDedxCosineCor
 dE/dx wire gain calibration constants More...
 
class  CDCDedxCosineEdge
 dE/dx special large cosine calibration to fix bending shoulder at large costh More...
 
class  CDCDedxDatabaseImporter
 dE/dx database importer. More...
 
class  CDCDedxHadronCor
 dE/dx hadron saturation parameterization constants More...
 
class  CDCDedxInjectionTime
 dE/dx injection time calibration constants More...
 
class  CDCDedxMeanPars
 dE/dx mean (curve versus beta-gamma) parameterization constants More...
 
class  CDCDedxMomentumCor
 dE/dx wire gain calibration constants More...
 
class  CDCdEdxPDFs
 Specialized class for holding the CDC dE/dx PDFs. More...
 
class  CDCDedxRunGain
 dE/dx run gain calibration constants More...
 
class  CDCDedxScaleFactor
 dE/dx run gain calibration constants More...
 
class  CDCDedxSigmaPars
 dE/dx sigma (versus beta-gamma) parameterization constants More...
 
class  CDCDedxWireGain
 dE/dx wire gain calibration constants More...
 
class  CDCDisplacement
 Database object for displacement of sense wire position. More...
 
class  CDCEDepToADCConversions
 Database object for energy-deposit to ADC-count conversion. More...
 
class  CDCFEElectronics
 Database object for Fron-endt electronics params. More...
 
class  CDCFEEParams
 Database object for FEE params. More...
 
class  CDCFudgeFactorsForSigma
 Database object for fudge factors for CDC space resol. More...
 
class  CDCGeometry
 The Class for CDC geometry. More...
 
class  CDCLayerAlignment
 CDC layers alignment constants. More...
 
class  CDClayerTimeCut
 Database object for timing offset (t0). More...
 
class  CDCMisalignment
 CDC misalignment constants. More...
 
class  CDCPropSpeeds
 Database object for signal propagation speed along the wire. More...
 
class  CDCSpaceResols
 Database object for space resolutions. More...
 
class  CDCTimeWalks
 Database object for time-walk. More...
 
class  CDCTimeZeros
 Database object for timing offset (t0). More...
 
class  CDCTriggerPlane
 Database object for timing offset (t0). More...
 
class  CDCWireHitRequirements
 Database object containing cut values to filter CDCWireHits. More...
 
class  CDCXtRelations
 Database object for xt-relations. More...
 
class  CDCDedxHadSat
 Class to hold the hadron saturation functions. More...
 
class  CDCDedxMeanPred
 Class to hold the prediction of mean as a function of beta-gamma (bg) More...
 
class  CDCDedxSigmaPred
 Class to hold the prediction of resolution depending dE/dx, nhit, and cos(theta) More...
 
class  CDCDedxWidgetCurve
 Class to hold the beta-gamma (bg) mean function. More...
 
class  CDCDedxWidgetSigma
 Class to hold the beta-gamma (bg) resolution function. More...
 

Typedefs

typedef array< asicChannel, 8 > asicChannels
 fixed sized array of ASIC channels
 

Functions

CDCDedxTrack const * getDedxFromParticle (Particle const *particle)
 CDC dEdx value from particle.
 
double D2I (double cosTheta, double D) const
 hadron saturation parameterization part 2
 
double I2D (double cosTheta, double I) const
 hadron saturation parameterization part 1
 
double meanCurve (double x, const double *par, int version) const
 parameterized beta-gamma curve for predicted means
 
double getMean (double bg) const
 Returns the predicted dE/dx mean at given beta-gamma.
 
double sigmaCurve (double x, const double *par, int version) const
 parameterized resolution for predictions
 
double getSigma (double dedx, double nhit, double cosTheta, double timeReso) const
 Returns predicted dE/dx sigma.
 
 CDCSensitiveDetector (G4String name, G4double thresholdEnergyDeposit, G4double thresholdKineticEnergy)
 Constructor.
 
void Initialize (G4HCofThisEvent *) override
 Register CDC hits collection into G4HCofThisEvent.
 
bool step (G4Step *aStep, G4TouchableHistory *history) override
 Process each step and calculate variables defined in CDCB4VHit.
 
void EndOfEvent (G4HCofThisEvent *) override
 Do what you want to do at the beginning of each event (why this is not called ?)
 
void saveSimHit (const G4int layerId, const G4int wireId, const G4int trackID, const G4int pid, const G4double distance, const G4double tof, const G4double edep, const G4double stepLength, const G4ThreeVector &mom, const G4ThreeVector &posW, const G4ThreeVector &posIn, const G4ThreeVector &posOut, const G4ThreeVector &posTrack, const G4int lr, const G4int NewLrRaw, const G4int NewLr, const G4double speed, const G4double hitWeight)
 Save CDCSimHit into datastore.
 
void CellBound (const G4int layerId, const G4int ic1, const G4int ic2, const G4double venter[6], const G4double vexit[6], const G4double s1, const G4double s2, G4double xint[6], G4double &sint, G4int &iflag)
 Calculate intersection of track with cell boundary.
 
void GCUBS (const G4double x, const G4double y, const G4double d1, const G4double d2, G4double a[4])
 
void for_Rotat (const G4double bfld[3])
 Calculates a rotation matrix.
 
void Rotat (G4double &x, G4double &y, G4double &z, const int mode)
 Translation method.
 
void Rotat (G4double x[3], const int mode)
 Overloaded translation method.
 
void HELWIR (const G4double xwb4, const G4double ywb4, const G4double zwb4, const G4double xwf4, const G4double ywf4, const G4double zwf4, const G4double xp, const G4double yp, const G4double zp, const G4double px, const G4double py, const G4double pz, const G4double B_kG[3], const G4double charge, const G4int ntryMax, G4double &distance, G4double q2[3], G4double q1[3], G4double q3[3], G4int &ntry)
 Calculate closest points between helix and wire.
 
void Mvopr (const G4int ndim, const G4double b[3], const G4double m[3][3], const G4double a[3], G4double c[3], const G4int mode)
 Calculate the result of a matrix times vector.
 
std::vector< int > WireId_in_hit_order (int id0, int id1, int nWires)
 Sort wire id.
 
G4double ClosestApproach (G4ThreeVector bwp, G4ThreeVector fwp, G4ThreeVector posIn, G4ThreeVector posOut, G4ThreeVector &hitPosition, G4ThreeVector &wirePosition)
 Assume line track to calculate distance between track and wire (drift length).
 
void setModifiedLeftRightFlag ()
 set left/right flag modified for tracking
 
void reAssignLeftRightInfo ()
 Re-assign left/right info.
 
unsigned short areNeighbors (const WireID &wireId, const WireID &otherWireId) const
 Check if neighboring cell in the same super-layer; essentially a copy from cdcLocalTracking/mclookup.
 
unsigned short areNeighbors (unsigned short iCLayer, unsigned short iSuperLayer, unsigned short iLayer, unsigned short iWire, const WireID &otherWireId) const
 Check if neighboring cell in the same super-layer; essentially a copy from cdcLocalTracking/mclookup.
 

Detailed Description

Typedef Documentation

◆ asicChannels

typedef array<asicChannel, 8> asicChannels

fixed sized array of ASIC channels

Definition at line 28 of file CDCCrossTalkClasses.h.

Function Documentation

◆ areNeighbors() [1/2]

unsigned short areNeighbors ( const WireID &  wireId,
const WireID &  otherWireId 
) const
private

Check if neighboring cell in the same super-layer; essentially a copy from cdcLocalTracking/mclookup.

Parameters
[in]wireIdwire-id. in question (reference)
[in]otherWireIdanother wire-id. in question

Definition at line 1609 of file CDCSensitiveDetector.cc.

1610 {
1611 //require within the same super-layer
1612 if (otherWireId.getISuperLayer() != wireId.getISuperLayer()) return 0;
1613
1614 const signed short iWire = wireId.getIWire();
1615 const signed short iOtherWire = otherWireId.getIWire();
1616 const signed short iCLayer = wireId.getICLayer();
1617 const signed short iOtherCLayer = otherWireId.getICLayer();
1618
1619 //require nearby wire
1620 if (iWire == iOtherWire) {
1621 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1622 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1623 } else {
1624 return 0;
1625 }
1626 // std::cout <<"iCLayer,iLayer,nShifts= " << iCLayer <<" "<< iLayer <<" "<< nShifts(iCLayer) << std::endl;
1627
1628 signed short iLayerDifference = otherWireId.getILayer() - wireId.getILayer();
1629 if (abs(iLayerDifference) > 1) return 0;
1630
1631 if (iLayerDifference == 0) {
1632 if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer))) return CW_NEIGHBOR;
1633 else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) return CCW_NEIGHBOR;
1634 else return 0;
1635 } else if (iLayerDifference == -1) {
1636 // const CCWInfo deltaShift = otherLayer.getShift() - layer.getShift();
1637 const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1638 m_cdcgp->getShiftInSuperLayer(wireId.getISuperLayer(), wireId.getILayer());
1639 // std::cout <<"in deltaShift,iOtherWire,iWire= " << deltaShift <<" "<< iOtherWire <<" "<< iWire << std::endl;
1640 if (iWire == iOtherWire) {
1641 if (deltaShift == CW) return CW_IN_NEIGHBOR;
1642 else if (deltaShift == CCW) return CCW_IN_NEIGHBOR;
1643 else return 0;
1644 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1645 if (deltaShift == CCW) return CW_IN_NEIGHBOR;
1646 else return 0;
1647 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1648 if (deltaShift == CW) return CCW_IN_NEIGHBOR;
1649 else return 0;
1650 } else return 0;
1651 } else if (iLayerDifference == 1) {
1652 // const CCWInfo deltaShift = otherLayer.getShift() - layer.getShift();
1653 const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1654 m_cdcgp->getShiftInSuperLayer(wireId.getISuperLayer(), wireId.getILayer());
1655 // std::cout <<"out deltaShift,iOtherWire,iWire= " << deltaShift <<" "<< iOtherWire <<" "<< iWire << std::endl;
1656 if (iWire == iOtherWire) {
1657 if (deltaShift == CW) return CW_OUT_NEIGHBOR;
1658 else if (deltaShift == CCW) return CCW_OUT_NEIGHBOR;
1659 else return 0;
1660 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1661 if (deltaShift == CCW) return CW_OUT_NEIGHBOR;
1662 else return 0;
1663 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1664 if (deltaShift == CW) return CCW_OUT_NEIGHBOR;
1665 else return 0;
1666 } else return 0;
1667 } else return 0;
1668
1669 }
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1254
Belle2::CDC::CDCGeometryPar::getShiftInSuperLayer
signed short getShiftInSuperLayer(unsigned short iSuperLayer, unsigned short iLayer) const
Returns shift in the super-layer.
Definition: CDCGeometryPar.cc:3010
Belle2::CDC::CDCSensitiveDetector::CW_NEIGHBOR
const signed short CW_NEIGHBOR
Constant for clockwise.
Definition: CDCSensitiveDetector.h:293
Belle2::CDC::CDCSensitiveDetector::CCW_NEIGHBOR
const signed short CCW_NEIGHBOR
Constant for counterclockwise.
Definition: CDCSensitiveDetector.h:296
Belle2::CDC::CDCSensitiveDetector::m_cdcgp
CDCGeometryPar * m_cdcgp
Pointer to CDCGeometryPar object.
Definition: CDCSensitiveDetector.h:259
Belle2::CDC::CDCSensitiveDetector::CW_IN_NEIGHBOR
const signed short CW_IN_NEIGHBOR
Constant for clockwise inwards.
Definition: CDCSensitiveDetector.h:294
Belle2::CDC::CDCSensitiveDetector::CCW_OUT_NEIGHBOR
const signed short CCW_OUT_NEIGHBOR
Constant for counterclockwise outwards.
Definition: CDCSensitiveDetector.h:297
Belle2::CDC::CDCSensitiveDetector::CW_OUT_NEIGHBOR
const signed short CW_OUT_NEIGHBOR
Constant for clockwise outwards.
Definition: CDCSensitiveDetector.h:292
Belle2::CDC::CDCSensitiveDetector::CCW
const signed short CCW
Constant for counterclockwise orientation.
Definition: CDCSensitiveDetector.h:290
Belle2::CDC::CDCSensitiveDetector::CCW_IN_NEIGHBOR
const signed short CCW_IN_NEIGHBOR
Constant for counterclockwise inwards.
Definition: CDCSensitiveDetector.h:295
Belle2::CDC::CDCSensitiveDetector::CW
const signed short CW
Constant for clockwise orientation.
Definition: CDCSensitiveDetector.h:291

◆ areNeighbors() [2/2]

unsigned short areNeighbors ( unsigned short  iCLayer,
unsigned short  iSuperLayer,
unsigned short  iLayer,
unsigned short  iWire,
const WireID &  otherWireId 
) const
private

Check if neighboring cell in the same super-layer; essentially a copy from cdcLocalTracking/mclookup.

Parameters
[in]iCLayerlater-id (continuous) in question (reference)
[in]iSuperLayersuper-later-id in question (reference)
[in]iLayerlater-id in the super-layer in question (reference)
[in]iWirewire-id in the layer in question (reference)
[in]otherWireIdanother wire-id. in question

Definition at line 1671 of file CDCSensitiveDetector.cc.

1673 {
1674 //require within the same super-layer
1675 if (otherWireId.getISuperLayer() != iSuperLayer) return 0;
1676
1677 const signed short iOtherWire = otherWireId.getIWire();
1678 const signed short iOtherCLayer = otherWireId.getICLayer();
1679
1680 //require nearby wire
1681 if (iWire == iOtherWire) {
1682 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1683 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1684 } else {
1685 return 0;
1686 }
1687
1688 // std::cout <<"iCLayer,iLayer,nShifts= " << iCLayer <<" "<< iLayer <<" "<< nShifts(iCLayer) << std::endl;
1689 signed short iLayerDifference = otherWireId.getILayer() - iLayer;
1690 if (abs(iLayerDifference) > 1) return 0;
1691
1692 if (iLayerDifference == 0) {
1693 if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer))) return CW_NEIGHBOR;
1694 else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) return CCW_NEIGHBOR;
1695 else return 0;
1696 } else if (iLayerDifference == -1) {
1697 // const CCWInfo deltaShift = otherLayer.getShift() - layer.getShift();
1698 const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1699 m_cdcgp->getShiftInSuperLayer(iSuperLayer, iLayer);
1700 // std::cout <<"in deltaShift,iOtherWire,iWire= " << deltaShift <<" "<< iOtherWire <<" "<< iWire << std::endl;
1701 if (iWire == iOtherWire) {
1702 if (deltaShift == CW) return CW_IN_NEIGHBOR;
1703 else if (deltaShift == CCW) return CCW_IN_NEIGHBOR;
1704 else return 0;
1705 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1706 if (deltaShift == CCW) return CW_IN_NEIGHBOR;
1707 else return 0;
1708 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1709 if (deltaShift == CW) return CCW_IN_NEIGHBOR;
1710 else return 0;
1711 } else return 0;
1712 } else if (iLayerDifference == 1) {
1713 // const CCWInfo deltaShift = otherLayer.getShift() - layer.getShift();
1714 const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1715 m_cdcgp->getShiftInSuperLayer(iSuperLayer, iLayer);
1716 // std::cout <<"out deltaShift,iOtherWire,iWire= " << deltaShift <<" "<< iOtherWire <<" "<< iWire << std::endl;
1717 if (iWire == iOtherWire) {
1718 if (deltaShift == CW) return CW_OUT_NEIGHBOR;
1719 else if (deltaShift == CCW) return CCW_OUT_NEIGHBOR;
1720 else return 0;
1721 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1722 if (deltaShift == CCW) return CW_OUT_NEIGHBOR;
1723 else return 0;
1724 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1725 if (deltaShift == CW) return CCW_OUT_NEIGHBOR;
1726 else return 0;
1727 } else return 0;
1728 } else return 0;
1729
1730 }

◆ CDCSensitiveDetector()

CDCSensitiveDetector ( G4String  name,
G4double  thresholdEnergyDeposit,
G4double  thresholdKineticEnergy 
)

Constructor.

Definition at line 48 of file CDCSensitiveDetector.cc.

48 :
49 SensitiveDetectorBase(name, Const::CDC),
50 // m_cdcgp(CDCGeometryPar::Instance()),
51 m_cdcgp(nullptr),
52 m_thresholdEnergyDeposit(thresholdEnergyDeposit),
53 m_thresholdKineticEnergy(thresholdKineticEnergy), m_hitNumber(0)
54 {
55 RelationArray cdcSimHitRel(m_MCParticles, m_CDCSimHits);
56 registerMCParticleRelation(cdcSimHitRel);
57 // registerMCParticleRelation(cdcSimHitRel, RelationArray::c_doNothing);
58 // registerMCParticleRelation(cdcSimHitRel, RelationArray::c_negativeWeight);
59 // registerMCParticleRelation(cdcSimHitRel, RelationArray::c_deleteElement);
60 // registerMCParticleRelation(cdcSimHitRel, RelationArray::c_zeroWeight);
61 m_CDCSimHits.registerInDataStore();
62 m_MCParticles.registerRelationTo(m_CDCSimHits);
63
64 const CDCSimControlPar& cntlp = CDCSimControlPar::getInstance();
65
66 m_thresholdEnergyDeposit = cntlp.getThresholdEnergyDeposit();
67 m_thresholdEnergyDeposit *= CLHEP::GeV; //GeV to MeV (=unit in G4)
68 B2DEBUG(150, "CDCSensitiveDetector: Threshold energy (MeV): " << m_thresholdEnergyDeposit);
69 m_thresholdKineticEnergy = 0.0; // Dummy to avoid a warning (tentative).
70
71 m_wireSag = cntlp.getWireSag();
72 // m_wireSag = false;
73 B2DEBUG(150, "CDCSensitiveDetector: Sense wire sag on(=1)/off(=0): " << m_wireSag);
74
75 m_modifiedLeftRightFlag = cntlp.getModLeftRightFlag();
76 B2DEBUG(150, "CDCSensitiveDetector: Set left/right flag modified for tracking (=1)/ not set (=0): " << m_modifiedLeftRightFlag);
77
78 m_minTrackLength = cntlp.getMinTrackLength();
79 m_minTrackLength *= CLHEP::cm; //cm to mm (=unit in G4)
80 B2DEBUG(150, "CDCSensitiveDetector: MinTrackLength (mm): " << m_minTrackLength);
81
82 //For activating Initialize and EndOfEvent functions
83 //but not work --> commented out for a while...
84 // if (m_modifiedLeftRightFlag) {
85 // G4SDManager::GetSDMpointer()->AddNewDetector(this);
86 // }
87 }
Belle2::CDC::CDCSensitiveDetector::m_thresholdEnergyDeposit
G4double m_thresholdEnergyDeposit
Threshold energy deposit to be stored.
Definition: CDCSensitiveDetector.h:264
Belle2::CDC::CDCSensitiveDetector::m_wireSag
G4bool m_wireSag
Switch to activate wire sag effect.
Definition: CDCSensitiveDetector.h:271
Belle2::CDC::CDCSensitiveDetector::m_CDCSimHits
StoreArray< CDCSimHit > m_CDCSimHits
CDC simulation hits.
Definition: CDCSensitiveDetector.h:243
Belle2::CDC::CDCSensitiveDetector::m_thresholdKineticEnergy
G4double m_thresholdKineticEnergy
Threshold kinetic energy to be stored.
Definition: CDCSensitiveDetector.h:269
Belle2::CDC::CDCSensitiveDetector::m_modifiedLeftRightFlag
G4bool m_modifiedLeftRightFlag
Switch for left/right flag modified for tracking.
Definition: CDCSensitiveDetector.h:273
Belle2::CDC::CDCSensitiveDetector::m_hitNumber
int m_hitNumber
The current number of created hits in an event.
Definition: CDCSensitiveDetector.h:277
Belle2::CDC::CDCSensitiveDetector::m_MCParticles
StoreArray< MCParticle > m_MCParticles
MC particles.
Definition: CDCSensitiveDetector.h:240
Belle2::CDC::CDCSimControlPar::getInstance
static CDCSimControlPar & getInstance()
Static method to get a reference to the CDCSimControlPar instance.
Definition: CDCSimControlPar.cc:16
Belle2::Simulation::SensitiveDetectorBase::registerMCParticleRelation
static void registerMCParticleRelation(const std::string &name, RelationArray::EConsolidationAction ignoreAction=RelationArray::c_negativeWeight)
Register an relation involving MCParticles.
Definition: SensitiveDetectorBase.cc:22
Belle2::Simulation::SensitiveDetectorBase::SensitiveDetectorBase
SensitiveDetectorBase(const std::string &name, Const::EDetector subdetector)
Create a new Sensitive detecor with a given name and belonging to a given subdetector.
Definition: SensitiveDetectorBase.h:35
Belle2::StoreArray::registerRelationTo
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:140

◆ CellBound()

void CellBound ( const G4int  layerId,
const G4int  ic1,
const G4int  ic2,
const G4double  venter[6],
const G4double  vexit[6],
const G4double  s1,
const G4double  s2,
G4double  xint[6],
G4double &  sint,
G4int &  iflag 
)
private

Calculate intersection of track with cell boundary.

Parameters
[in]layerIdId of the layer.
[in]ic1serial cell number (start w/ one) of entrance.
[in]ic2serial cell number (start w/ one) of exit.
[in]venter(x,y,z,px/p,py/p,pz/p) at entrance.
[in]vexit(x,y,z,px/p,py/p,pz/p) at exit.
[in]s1track length at entrance.
[in]s2track length at exit.
[out]xint(x,y,z,px/p,py/p,pz/p) at intersection of cell boundary.
[out]sinttrack length at intersection of cell boundary.
[out]iflagreturn code from GIPLAN.

Definition at line 689 of file CDCSensitiveDetector.cc.

697 {
698 //---------------------------------------------------------------------------
699 // (Purpose)
700 // calculate an intersection of track with cell boundary.
701 //
702 // (Relations)
703 // Calls GCUBS
704 //
705 // (Arguments)
706 // input
707 // ic1 serial cell# (start w/ one) of entrance.
708 // ic2 serial cell# (start w/ one) of exit.
709 // venter(6) (x,y,z,px/p,py/p,pz/p) at entrance.
710 // vexit(6) (x,y,z,px/p,py/p,pz/p) at exit.
711 // s1 track length at entrance.
712 // s2 track length at exit.
713 // output
714 // xint(6) (x,y,z,px/p,py/p,pz/p) at intersection of cell boundary.
715 // sint track length at intersection of cell boundary.
716 // iflag return code.
717 //
718 // N.B.(TODO ?) CDC misalignment wrt Belle2 coordinate system is ignored
719 // when calculating the cell-boundary assuming misalign. is small.
720 //--------------------------------------------------------------------------
721
722 G4double div = m_cdcgp->nWiresInLayer(layerId);
723
724 //Check if s1, s2, ic1 and ic2 are ok
725 if (s1 >= s2) {
726 B2ERROR("CDCSensitiveDetector: s1(=" << s1 << ") > s2(=" << s2 << ")");
727 }
728 if (std::abs(ic1 - ic2) != 1) {
729 if (ic1 == 1 && ic2 == div) {
730 } else if (ic1 == div && ic2 == 1) {
731 } else {
732 B2ERROR("CDCSensitiveDetector: |ic1 - ic2| != 1 in CellBound; " << "ic1=" << ic1 << " " << "ic2=" << ic2);
733 }
734 }
735
736 //get wire positions for the entrance cell
737 G4double xwb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).x();
738 G4double ywb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).y();
739 G4double zwb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).z();
740 G4double xwf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).x();
741 G4double ywf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).y();
742 G4double zwf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).z();
743
744 /*
745 G4double pathl = sqrt((vexit[0] - venter[0]) * (vexit[0] - venter[0])
746 + (vexit[1] - venter[1]) * (vexit[1] - venter[1])
747 + (vexit[2] - venter[2]) * (vexit[2] - venter[2]));
748 std::cout << "app pathl= " << pathl << std::endl;
749 G4double dot = venter[3] * vexit[3] + venter[4] * vexit[4];
750 dot /= sqrt(venter[3] * venter[3] + venter[4] * venter[4]);
751 dot /= sqrt( vexit[3] * vexit[3] + vexit[4] * vexit[4]);
752 if (dot < 0.) std::cout <<"negativedot= " << dot << std::endl;
753 */
754
755 //copy arrays
756 G4double xx1[6], xx2[6];
757 for (int i = 0; i < 6; ++i) {
758 xx1[i] = venter[i];
759 xx2[i] = vexit [i];
760 }
761
762 //calculate the field wire position betw. cell#1 and #2
763 G4double psi = double(ic2 - ic1) * CLHEP::pi / div;
764 if (ic1 == 1 && ic2 == div) {
765 psi = -CLHEP::pi / div;
766 } else if (ic1 == div && ic2 == 1) {
767 psi = CLHEP::pi / div;
768 }
769 G4double cospsi = cos(psi);
770 G4double sinpsi = sin(psi);
771
772 G4double xfwb = cospsi * xwb - sinpsi * ywb;
773 G4double yfwb = sinpsi * xwb + cospsi * ywb;
774 G4double xfwf = cospsi * xwf - sinpsi * ywf;
775 G4double yfwf = sinpsi * xwf + cospsi * ywf;
776 G4double zfwb = zwb;
777 G4double zfwf = zwf;
778
779 //prepare quantities related to the cell-boundary
780 G4double vx = xfwf - xfwb;
781 G4double vy = yfwf - yfwb;
782 G4double vz = zfwf - zfwb;
783 G4double vv = sqrt(vx * vx + vy * vy + vz * vz);
784 vx /= vv; vy /= vv; vz /= vv;
785
786 //translate to make the cubic description easier
787 G4double shiftx = (xx1[0] + xx2[0]) * 0.5;
788 G4double shifty = (xx1[1] + xx2[1]) * 0.5;
789 G4double shiftz = (xx1[2] + xx2[2]) * 0.5;
790 G4double shifts = (s1 + s2) * 0.5;
791 G4double xshft = xx1[0] - shiftx;
792 G4double yshft = xx1[1] - shifty;
793 G4double zshft = xx1[2] - shiftz;
794 G4double sshft = s1 - shifts;
795
796 //approximate the trajectroy by cubic curves
797 G4double pabs1 = sqrt(xx1[3] * xx1[3] + xx1[4] * xx1[4] + xx1[5] * xx1[5]);
798 G4double pabs2 = sqrt(xx2[3] * xx2[3] + xx2[4] * xx2[4] + xx2[5] * xx2[5]);
799 // std::cout << "pabs1,2= " << pabs1 <<" "<< pabs2 << std::endl;
800
801 G4double a[4] = {0.}, b[4] = {0.}, c[4] = {0.};
802
803 if (m_magneticField) {
804 GCUBS(sshft, xshft, xx1[3] / pabs1, xx2[3] / pabs2, a);
805 GCUBS(sshft, yshft, xx1[4] / pabs1, xx2[4] / pabs2, b);
806 GCUBS(sshft, zshft, xx1[5] / pabs1, xx2[5] / pabs2, c);
807 // std::cout <<"a= " << a[0] <<" "<< a[1] <<" "<< a[2] <<" "<< a[3] << std::endl;
808 // std::cout <<"b= " << b[0] <<" "<< b[1] <<" "<< b[2] <<" "<< b[3] << std::endl;
809 // std::cout <<"c= " << c[0] <<" "<< c[1] <<" "<< c[2] <<" "<< c[3] << std::endl;
810 } else {
811 //n.b. following is really better ?
812 a[1] = xshft / sshft;
813 b[1] = yshft / sshft;
814 c[1] = zshft / sshft;
815 }
816
817 //calculate an int. point betw. the trajectory and the cell-boundary
818 G4double stry(0.), xtry(0.), ytry(0.), ztry(0.);
819 G4double beta(0.), xfw(0.), yfw(0.);
820 G4double sphi(0.), cphi(0.), dphil(0.), dphih(0.);
821 const G4int maxTrials = 100;
822 const G4double eps = 5.e-4;
823 G4double sl = sshft; // negative value
824 G4double sh = -sshft; // positive value
825 G4int i = 0;
826
827 //set initial value (dphil) for the 1st iteration
828 stry = sl;
829 xtry = shiftx + a[0] + stry * (a[1] + stry * (a[2] + stry * a[3]));
830 ytry = shifty + b[0] + stry * (b[1] + stry * (b[2] + stry * b[3]));
831 ztry = shiftz + c[0] + stry * (c[1] + stry * (c[2] + stry * c[3]));
832 beta = (ztry - zfwb) / vz;
833 xfw = xfwb + beta * vx;
834 yfw = yfwb + beta * vy;
835 sphi = (xtry * yfw - ytry * xfw);
836 cphi = (xtry * xfw + ytry * yfw);
837 dphil = atan2(sphi, cphi); //n.b. no need to conv. to dphi...
838
839 iflag = 1;
840
841 while (((sh - sl) > eps) && (i < maxTrials)) {
842 stry = 0.5 * (sl + sh);
843 xtry = shiftx + a[0] + stry * (a[1] + stry * (a[2] + stry * a[3]));
844 ytry = shifty + b[0] + stry * (b[1] + stry * (b[2] + stry * b[3]));
845 ztry = shiftz + c[0] + stry * (c[1] + stry * (c[2] + stry * c[3]));
846 beta = (ztry - zfwb) / vz;
847 xfw = xfwb + beta * vx;
848 yfw = yfwb + beta * vy;
849
850 sphi = (xtry * yfw - ytry * xfw);
851 cphi = (xtry * xfw + ytry * yfw);
852 dphih = atan2(sphi, cphi); //n.b. no need to conv. to dphi...
853
854 if (dphil * dphih > 0.) {
855 sl = stry;
856 } else {
857 sh = stry;
858 }
859 ++i;
860 }
861
862 // std::cout << "itry= " << i << std::endl;
863 if (i >= maxTrials - 1) {
864 iflag = 0;
865 B2WARNING("CDCSensitiveDetector: No intersection ?");
866
867 /* G4double ds = 1.e-4;
868 G4int imax = (s2 - s1) / ds + 1;
869 G4double rdphimin = DBL_MAX;
870
871 for (i=0; i <= imax; ++i) {
872 stry = sshft + i * ds;
873 xtry = shiftx + a[0] + stry * (a[1] + stry * (a[2] + stry * a[3]));
874 ytry = shifty + b[0] + stry * (b[1] + stry * (b[2] + stry * b[3]));
875 ztry = shiftz + c[0] + stry * (c[1] + stry * (c[2] + stry * c[3]));
876 beta = (ztry - zfwb) / vz;
877 xfw = xfwb + beta * vx;
878 yfw = yfwb + beta * vy;
879
880 sphi = (xtry * yfw - ytry * xfw);
881 cphi = (xtry * xfw + ytry * yfw);
882 dphi = atan2(sphi, cphi);
883 rdphi = sqrt(xfw * xfw + yfw * yfw) * dphi;
884
885 if ( std::abs(rdphi) < std::abs(rdphimin)) {
886 rdphimin = rdphi;
887 imin = i;
888 }
889 }
890 */
891 }
892 // sint = sshft + imin * ds;
893 sint = stry;
894
895 // std::cout <<"i,dphil,dphih,sint= " << i <<" "<< dphil <<" "<< dphih <<" "<< sint << std::endl;
896 //get the trajectory at the int. point
897 xint[0] = a[0] + sint * (a[1] + sint * (a[2] + sint * a[3]));
898 xint[1] = b[0] + sint * (b[1] + sint * (b[2] + sint * b[3]));
899 xint[2] = c[0] + sint * (c[1] + sint * (c[2] + sint * c[3]));
900 xint[3] = a[1] + sint * (2. * a[2] + 3. * sint * a[3]);
901 xint[4] = b[1] + sint * (2. * b[2] + 3. * sint * b[3]);
902 xint[5] = c[1] + sint * (2. * c[2] + 3. * sint * c[3]);
903
904 //translate back to the lab. frame
905 xint[0] += shiftx;
906 xint[1] += shifty;
907 xint[2] += shiftz;
908 sint += shifts;
909 /*
910 std::cout <<"s1,s2,sint= " << s1 <<" "<< s2 <<" "<< sint << std::endl;
911 std::cout <<" xx1= " << xx1[0] <<" "<< xx1[1] <<" "<< xx1[2] << std::endl;
912 std::cout <<" xx2= " << xx2[0] <<" "<< xx2[1] <<" "<< xx2[2] << std::endl;
913 std::cout <<"xint= " << xint[0] <<" "<< xint[1] <<" "<< xint[2] << std::endl;
914 */
915
916 /* if (((xx1[0] <= xint[0] && xint[0] <= xx2[0]) ||
917 (xx2[0] <= xint[0] && xint[0] <= xx1[0])) &&
918 ((xx1[1] <= xint[1] && xint[1] <= xx2[1]) ||
919 (xx2[1] <= xint[1] && xint[1] <= xx1[1])) &&
920 ((xx1[2] <= xint[2] && xint[2] <= xx2[2]) ||
921 (xx2[2] <= xint[2] && xint[2] <= xx1[2])) &&
922 (s1 <= sint && sint <= s2)) {
923 } else {
924 std::cout << "strangeinttersection" << std::endl;
925 }
926 */
927 //re-normalize to one since abs=1 is not guearanteed in the cubic approx.
928 G4double p = sqrt(xint[3] * xint[3] + xint[4] * xint[4] + xint[5] * xint[5]);
929 xint[3] /= p; xint[4] /= p; xint[5] /= p;
930 // std::cout << "norm= " << p << std::endl;
931 // std::cout <<"s1,s2,sint= " << s1 <<" "<< s2 <<" "<< sint << std::endl;
932 // std::cout <<"xint= " << xint[0] <<" "<< xint[1] <<" "<< xint[2] << std::endl;
933 // std::cout <<"xint= " << xint[3] <<" "<< xint[4] <<" "<< xint[5] << std::endl;
934 }
Belle2::CDC::CDCGeometryPar::wireForwardPosition
const B2Vector3D wireForwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the forward position of the input sense wire.
Definition: CDCGeometryPar.cc:1699
Belle2::CDC::CDCGeometryPar::wireBackwardPosition
const B2Vector3D wireBackwardPosition(uint layerId, int cellId, EWirePosition set=c_Base) const
Returns the backward position of the input sense wire.
Definition: CDCGeometryPar.cc:1748
Belle2::CDC::CDCSensitiveDetector::m_magneticField
G4int m_magneticField
Magnetic field is on or off.
Definition: CDCSensitiveDetector.h:248
Belle2::CDC::CDCSensitiveDetector::GCUBS
void GCUBS(const G4double x, const G4double y, const G4double d1, const G4double d2, G4double a[4])
Definition: CDCSensitiveDetector.cc:936
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28

◆ ClosestApproach()

G4double ClosestApproach ( G4ThreeVector  bwp,
G4ThreeVector  fwp,
G4ThreeVector  posIn,
G4ThreeVector  posOut,
G4ThreeVector &  hitPosition,
G4ThreeVector &  wirePosition 
)
private

Assume line track to calculate distance between track and wire (drift length).

Definition at line 1412 of file CDCSensitiveDetector.cc.

1414 {
1415
1416 B2Vector3D tbwp(bwp.x(), bwp.y(), bwp.z());
1417 B2Vector3D tfwp(fwp.x(), fwp.y(), fwp.z());
1418 B2Vector3D tposIn(posIn.x(), posIn.y(), posIn.z());
1419 B2Vector3D tposOut(posOut.x(), posOut.y(), posOut.z());
1420 B2Vector3D thitPosition(0., 0., 0.);
1421 B2Vector3D twirePosition(0., 0., 0.);
1422
1423 // G4double distance = m_cdcgp.ClosestApproach(tbwp, tfwp, tposIn, tposOut, thitPosition, twirePosition);
1424 G4double distance = CDC::ClosestApproach(tbwp, tfwp, tposIn, tposOut, thitPosition, twirePosition);
1425
1426 hitPosition.setX(thitPosition.x());
1427 hitPosition.setY(thitPosition.y());
1428 hitPosition.setZ(thitPosition.z());
1429
1430 wirePosition.setX(twirePosition.x());
1431 wirePosition.setY(twirePosition.y());
1432 wirePosition.setZ(twirePosition.z());
1433
1434 return distance;
1435 }
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Belle2::CDC::ClosestApproach
double ClosestApproach(const B2Vector3D &bwp, const B2Vector3D &fwp, const B2Vector3D &posIn, const B2Vector3D &posOut, B2Vector3D &hitPosition, B2Vector3D &wirePosition)
Returns a closest distance between a track and a wire.
Definition: ClosestApproach.cc:17

◆ D2I()

double D2I ( double  cosTheta,
double  D 
) const
private

hadron saturation parameterization part 2

Definition at line 19 of file CDCDedxHadronCor.cc.

20 {
21 const auto& params = m_hadronpars;
22 if (params.size() < 5) {
23 B2WARNING("Vector of dE/dx hadron constants too short!");
24 return D;
25 }
26
27 double projection = std::pow(std::abs(cosTheta), params[3]) + params[2];
28 if (projection == 0) {
29 B2WARNING("Something wrong with dE/dx hadron constants!");
30 return D;
31 }
32
33 double chargeDensity = D / projection;
34 double numerator = 1 + params[0] * chargeDensity;
35 double denominator = 1 + params[1] * chargeDensity;
36
37 if (denominator == 0) {
38 B2WARNING("Something wrong with dE/dx hadron constants!");
39 return D;
40 }
41
42 double I = D * params[4] * numerator / denominator;
43 return I;
44 }
Belle2::CDCDedxHadronCor::m_hadronpars
std::vector< double > m_hadronpars
dE/dx hadron correction parameters
Definition: CDCDedxHadronCor.h:88

◆ EndOfEvent()

void EndOfEvent ( G4HCofThisEvent *  )
override

Do what you want to do at the beginning of each event (why this is not called ?)

Do what you want to do at the end of each event

Definition at line 580 of file CDCSensitiveDetector.cc.

581 {
582 setModifiedLeftRightFlag();
583 }
Belle2::CDC::CDCSensitiveDetector::setModifiedLeftRightFlag
void setModifiedLeftRightFlag()
set left/right flag modified for tracking
Definition: CDCSensitiveDetector.cc:1439

◆ for_Rotat()

void for_Rotat ( const G4double  bfld[3])
private

Calculates a rotation matrix.

Calculates a rotation matrix. in advance at a local position in lab. frame. The rotation is done about the coord. origin; lab.-frame to B-field frame in which only Bz-comp. is non-zero.

Definition at line 971 of file CDCSensitiveDetector.cc.

972 {
973 //Calculates a rotation matrix in advance at a local position in lab.
974 //The rotation is done about the coord. origin; lab.-frame to B-field
975 //frame in which only Bz-comp. is non-zero.
976 //~dead copy of gsim_cdc_for_rotat.F in gsim-cdc for Belle (for tentative use)
977
978 if (m_nonUniformField == 0) return;
979
980 G4double bx, by, bz;
981 bx = bfld[0];
982 by = bfld[1];
983 bz = bfld[2];
984
985 //cal. rotation matrix
986 G4double bxz, bfield;
987 bxz = bx * bx + bz * bz;
988 bfield = bxz + by * by;
989 bxz = sqrt(bxz);
990 bfield = sqrt(bfield);
991
992 m_brot[0][0] = bz / bxz;
993 m_brot[1][0] = 0.;
994 m_brot[2][0] = -bx / bxz;
995 m_brot[0][1] = -by * bx / bxz / bfield;
996 m_brot[1][1] = bxz / bfield;
997 m_brot[2][1] = -by * bz / bxz / bfield;
998 m_brot[0][2] = bx / bfield;
999 m_brot[1][2] = by / bfield;
1000 m_brot[2][2] = bz / bfield;
1001
1002 return;
1003
1004 }
Belle2::CDC::CDCSensitiveDetector::m_brot
G4double m_brot[3][3]
a rotation matrix.
Definition: CDCSensitiveDetector.h:256
Belle2::CDC::CDCSensitiveDetector::m_nonUniformField
G4int m_nonUniformField
Magnetic field is uniform or non-uniform.
Definition: CDCSensitiveDetector.h:254

◆ GCUBS()

void GCUBS ( const G4double  x,
const G4double  y,
const G4double  d1,
const G4double  d2,
G4double  a[4] 
)
private 
                                                        *

Calculates a cubic through P1,(-X,Y1),(X,Y2),P2 * where Y2=-Y1 * Y=A(1)+A(2)*X+A(3)*X**2+A(4)*X**3 * The coordinate system is assumed to be the cms system * of P1,P2. *

  • ==>Called by : GIPLAN,GICYL * Author H.Boerner ********* *

Definition at line 936 of file CDCSensitiveDetector.cc.

937 {
938 //Original: GCUBS in Geant3
939 // ******************************************************************
940 // * *
941 // * Calculates a cubic through P1,(X,Y),(-X,-Y),P2 *
942 // * Y=A(1)+A(2)*X+A(3)*X**2+A(4)*X**3 *
943 // * The coordinate system is assumed to be the cms system *
944 // * of P1,P2. *
945 // * d1(2): directional cosine at P1(2). *
946 // * *
947 // * ==>Called by : GIPLAN,GICYL *
948 // * Author H.Boerner ********* *
949 // * *
950 // ******************************************************************
951
952 G4double fact(0);
953
954 if (x == 0.) goto L10;
955
956 fact = (d1 - d2) * 0.25;
957 a[0] = - 1. * fact * x;
958 a[2] = fact / x;
959 a[1] = (6. * y - (d1 + d2) * x) / (4. * x);
960 a[3] = ((d1 + d2) * x - 2.*y) / (4.*x * x * x);
961 return;
962
963L10:
964 a[0] = 0.;
965 a[1] = 1.;
966 a[2] = 0.;
967 a[3] = 0.;
968 }

◆ getDedxFromParticle()

CDCDedxTrack const * getDedxFromParticle ( Particle const *  particle)

CDC dEdx value from particle.

Definition at line 36 of file DedxVariables.cc.

37 {
38 const Track* track = particle->getTrack();
39 if (!track) {
40 return nullptr;
41 }
42
43 const CDCDedxTrack* dedxTrack = track->getRelatedTo<CDCDedxTrack>();
44 if (!dedxTrack) {
45 return nullptr;
46 }
47
48 return dedxTrack;
49 }

◆ getMean()

double getMean ( double  bg) const

Returns the predicted dE/dx mean at given beta-gamma.

Parameters
bgbeta-gamma
Returns
predicted dE/dx mean

Definition at line 38 of file CDCDedxMeanPars.cc.

39 {
40 // define the section of the mean to use
41 double A = 0, B = 0, C = 0;
42 if (bg < 4.5)
43 A = 1;
44 else if (bg < 10)
45 B = 1;
46 else
47 C = 1;
48
49 double parsA[9];
50 double parsB[5];
51 double parsC[5];
52
53 const auto& params = m_meanpars;
54 if (params.size() < 15) B2FATAL("CDCDedxMeanPars: vector of parameters too short");
55
56 parsA[0] = 1; parsB[0] = 2; parsC[0] = 3;
57 for (int i = 0; i < 15; ++i) {
58 if (i < 7) parsA[i + 1] = params[i];
59 else if (i < 11) parsB[i % 7 + 1] = params[i];
60 else parsC[i % 11 + 1] = params[i];
61 }
62
63 // calculate dE/dx from the Bethe-Bloch mean
64 double partA = meanCurve(bg, parsA, 0);
65 double partB = meanCurve(bg, parsB, 0);
66 double partC = meanCurve(bg, parsC, 0);
67
68 return (A * partA + B * partB + C * partC);
69 }
Belle2::CDCDedxMeanPars::m_meanpars
std::vector< double > m_meanpars
dE/dx mean parameters
Definition: CDCDedxMeanPars.h:84
Belle2::CDCDedxMeanPars::meanCurve
double meanCurve(double x, const double *par, int version) const
parameterized beta-gamma curve for predicted means
Definition: CDCDedxMeanPars.cc:19

◆ getSigma()

double getSigma ( double  dedx,
double  nhit,
double  cosTheta,
double  timeReso 
) const

Returns predicted dE/dx sigma.

Parameters
dedxpredicted mean
nhitnumber of measurements
cosThetacosine of polar angle
timeResoinjection time resolution scaling factor
Returns
predicted sigma

Definition at line 40 of file CDCDedxSigmaPars.cc.

41 {
42 const auto& params = m_sigmapars;
43 if (params.size() < 17) B2FATAL("CDCDedxSigmaPars: vector of parameters too short");
44
45 double dedxpar[3];
46 double nhitpar[6];
47 double cospar[11];
48 dedxpar[0] = 1; nhitpar[0] = 2; cospar[0] = 3;
49 for (int i = 0; i < 10; ++i) {
50 if (i < 2) dedxpar[i + 1] = params[i];
51 if (i < 5) nhitpar[i + 1] = params[i + 2];
52 cospar[i + 1] = params[i + 7];
53 }
54
55 // determine sigma from the parameterization
56 double corDedx = sigmaCurve(dedx, dedxpar, 0);
57
58 double nhit_min = 8, nhit_max = 37;
59 double corNHit = 0;
60 if (nhit < nhit_min) {
61 corNHit = sigmaCurve(nhit_min, nhitpar, 0) * std::sqrt(nhit_min / nhit);
62 } else if (nhit > nhit_max) {
63 corNHit = sigmaCurve(nhit_max, nhitpar, 0) * std::sqrt(nhit_max / nhit);
64 } else {
65 corNHit = sigmaCurve(nhit, nhitpar, 0);
66 }
67
68 double corCos = sigmaCurve(cosTheta, cospar, 0);
69
70 return (corDedx * corCos * corNHit * timereso);
71 }
Belle2::CDCDedxSigmaPars::m_sigmapars
std::vector< double > m_sigmapars
dE/dx resolution parameters
Definition: CDCDedxSigmaPars.h:87
Belle2::CDCDedxSigmaPars::sigmaCurve
double sigmaCurve(double x, const double *par, int version) const
parameterized resolution for predictions
Definition: CDCDedxSigmaPars.cc:19

◆ HELWIR()

void HELWIR ( const G4double  xwb4,
const G4double  ywb4,
const G4double  zwb4,
const G4double  xwf4,
const G4double  ywf4,
const G4double  zwf4,
const G4double  xp,
const G4double  yp,
const G4double  zp,
const G4double  px,
const G4double  py,
const G4double  pz,
const G4double  B_kG[3],
const G4double  charge,
const G4int  ntryMax,
G4double &  distance,
G4double  q2[3],
G4double  q1[3],
G4double  q3[3],
G4int &  ntry 
)
private

Calculate closest points between helix and wire.

Input xwb4 : x of wire at backward endplate in lab. ywb4 : y of wire at backward endplate " zwb4 : z of wire at backward endplate " xwf4 : x of wire at forward endplate " ywf4 : y of wire at forward endplate " zwf4 : z of wire at forward endplate " xp : x of helix in lab. yp : y of helix " zp : z of helix " px : px of helix in lab. py : py of helix " pz : pz of helix " Output q2(1) : x of wire at closest point in lab. q2(2) : y of wire at closest point " q2(3) : z of wire at closest point " q1(1) : x of helix at closest point " q1(2) : y of helix at closest point " q1(3) : z of helix at closest point " q3 : momentum of helix at closest point in lab. ntry :

Definition at line 1062 of file CDCSensitiveDetector.cc.

1076 {
1077 //~dead copy of gsim_cdc_hit.F in gsim-cdc for Belle (for tentative use)
1078 // ---------------------------------------------------------------------
1079 // Purpose : Calculate closest points between helix and wire.
1080 //
1081 // Input
1082 // xwb4 : x of wire at backward endplate in lab.
1083 // ywb4 : y of wire at backward endplate "
1084 // zwb4 : z of wire at backward endplate "
1085 // xwf4 : x of wire at forward endplate "
1086 // ywf4 : y of wire at forward endplate "
1087 // zwf4 : z of wire at forward endplate "
1088 //
1089 // Output
1090 // q2(1) : x of wire at closest point in lab.
1091 // q2(2) : y of wire at closest point "
1092 // q2(3) : z of wire at closest point "
1093 // q1(1) : x of helix at closest point "
1094 // q1(2) : y of helix at closest point "
1095 // q1(3) : z of helix at closest point "
1096 // ntry :
1097 // ---------------------------------------------------------------------
1098
1099 const G4int ndim = 3;
1100 const G4double delta = 1.e-5;
1101
1102
1103 G4double xwb, ywb, zwb, xwf, ywf, zwf;
1104 G4double xw, yw, zw, xh, yh, zh, pxh, pyh, pzh;
1105 G4double fi, fi_corr;
1106
1107 G4double dr, fi0, cpa, dz, tanl;
1108 G4double x0, y0, z0;
1109 // "chrg" removed by M. U. June, 2nd, 2013
1110 // G4double xc, yc, r, chrg;
1111 G4double xc, yc, r;
1112 G4double xwm, ywm;
1113 G4double sinfi0, cosfi0, sinfi0fi, cosfi0fi;
1114
1115 G4double vx, vy, vz, vv, cx, cy, cz, tt[3][3];
1116 G4double tmp[3];
1117
1118 G4double xx[3], dxx[3], ddxx[3], pp[3];
1119 G4double xxtdxx, dxxtdxx, xxtddxx;
1120
1121
1122 G4double fst = 0.0;
1123 G4double f, fderiv, deltafi, fact, eval;
1124 G4double dx1, dy1, dx2, dy2, crs, dot;
1125
1126 G4int iflg;
1127
1128 //set parameters
1129 xwb = xwb4; ywb = ywb4; zwb = zwb4;
1130 xwf = xwf4; ywf = ywf4; zwf = zwf4;
1131
1132 G4double xxx(xp), yyy(yp), zzz(zp);
1133 G4double pxx(px), pyy(py), pzz(pz);
1134
1135 //rotate z-axis to be parallel to B-field in case of non-uniform B
1136 Rotat(xwb, ywb, zwb, 1);
1137 Rotat(xwf, ywf, zwf, 1);
1138 Rotat(xxx, yyy, zzz, 1);
1139 Rotat(pxx, pyy, pzz, 1);
1140
1141 G4double a[8] = {0.};
1142 G4double pt = sqrt(pxx * pxx + pyy * pyy);
1143 a[1] = atan2(-pxx, pyy);
1144 a[2] = charge / pt;
1145 a[4] = pzz / pt;
1146 a[5] = xxx; a[6] = yyy; a[7] = zzz;
1147
1148 //calculate unit direction vector of the sense wire
1149 vx = xwf - xwb; vy = ywf - ywb; vz = zwf - zwb;
1150 vv = sqrt(vx * vx + vy * vy + vz * vz);
1151 vx /= vv; vy /= vv; vz /= vv;
1152
1153 //flag for distinguishing between stereo and axial wire
1154 iflg = 0;
1155 if (vx == 0. && vy == 0.) iflg = 1;
1156 // std::cout << "iflg= " << iflg << std::endl;
1157 //write(6,*) ' hlx2wir ', xwb, ywb, zwb, vx, vy, vz
1158
1159 //calculate coefficients of f
1160 cx = xwb - vx * (vx * xwb + vy * ywb + vz * zwb);
1161 cy = ywb - vy * (vx * xwb + vy * ywb + vz * zwb);
1162 cz = zwb - vz * (vx * xwb + vy * ywb + vz * zwb);
1163
1164 //calculate tensor for f
1165 tt[0][0] = vx * vx - 1.; tt[1][0] = vx * vy; tt[2][0] = vx * vz;
1166 tt[0][1] = vy * vx; tt[1][1] = vy * vy - 1.; tt[2][1] = vy * vz;
1167 tt[0][2] = vz * vx; tt[1][2] = vz * vy; tt[2][2] = vz * vz - 1.;
1168
1169 //set helix parameters
1170 dr = a[0]; fi0 = a[1]; cpa = a[2];
1171 dz = a[3]; tanl = a[4];
1172 x0 = a[5]; y0 = a[6]; z0 = a[7];
1173
1174 //
1175 // set initial value for phi
1176 //
1177
1178 xwm = xxx;
1179 ywm = yyy;
1180 //r(cm) = alpha/cpa = alpha * pt(GeV); bfield(kG)
1181 G4double bfield = sqrt(B_kG[0] * B_kG[0] +
1182 B_kG[1] * B_kG[1] +
1183 B_kG[2] * B_kG[2]);
1184 G4double alpha = 1.e4 / 2.99792458 / bfield;
1185 r = alpha / cpa;
1186 cosfi0 = cos(fi0);
1187 sinfi0 = sin(fi0);
1188 xc = x0 + (dr + r) * cosfi0;
1189 yc = y0 + (dr + r) * sinfi0;
1190 dx1 = x0 - xc;
1191 dy1 = y0 - yc;
1192 dx2 = xwm - xc;
1193 dy2 = ywm - yc;
1194 crs = dx1 * dy2 - dy1 * dx2;
1195 dot = dx1 * dx2 + dy1 * dy2;
1196 fi = atan2(crs, dot);
1197
1198 //begin iterative procedure for newton 's method '
1199 fact = 1.;
1200 ntry = 0;
1201line1:
1202 ntry += 1;
1203 cosfi0fi = cos(fi0 + fi);
1204 sinfi0fi = sin(fi0 + fi);
1205
1206 //calculate spatial point Q(x,y,z) along the helix
1207 xx[0] = x0 + dr * cosfi0 + r * (cosfi0 - cosfi0fi);
1208 xx[1] = y0 + dr * sinfi0 + r * (sinfi0 - sinfi0fi);
1209 xx[2] = z0 + dz - r * tanl * fi;
1210 pp[0] = -pt * sinfi0fi;
1211 pp[1] = pt * cosfi0fi;
1212 pp[2] = pt * tanl;
1213
1214 if (iflg == 1) {
1215 q2[0] = xwb; q2[1] = ywb; q2[2] = xx[2];
1216 q1[0] = xx[0]; q1[1] = xx[1]; q1[2] = xx[2];
1217 q3[0] = pp[0]; q3[1] = pp[1]; q3[2] = pp[2];
1218 //inverse rotation to lab. frame in case of non-uniform B
1219 Rotat(q1, -1);
1220 Rotat(q2, -1);
1221 Rotat(q3, -1);
1222 distance = sqrt((q2[0] - q1[0]) * (q2[0] - q1[0]) +
1223 (q2[1] - q1[1]) * (q2[1] - q1[1]) +
1224 (q2[2] - q1[2]) * (q2[2] - q1[2]));
1225 return;
1226 }
1227
1228 //calculate direction vector (dx/dphi,dy/dphi,dz/dphi)
1229 //on a point along the helix.
1230 dxx[0] = r * sinfi0fi; dxx[1] = - r * cosfi0fi; dxx[2] = - r * tanl;
1231
1232 // In order to derive the closest point between straight line and helix,
1233 // we can put following two conditions:
1234 // (i) A point H(xh,yh,zh) on the helix given should be on
1235 // the plane which is perpendicular to the straight line.
1236 // (ii) A line HW from W(xw,yw,zw) which is a point on the straight
1237 // line to H(xh,yh,zh) should normal to the direction vector
1238 // on the point H.
1239 //
1240 // Thus, we can make a equation from above conditions.
1241 // f(phi) = cx*(dx/dphi) + cy*(dy/dphi) + cz*(dz/dphi)
1242 // + (x,y,z)*tt(i,j)*(dx/dphi,dy/dphi,dz/dphi)
1243 // = 0,
1244 // where
1245 // cx = xwb - vx*( vx*xwb + vy*ywb + vz*zwb )
1246 // cy = ywb - vy*( vx*xwb + vy*ywb + vz*zwb )
1247 // cz = zwb - vz*( vx*xwb + vy*ywb + vz*zwb )
1248 //
1249 // tt(1,1) = vx*vx - 1 tt(1,2) = vx*vy tt(1,3) = vx*vz
1250 // tt(2,1) = vy*vx tt(2,2) = vy*vy - 1 tt(2,3) = vy*vz
1251 // tt(3,1) = vz*vx tt(3,2) = vz*vy tt(3,3) = vz*vz - 1
1252 //
1253 // and the equation of straight line(stereo wire) is written by
1254 // (x,y,z) = (xwb,ywb,zwb) + beta*(vx,vy,vz), beta is free parameter.
1255
1256 //Now calculate f
1257 Mvopr(ndim, xx, tt, dxx, tmp, 1);
1258 xxtdxx = tmp[0];
1259 f = cx * dxx[0] + cy * dxx[1] + cz * dxx[2] + xxtdxx;
1260 if (std::abs(f) < delta) goto line100;
1261
1262 //evaluate fitting result and prepare some factor to multiply to 1/derivative
1263 if (ntry > 1) {
1264 eval = (1.0 - 0.25 * fact) * std::abs(fst) - std::abs(f);
1265 if (eval <= 0.) fact *= 0.5;
1266 }
1267
1268 //calculate derivative of f
1269 ddxx[0] = r * cosfi0fi; ddxx[1] = r * sinfi0fi; ddxx[2] = 0.;
1270
1271 //Now we have derivative of f
1272 Mvopr(ndim, dxx, tt, dxx, tmp, 1);
1273 dxxtdxx = tmp[0];
1274 Mvopr(ndim, xx, tt, ddxx, tmp, 1);
1275 xxtddxx = tmp[0];
1276 fderiv = cx * ddxx[0] + cy * ddxx[1] + cz * ddxx[2] + dxxtdxx + xxtddxx;
1277 // Commented by M. U. June, 2nd, 2013
1278 // fist = fi;
1279 deltafi = f / fderiv;
1280 fi -= fact * deltafi;
1281 fst = f;
1282
1283 if (ntry > ntryMax) {
1284 //B2DEBUG(" Exceed max. trials HelWir ");
1285 goto line100;
1286 }
1287 //write(6,*) ntry, fist, deltafi
1288 goto line1;
1289
1290 //check if zh is btw zwb and zwf; if not, set zh=zwb or zh=zwf.
1291 //dead regions due to feed-throughs should be considered later.
1292line100:
1293 zh = z0 + dz - r * tanl * fi;
1294 fi_corr = 0.;
1295 if (zh < zwb) fi_corr = (zwb - zh) / (-r * tanl);
1296 if (zh > zwf) fi_corr = (zwf - zh) / (-r * tanl);
1297 fi += fi_corr;
1298
1299 cosfi0fi = cos(fi0 + fi);
1300 sinfi0fi = sin(fi0 + fi);
1301
1302 xh = x0 + dr * cosfi0 + r * (cosfi0 - cosfi0fi);
1303 yh = y0 + dr * sinfi0 + r * (sinfi0 - sinfi0fi);
1304 zh = z0 + dz - r * tanl * fi;
1305 pxh = -pt * sinfi0fi;
1306 pyh = pt * cosfi0fi;
1307 pzh = pt * tanl;
1308
1309 //write(6,*) 'fi_corr, zh, zwb, zwf=', fi_corr, zh, zwb, zwf
1310 //write(6,*) 'zh = ', z0, dz, r, tanl, fi
1311
1312 zw = vx * vz * xh + vy * vz * yh + vz * vz * zh + zwb - vz * (vx * xwb + vy * ywb + vz * zwb);
1313 xw = xwb + vx * (zw - zwb) / vz;
1314 yw = ywb + vy * (zw - zwb) / vz;
1315
1316 q2[0] = xw; q2[1] = yw; q2[2] = zw;
1317 q1[0] = xh; q1[1] = yh; q1[2] = zh;
1318 q3[0] = pxh; q3[1] = pyh; q3[2] = pzh;
1319
1320 //inverse rotation to lab. frame in case of non-uniform B
1321 Rotat(q1, -1);
1322 Rotat(q2, -1);
1323 Rotat(q3, -1);
1324 distance = sqrt((q2[0] - q1[0]) * (q2[0] - q1[0]) +
1325 (q2[1] - q1[1]) * (q2[1] - q1[1]) +
1326 (q2[2] - q1[2]) * (q2[2] - q1[2]));
1327 return;
1328
1329 }
Belle2::CDC::CDCSensitiveDetector::Rotat
void Rotat(G4double &x, G4double &y, G4double &z, const int mode)
Translation method.
Definition: CDCSensitiveDetector.cc:1007
Belle2::CDC::CDCSensitiveDetector::Mvopr
void Mvopr(const G4int ndim, const G4double b[3], const G4double m[3][3], const G4double a[3], G4double c[3], const G4int mode)
Calculate the result of a matrix times vector.
Definition: CDCSensitiveDetector.cc:1332
Belle2::dot
T dot(GeneralVector< T > a, GeneralVector< T > b)
dot product of two general vectors
Definition: beamHelpers.h:163
Belle2::eval
double eval(const std::vector< double > &spl, const std::vector< double > &vals, double x)
Evaluate spline (zero order or first order) in point x.
Definition: tools.h:115
Belle2::EvtPDLUtil::charge
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44

◆ I2D()

double I2D ( double  cosTheta,
double  I 
) const
private

hadron saturation parameterization part 1

Definition at line 47 of file CDCDedxHadronCor.cc.

48 {
49 const auto& params = m_hadronpars;
50 if (params.size() < 5) {
51 B2WARNING("Vector of dE/dx hadron constants too short!");
52 return I;
53 }
54
55 double projection = std::pow(std::abs(cosTheta), params[3]) + params[2];
56 if (projection == 0 or params[4] == 0) {
57 B2WARNING("Something wrong with dE/dx hadron constants!");
58 return I;
59 }
60
61 double a = params[0] / projection;
62 double b = 1 - params[1] / projection * (I / params[4]);
63 double c = -1.0 * I / params[4];
64
65 if (b == 0 and a == 0) {
66 B2WARNING("both a and b coefficiants for hadron correction are 0");
67 return I;
68 }
69
70 double discr = b * b - 4.0 * a * c;
71 if (discr < 0) {
72 B2WARNING("negative discriminant; return uncorrectecd value");
73 return I;
74 }
75
76 double D = (a != 0) ? (-b + std::sqrt(discr)) / (2.0 * a) : -c / b;
77 if (D < 0) {
78 D = (a != 0) ? (-b - std::sqrt(discr)) / (2.0 * a) : -c / b;
79 if (D < 0) {
80 B2WARNING("D is less than 0; return uncorrectecd value");
81 return I;
82 }
83 }
84
85 return D;
86 }

◆ Initialize()

void Initialize ( G4HCofThisEvent *  )
override

Register CDC hits collection into G4HCofThisEvent.

Definition at line 89 of file CDCSensitiveDetector.cc.

90 {
91 /*
92 m_cdcgp = &CDCGeometryPar::Instance();
93
94 m_thresholdEnergyDeposit = m_cdcgp->getThresholdEnergyDeposit();
95 m_thresholdEnergyDeposit *= CLHEP::GeV; //GeV to MeV (=unit in G4)
96 B2INFO("CDCSensitiveDetector: Threshold energy (MeV): " << m_thresholdEnergyDeposit);
97 m_modifiedLeftRightFlag = m_cdcgp->isModifiedLeftRightFlagOn();
98 B2INFO("CDCSensitiveDetector: Set left/right flag modified for tracking (=1)/ not set (=0): " << m_modifiedLeftRightFlag);
99
100 m_minTrackLength = m_cdcgp->getMinTrackLength();
101 m_minTrackLength *= CLHEP::cm; //cm to mm (=unit in G4)
102 B2INFO("CDCSensitiveDetector: MinTrackLength (mm): " << m_minTrackLength);
103 */
104
105 // Initialize
106 m_nonUniformField = 0;
107 // std::cout << "Initialize called" << std::endl;
108 // m_nPosHits = 0;
109 // m_nNegHits = 0;
110 }

◆ meanCurve()

double meanCurve ( double  x,
const double *  par,
int  version 
) const
private

parameterized beta-gamma curve for predicted means

Definition at line 19 of file CDCDedxMeanPars.cc.

20 {
21 // calculate the predicted mean value as a function of beta-gamma (bg)
22 // this is done with a different function depending on the value of bg
23 double f = 0;
24
25 if (version == 0) {
26 if (par[0] == 1)
27 f = par[1] * std::pow(std::sqrt(x * x + 1), par[3]) / std::pow(x, par[3]) *
28 (par[2] - par[5] * std::log(1 / x)) - par[4] + std::exp(par[6] + par[7] * x);
29 else if (par[0] == 2)
30 f = par[1] * x * x * x + par[2] * x * x + par[3] * x + par[4];
31 else if (par[0] == 3)
32 f = -1.0 * par[1] * std::log(par[4] + std::pow(1 / x, par[2])) + par[3];
33 }
34
35 return f;
36 }

◆ Mvopr()

void Mvopr ( const G4int  ndim,
const G4double  b[3],
const G4double  m[3][3],
const G4double  a[3],
G4double  c[3],
const G4int  mode 
)
private

Calculate the result of a matrix times vector.

Input ndim : dimension b(1-ndim) : vector m(1-ndim,1-ndim) : matrix a(1-ndim) : vector c(1-ndim) : vector mode : c = m * a for mode=0 c = b * m * a for mode=1 Output c(1-ndim) : for mode 1, solution is put on c[0]

Definition at line 1332 of file CDCSensitiveDetector.cc.

1334 {
1335 //~dead copy of UtilCDC_mvopr in com-cdc for Belle (for tentative use)
1336 //-----------------------------------------------------------------------
1337 // Input
1338 // ndim : dimension
1339 // b(1-ndim) : vector
1340 // m(1-ndim,1-ndim) : matrix
1341 // a(1-ndim) : vector
1342 // c(1-ndim) : vector
1343 // mode : c = m * a for mode=0
1344 // c = b * m * a for mode=1
1345 // Output
1346 // c(1-ndim) : for mode 1, solution is put on c[0]
1347 //-----------------------------------------------------------------------
1348
1349 if (ndim != 3) {
1350 //B2ERROR("invalid ndim " << ndim << " specified");
1351 return;
1352 }
1353
1354 for (int i = 0; i < ndim; ++i) c[i] = 0.;
1355 G4double tmp[3];
1356 for (int i = 0; i < ndim; ++i) tmp[i] = 0.;
1357
1358 if (mode == 0) {
1359 for (int i = 0; i < ndim; ++i) {
1360 for (int j = 0; j < ndim; ++j) {
1361 c[i] += m[j][i] * a[j];
1362 }
1363 }
1364 return;
1365 } else if (mode == 1) {
1366 for (int i = 0; i < ndim; ++i) {
1367 for (int j = 0; j < ndim; ++j) {
1368 tmp[i] += m[j][i] * a[j];
1369 }
1370 c[0] += b[i] * tmp[i];
1371 }
1372 } else {
1373 //B2ERROR("Error, you specified invalid mode= " << mode);
1374 }
1375
1376 return;
1377
1378 }

◆ reAssignLeftRightInfo()

void reAssignLeftRightInfo ( )
private

Re-assign left/right info.

Definition at line 1523 of file CDCSensitiveDetector.cc.

1524 {
1525 CDCSimHit* sHit = nullptr;
1526 WireID sWireId; // = WireID();
1527 B2Vector3D sPos; // = B2Vector3D();
1528
1529 CDCSimHit* pHit = nullptr;
1530 WireID pWireId; // = WireID();
1531 // double minDistance2 = DBL_MAX;
1532 // double distance2 = DBL_MAX;
1533 // unsigned short bestNeighb = 0;
1534 // unsigned short neighb = 0;
1535
1536 // std::multimap<unsigned short, CDCSimHit*>::iterator pItBegin = m_hitWithPosWeight.begin();
1537 // std::multimap<unsigned short, CDCSimHit*>::iterator pItEnd = m_hitWithPosWeight.end();
1538
1539 // unsigned short sClayer = 0;
1540 // unsigned short sSuperLayer = 0;
1541 // unsigned short sLayer = 0;
1542 // unsigned short sWire = 0;
1543 // CDCSimHit* fHit = nullptr;
1544
1545 //Find a primary track close to the input 2'ndary hit in question
1546 for (std::vector<CDCSimHit*>::iterator nIt = m_hitWithNegWeight.begin(), nItEnd = m_hitWithNegWeight.end(); nIt != nItEnd; ++nIt) {
1547
1548 sHit = *nIt;
1549 sPos = sHit->getPosTrack();
1550 sWireId = sHit->getWireID();
1551 // sClayer = sWireId.getICLayer();
1552 // sSuperLayer = sWireId.getISuperLayer();
1553 // sLayer = sWireId.getILayer();
1554 // sWire = sWireId.getIWire();
1555 // fHit = sHit;
1556 unsigned short sClayer = sWireId.getICLayer();
1557 unsigned short sSuperLayer = sWireId.getISuperLayer();
1558 unsigned short sLayer = sWireId.getILayer();
1559 unsigned short sWire = sWireId.getIWire();
1560 CDCSimHit* fHit = sHit;
1561
1562 std::multimap<unsigned short, CDCSimHit*>::iterator pItBegin = m_hitWithPosWeight.find(sSuperLayer);
1563 std::multimap<unsigned short, CDCSimHit*>::iterator pItEnd = m_hitWithPosWeight.find(sSuperLayer + 1);
1564 /*
1565 if (sSuperLayer <= 8) {
1566 pItBegin = m_posWeightMapItBegin.at(sSuperLayer);
1567 pItEnd = m_posWeightMapItEnd.at(sSuperLayer);
1568 } else {
1569 B2FATAL("CDCSensitiveDetector::EndOfEvent: invalid super-layer id ! " << sSuperLayer);
1570 }
1571 */
1572
1573 double minDistance2 = DBL_MAX;
1574 // bestNeighb = 0;
1575
1576 /* for (std::multimap<unsigned short, CDCSimHit*>::iterator pIt = m_hitWithPosWeight.begin(); pIt != m_hitWithPosWeight.end(); ++pIt) {
1577 std::cout <<"superLyr#= " << pIt->first << std::endl;
1578 }
1579 */
1580
1581 for (std::multimap<unsigned short, CDCSimHit*>::iterator pIt = pItBegin; pIt != pItEnd; ++pIt) {
1582
1583 //scan hits in the same/neighboring cells
1584 pHit = pIt->second;
1585 pWireId = pHit->getWireID();
1586 // neigh = areNeighbors(sWireId, pWireId);
1587 unsigned short neighb = areNeighbors(sClayer, sSuperLayer, sLayer, sWire, pWireId);
1588 if (neighb != 0 || pWireId == sWireId) {
1589 double distance2 = (pHit->getPosTrack() - sPos).Mag2();
1590 if (distance2 < minDistance2) {
1591 fHit = pHit;
1592 minDistance2 = distance2;
1593 // bestNeighb = neighb;
1594 }
1595 }
1596 }
1597
1598 //reassign LR using the momentum-direction of the primary particle found
1599 unsigned short lR = m_cdcgp->getNewLeftRightRaw(sHit->getPosWire(),
1600 sHit->getPosTrack(),
1601 fHit->getMomentum());
1602 // unsigned short bflr = sHit->getLeftRightPassage();
1603 sHit->setLeftRightPassage(lR);
1604 // std::cout <<"neighb, bfaf lrs, minDistance= " << bestNeighb <<" "<<" "<< bflr <<" "<< sHit->getLeftRightPassage() <<" "<< std::scientific << sqrt(minDistance2) << std::endl;
1605 }
1606 }
Belle2::CDC::CDCGeometryPar::getNewLeftRightRaw
unsigned short getNewLeftRightRaw(const B2Vector3D &posOnWire, const B2Vector3D &posOnTrack, const B2Vector3D &momentum) const
Returns new left/right_raw.
Definition: CDCGeometryPar.cc:2786
Belle2::CDC::CDCSensitiveDetector::m_hitWithPosWeight
std::multimap< unsigned short, CDCSimHit * > m_hitWithPosWeight
Map containing hits with positive weight.
Definition: CDCSensitiveDetector.h:280
Belle2::CDC::CDCSensitiveDetector::m_hitWithNegWeight
std::vector< CDCSimHit * > m_hitWithNegWeight
Vector containing hits with negative weight.
Definition: CDCSensitiveDetector.h:282
Belle2::CDC::CDCSensitiveDetector::areNeighbors
unsigned short areNeighbors(const WireID &wireId, const WireID &otherWireId) const
Check if neighboring cell in the same super-layer; essentially a copy from cdcLocalTracking/mclookup.
Definition: CDCSensitiveDetector.cc:1609

◆ Rotat() [1/2]

void Rotat ( G4double &  x,
G4double &  y,
G4double &  z,
const int  mode 
)
private

Translation method.

Translates (x,y,z) in lab. to (x,y,z) in B-field frame (mode=1), or reverse translation (mode=-1).

Definition at line 1007 of file CDCSensitiveDetector.cc.

1009 {
1010 //Translates (x,y,z) in lab. to (x,y,z) in B-field frame (mode=1), or reverse
1011 // translation (mode=-1).
1012 //~dead copy (for tentative use) of gsim_cdc_rotat/irotat.F in gsim-cdc
1013 //for Belle
1014
1015 if (m_nonUniformField == 0) return;
1016
1017 G4double x0(x), y0(y), z0(z);
1018
1019 if (mode == 1) {
1020 x = m_brot[0][0] * x0 + m_brot[1][0] * y0 + m_brot[2][0] * z0;
1021 y = m_brot[0][1] * x0 + m_brot[1][1] * y0 + m_brot[2][1] * z0;
1022 z = m_brot[0][2] * x0 + m_brot[1][2] * y0 + m_brot[2][2] * z0;
1023 } else if (mode == -1) {
1024 x = m_brot[0][0] * x0 + m_brot[0][1] * y0 + m_brot[0][2] * z0;
1025 y = m_brot[1][0] * x0 + m_brot[1][1] * y0 + m_brot[1][2] * z0;
1026 z = m_brot[2][0] * x0 + m_brot[2][1] * y0 + m_brot[2][2] * z0;
1027 } else {
1028 //B2ERROR("SensitiveDetector " <<"invalid mode " << mode << "specified");
1029 }
1030 return;
1031
1032 }

◆ Rotat() [2/2]

void Rotat ( G4double  x[3],
const int  mode 
)
private

Overloaded translation method.

Definition at line 1035 of file CDCSensitiveDetector.cc.

1036 {
1037 //Translates (x,y,z) in lab. to (x,y,z) in B-field frame (mode=1), or reverse
1038 // translation (mode=-1).
1039 //~dead copy (for tentative use) of gsim_cdc_rotat/irotat.F in gsim-cdc
1040 //for Belle
1041
1042 if (m_nonUniformField == 0) return;
1043
1044 G4double x0(x[0]), y0(x[1]), z0(x[2]);
1045
1046 if (mode == 1) {
1047 x[0] = m_brot[0][0] * x0 + m_brot[1][0] * y0 + m_brot[2][0] * z0;
1048 x[1] = m_brot[0][1] * x0 + m_brot[1][1] * y0 + m_brot[2][1] * z0;
1049 x[2] = m_brot[0][2] * x0 + m_brot[1][2] * y0 + m_brot[2][2] * z0;
1050 } else if (mode == -1) {
1051 x[0] = m_brot[0][0] * x0 + m_brot[0][1] * y0 + m_brot[0][2] * z0;
1052 x[1] = m_brot[1][0] * x0 + m_brot[1][1] * y0 + m_brot[1][2] * z0;
1053 x[2] = m_brot[2][0] * x0 + m_brot[2][1] * y0 + m_brot[2][2] * z0;
1054 } else {
1055 //B2ERROR("SensitiveDetector " <<"invalid mode " << mode << "specified");
1056 }
1057 return;
1058
1059 }

◆ saveSimHit()

void saveSimHit ( const G4int  layerId,
const G4int  wireId,
const G4int  trackID,
const G4int  pid,
const G4double  distance,
const G4double  tof,
const G4double  edep,
const G4double  stepLength,
const G4ThreeVector &  mom,
const G4ThreeVector &  posW,
const G4ThreeVector &  posIn,
const G4ThreeVector &  posOut,
const G4ThreeVector &  posTrack,
const G4int  lr,
const G4int  NewLrRaw,
const G4int  NewLr,
const G4double  speed,
const G4double  hitWeight 
)

Save CDCSimHit into datastore.

Definition at line 586 of file CDCSensitiveDetector.cc.

604 {
605
606 // Discard the hit below Edep_th
607 // if (edep <= m_thresholdEnergyDeposit) return 0;
608 if (edep <= m_thresholdEnergyDeposit) return;
609
610 //compute tof at the closest point; linear approx.
611 const G4double sign = (posTrack - posIn).dot(mom) < 0. ? -1. : 1.;
612 const G4double CorrectTof = tof + sign * (posTrack - posIn).mag() / speed;
613 // if (sign < 0.) std::cout <<"deltatof= "<< sign * (posTrack - posIn).mag() / speed << std::endl;
614#if defined(CDC_DEBUG)
615 std::cout << "posIn= " << posIn.x() << " " << posIn.y() << " " << posIn.z() << std::endl;
616 std::cout << "posOut= " << posOut.x() << " " << posOut.y() << " " << posOut.z() << std::endl;
617 std::cout << "posTrack= " << posTrack.x() << " " << posTrack.y() << " " << posTrack.z() << std::endl;
618 std::cout << "posW= " << posW.x() << " " << posW.y() << " " << posW.z() << std::endl;
619 std::cout << "tof = " << tof << std::endl;
620 std::cout << "deltaTof = " << (posTrack - posIn).mag() / speed << std::endl;
621 std::cout << "CorrectTof= " << CorrectTof << std::endl;
622 if (CorrectTof > 95) {
623 std::cout << "toolargecorrecttof" << std::endl;
624 }
625#endif
626
627 RelationArray cdcSimHitRel(m_MCParticles, m_CDCSimHits);
628
629 m_hitNumber = m_CDCSimHits.getEntries();
630
631 CDCSimHit* simHit = m_CDCSimHits.appendNew();
632
633 simHit->setWireID(layerId, wireId);
634 simHit->setTrackId(trackID);
635 simHit->setPDGCode(pid);
636 simHit->setDriftLength(distance / CLHEP::cm);
637 simHit->setFlightTime(CorrectTof / CLHEP::ns);
638 simHit->setGlobalTime(CorrectTof / CLHEP::ns);
639 simHit->setEnergyDep(edep / CLHEP::GeV);
640 simHit->setStepLength(stepLength / CLHEP::cm);
641 B2Vector3D momentum(mom.getX() / CLHEP::GeV, mom.getY() / CLHEP::GeV, mom.getZ() / CLHEP::GeV);
642 simHit->setMomentum(momentum);
643 B2Vector3D posWire(posW.getX() / CLHEP::cm, posW.getY() / CLHEP::cm, posW.getZ() / CLHEP::cm);
644 simHit->setPosWire(posWire);
645 B2Vector3D positionIn(posIn.getX() / CLHEP::cm, posIn.getY() / CLHEP::cm, posIn.getZ() / CLHEP::cm);
646 simHit->setPosIn(positionIn);
647 B2Vector3D positionOut(posOut.getX() / CLHEP::cm, posOut.getY() / CLHEP::cm, posOut.getZ() / CLHEP::cm);
648 simHit->setPosOut(positionOut);
649 B2Vector3D positionTrack(posTrack.getX() / CLHEP::cm, posTrack.getY() / CLHEP::cm, posTrack.getZ() / CLHEP::cm);
650 simHit->setPosTrack(positionTrack);
651 simHit->setPosFlag(lr);
652 simHit->setLeftRightPassageRaw(newLrRaw);
653 simHit->setLeftRightPassage(newLr);
654#if defined(CDC_DEBUG)
655 std::cout << "sensitived,oldlr,newlrRaw,newlr= " << lr << " " << newLrRaw << " " << newLr << std::endl;
656#endif
657
658 B2DEBUG(150, "HitNumber: " << m_hitNumber);
659 if (m_modifiedLeftRightFlag) {
660 //N.B. Negative hitWeight is allowed intentionally here; all weights are to be reset to positive in EndOfEvent
661 cdcSimHitRel.add(trackID, m_hitNumber, hitWeight);
662 } else {
663 cdcSimHitRel.add(trackID, m_hitNumber);
664 }
665
666 // if (hitWeight > 0) m_nPosHits++;
667 // if (hitWeight < 0) m_nNegHits++;
668 // std::cout <<"trackID,HitNumber,weight,driftL,edep= "<< trackID <<" "<< m_hitNumber <<" "<< hitWeight <<" "<< distance <<" "<< edep << std::endl;
669 // return (m_hitNumber);
670 }

◆ setModifiedLeftRightFlag()

void setModifiedLeftRightFlag ( )
private

set left/right flag modified for tracking

Definition at line 1439 of file CDCSensitiveDetector.cc.

1440 {
1441 if (!m_modifiedLeftRightFlag) return;
1442
1443 // std::cout <<"#posHits,#negHits= " << m_nPosHits <<" "<< m_nNegHits << std::endl;
1444
1445 // Get SimHit array and relation betw. MC and SimHit
1446 // N.B. MCParticle is incomplete at this stage; the relation betw it and
1447 // simHit is Okay.
1448 // MCParticle will be completed after all sub-detectors' EndOfEvent calls.
1449 RelationArray mcPartToSimHits(m_MCParticles, m_CDCSimHits);
1450 int nRelationsMinusOne = mcPartToSimHits.getEntries() - 1;
1451
1452 if (nRelationsMinusOne == -1) return;
1453
1454 // std::cout <<"#simHits= " << m_CDCSimHits.getEntries() << std::endl;
1455 // std::cout <<"#MCParticles= " << m_MCParticles.getEntries() << std::endl;
1456 // std::cout <<"#mcPartToSimHits= " << mcPartToSimHits.getEntries() << std::endl;
1457
1458 //reset some of negative weights to positive; this is needed for the hits
1459 //created by secondary particles whose track-lengths get larger than the
1460 //threshold (set by the user) during G4 swimming (i.e. the weights are
1461 //first set to negative as far as the track-lengths are shorther than the
1462 //threshold; set to positive when the track-lengths exceed the threshold).
1463
1464 size_t iRelation = 0;
1465 int trackIdOld = INT_MAX;
1466 // std::cout << "INT_MAX= " << INT_MAX << std::endl;
1467 m_hitWithPosWeight.clear();
1468 m_hitWithNegWeight.clear();
1469
1470 for (int it = nRelationsMinusOne; it >= 0; --it) {
1471 RelationElement& mcPartToSimHit = const_cast<RelationElement&>(mcPartToSimHits[it]);
1472 size_t nRelatedHits = mcPartToSimHit.getSize();
1473 if (nRelatedHits > 1) B2FATAL("CDCSensitiveDetector::EndOfEvent: MCParticle<-> CDCSimHit relation is not one-to-one !");
1474
1475 unsigned short trackId = mcPartToSimHit.getFromIndex();
1476 RelationElement::weight_type weight = mcPartToSimHit.getWeight(iRelation);
1477 if (weight > 0.) {
1478 trackIdOld = trackId;
1479 } else if (weight <= 0. && trackId == trackIdOld) {
1480 // RelationElement::index_type iSimHit = mcPartToSimHit.getToIndex(iRelation);
1481 weight *= -1.;
1482 mcPartToSimHit.setToIndex(mcPartToSimHit.getToIndex(iRelation), weight);
1483 trackIdOld = trackId;
1484 // std::cout <<"trackId,,iSimHit,wgtafterreset= "<< trackId <<" "<< iSimHit <<" "<< mcPartToSimHit.getWeight(iRelation) << std::endl;
1485 }
1486
1487 CDCSimHit* sHit = m_CDCSimHits[mcPartToSimHit.getToIndex(iRelation)];
1488
1489 if (weight > 0.) {
1490 m_hitWithPosWeight.insert(std::pair<unsigned short, CDCSimHit*>(sHit->getWireID().getISuperLayer(), sHit));
1491 } else {
1492 m_hitWithNegWeight.push_back(sHit);
1493 }
1494 }
1495
1496 /*
1497 // std::cout <<"m_hitWithPosWeight.size= " << m_hitWithPosWeight.size() << std::endl;
1498 for(int i=0; i<9; ++i) {
1499 // if (m_hitWithPosWeight.find(i) != m_hitWithPosWeight.end()) {
1500 // std::cout << i << " found" << std::endl;
1501 // }
1502 m_posWeightMapItBegin.push_back(m_hitWithPosWeight.find(i));
1503 m_posWeightMapItEnd.push_back(m_hitWithPosWeight.find(i+1));
1504 }
1505 */
1506
1507 //reassign L/R flag
1508 reAssignLeftRightInfo();
1509
1510 //reset all weights positive; this is required for completing MCParticle object at the EndOfEvent action of FullSim
1511 // is this part really needed ??? check again !
1512 for (int it = 0; it <= nRelationsMinusOne; ++it) {
1513 RelationElement& mcPartToSimHit = const_cast<RelationElement&>(mcPartToSimHits[it]);
1514 RelationElement::weight_type weight = mcPartToSimHit.getWeight(iRelation);
1515 if (weight < 0.) {
1516 mcPartToSimHit.setToIndex(mcPartToSimHit.getToIndex(iRelation), -1.*weight);
1517 }
1518 }
1519
1520 }
Belle2::RelationElement::weight_type
float weight_type
type used for weights.
Definition: RelationElement.h:32
Belle2::CDC::CDCSensitiveDetector::reAssignLeftRightInfo
void reAssignLeftRightInfo()
Re-assign left/right info.
Definition: CDCSensitiveDetector.cc:1523

◆ sigmaCurve()

double sigmaCurve ( double  x,
const double *  par,
int  version 
) const
private

parameterized resolution for predictions

Definition at line 19 of file CDCDedxSigmaPars.cc.

20 {
21 // calculate the predicted mean value as a function of beta-gamma (bg)
22 // this is done with a different function depending dE/dx, nhit, and sin(theta)
23 double f = 0;
24
25 if (version == 0) {
26 if (par[0] == 1) { // return dedx parameterization
27 f = par[1] + par[2] * x;
28 } else if (par[0] == 2) { // return nhit or sin(theta) parameterization
29 f = par[1] * x * x * x * x + par[2] * x * x * x + par[3] * x * x + par[4] * x + par[5];
30 } else if (par[0] == 3) { // return cos(theta) parameterization
31 f = par[1] * std::exp(-0.5 * std::pow(((x - par[2]) / par[3]), 2)) +
32 par[4] * std::pow(x, 6) + par[5] * std::pow(x, 5) + par[6] * std::pow(x, 4) +
33 par[7] * x * x * x + par[8] * x * x + par[9] * x + par[10];
34 }
35 }
36
37 return f;
38 }

◆ step()

bool step ( G4Step *  aStep,
G4TouchableHistory *  history 
)
overridevirtual

Process each step and calculate variables defined in CDCB4VHit.

Implements SensitiveDetectorBase.

Definition at line 115 of file CDCSensitiveDetector.cc.

116 {
117 // static bool firstCall = true;
118 // if (firstCall) {
119 // firstCall = false;
120 m_cdcgp = &CDCGeometryPar::Instance();
121 // CDCSimControlPar & m_cntlp = CDCSimControlPar::getInstance();
122
123 // // m_thresholdEnergyDeposit = m_cdcgp->getThresholdEnergyDeposit();
124 // m_thresholdEnergyDeposit = m_cntlp.getThresholdEnergyDeposit();
125 // m_thresholdEnergyDeposit *= CLHEP::GeV; //GeV to MeV (=unit in G4)
126 // B2INFO("CDCSensitiveDetector: Threshold energy (MeV): " << m_thresholdEnergyDeposit);
127
128 // // m_modifiedLeftRightFlag = m_cdcgp->isModifiedLeftRightFlagOn();
129 // m_modifiedLeftRightFlag = m_cntlp.getModLeftRightFlag();
130 // B2INFO("CDCSensitiveDetector: Set left/right flag modified for tracking (=1)/ not set (=0): " << m_modifiedLeftRightFlag);
131
132 // // m_minTrackLength = m_cdcgp->getMinTrackLength();
133 // m_minTrackLength = m_cntlp.getMinTrackLength();
134 // m_minTrackLength *= CLHEP::cm; //cm to mm (=unit in G4)
135 // B2INFO("CDCSensitiveDetector: MinTrackLength (mm): " << m_minTrackLength);
136
137 // // m_wireSag = m_cdcgp->isWireSagOn();
138 // m_wireSag = m_cntlp.getWireSag();
139
140 m_nonUniformField = 0;
141 // }
142
143#if defined(CDC_DEBUG)
144 std::cout << " " << std::endl;
145 std::cout << "********* step in ********" << std::endl;
146#endif
147 // Get deposited energy
148 const G4double edep = aStep->GetTotalEnergyDeposit();
149
150
151 // Discard the hit below Edep_th
152 if (edep <= m_thresholdEnergyDeposit) return false;
153
154 // Get step length
155 const G4double stepLength = aStep->GetStepLength();
156 if (stepLength == 0.) return false;
157
158 // Get step information
159 const G4Track& t = * aStep->GetTrack();
160
161 G4double hitWeight = Simulation::TrackInfo::getInfo(t).getIgnore() ? -1 : 1;
162 // save in MCParticle if track-length is enough long
163 if (t.GetTrackLength() > m_minTrackLength) {
164 Simulation::TrackInfo::getInfo(t).setIgnore(false);
165 hitWeight = 1.;
166 // std::cout <<"setignore=false for track= "<< t.GetTrackID() << std::endl;
167 // } else {
168 // std::cout <<"setignore=true for track= "<< t.GetTrackID() << std::endl;
169 }
170
171 // const G4double tof = t.GetGlobalTime(); //tof at post step point
172 // if (isnan(tof)) {
173 // B2ERROR("SensitiveDetector: global time is nan");
174 // return false;
175 // }
176
177 const G4int pid = t.GetDefinition()->GetPDGEncoding();
178 const G4double charge = t.GetDefinition()->GetPDGCharge();
179 const G4int trackID = t.GetTrackID();
180 // std::cout << "pid,stepl,trackID,trackl,weight= " << pid <<" "<< stepLength <<" "<< trackID <<" "<< t.GetTrackLength() <<" "<< hitWeight << std::endl;
181
182 const G4VPhysicalVolume& v = * t.GetVolume();
183 const G4StepPoint& in = * aStep->GetPreStepPoint();
184 const G4StepPoint& out = * aStep->GetPostStepPoint();
185 const G4ThreeVector& posIn = in.GetPosition();
186 const G4ThreeVector& posOut = out.GetPosition();
187 const G4ThreeVector momIn(in.GetMomentum().x(), in.GetMomentum().y(),
188 in.GetMomentum().z());
189#if defined(CDC_DEBUG)
190 std::cout << "pid = " << pid << std::endl;
191 std::cout << "mass = " << t.GetDefinition()->GetPDGMass() << std::endl;
192 std::cout << "posIn = " << posIn << std::endl;
193 std::cout << "posOut= " << posOut << std::endl;
194 std::cout << "tof at post-step = " << out.GetGlobalTime() << std::endl;
195 std::cout << "stepl = " << stepLength << std::endl;
196#endif
197
198 // Get layer ID
199 const unsigned layerId = v.GetCopyNo();
200 const unsigned layerIDWithLayerOffset = layerId + m_cdcgp->getOffsetOfFirstLayer();
201 B2DEBUG(150, "LayerID in continuous counting method: " << layerId);
202
203 // If neutral particles, ignore them, unless monopoles.
204
205 if ((charge == 0.) && (abs(pid) != 99666)) return false;
206
207 // Calculate cell ID
208 B2Vector3D tposIn(posIn.x() / CLHEP::cm, posIn.y() / CLHEP::cm, posIn.z() / CLHEP::cm);
209 B2Vector3D tposOut(posOut.x() / CLHEP::cm, posOut.y() / CLHEP::cm, posOut.z() / CLHEP::cm);
210 const unsigned idIn = m_cdcgp->cellId(layerIDWithLayerOffset, tposIn);
211 const unsigned idOut = m_cdcgp->cellId(layerIDWithLayerOffset, tposOut);
212#if defined(CDC_DEBUG)
213 std::cout << "edep= " << edep << std::endl;
214 std::cout << "idIn,idOut= " << idIn << " " << idOut << std::endl;
215#endif
216
217 // Calculate drift length
218 std::vector<int> wires = WireId_in_hit_order(idIn, idOut, m_cdcgp->nWiresInLayer(layerIDWithLayerOffset));
219 G4double sint(0.);
220 const G4double s_in_layer = stepLength / CLHEP::cm;
221 G4double xint[6] = {0};
222
223 const G4ThreeVector momOut(out.GetMomentum().x(), out.GetMomentum().y(),
224 out.GetMomentum().z());
225 const G4double speedIn = in.GetVelocity();
226 const G4double speedOut = out.GetVelocity();
227 const G4double speed = 0.5 * (speedIn + speedOut);
228 const G4double speedInCmPerNs = speed / CLHEP::cm;
229
230 const unsigned int nWires = wires.size();
231 G4double tofBefore = in.GetGlobalTime();
232 G4double kinEnergyBefore = in.GetKineticEnergy();
233 G4double momBefore = momIn.mag();
234 const G4double eLoss = kinEnergyBefore - out.GetKineticEnergy(); //n.b. not always equal to edep
235 const G4double mass = t.GetDefinition()->GetPDGMass();
236#if defined(CDC_DEBUG)
237 std::cout << "momBefore = " << momBefore << std::endl;
238 std::cout << "momIn = " << momIn.x() << " " << momIn.y() << " " << momIn.z() << std::endl;
239 std::cout << "momOut= " << momOut.x() << " " << momOut.y() << " " << momOut.z() << std::endl;
240 std::cout << "speedIn,speedOut= " << speedIn << " " << speedOut << std::endl;
241 std::cout << " speedInCmPerNs= " << speedInCmPerNs << std::endl;
242 std::cout << "tofBefore= " << tofBefore << std::endl;
243#endif
244
245 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
246
247 for (unsigned i = 0; i < nWires; ++i) {
248#if defined(CDC_DEBUG)
249 std::cout << "============ i,wires[i]= " << i << " " << wires[i] << std::endl;
250#endif
251
252 const G4double pos[3] = {posIn.x(), posIn.y(), posIn.z()};
253 G4double Bfield[3];
254 field->GetFieldValue(pos, Bfield);
255 m_magneticField = (Bfield[0] == 0. && Bfield[1] == 0. &&
256 Bfield[2] == 0.) ? false : true;
257#if defined(CDC_DEBUG)
258 std::cout << "Bfield= " << Bfield[0] << " " << Bfield[1] << " " << Bfield[2] << std::endl;
259 std::cout << "magneticField= " << m_magneticField << std::endl;
260#endif
261
262 double distance = 0;
263 G4ThreeVector posW(0, 0, 0);
264 HepPoint3D onTrack;
265 HepPoint3D pOnTrack;
266
267 // Calculate forward/backward position of current wire
268 const B2Vector3D tfw3v = m_cdcgp->wireForwardPosition(layerIDWithLayerOffset, wires[i]);
269 const B2Vector3D tbw3v = m_cdcgp->wireBackwardPosition(layerIDWithLayerOffset, wires[i]);
270
271 const HepPoint3D fwd(tfw3v.x(), tfw3v.y(), tfw3v.z());
272 const HepPoint3D bck(tbw3v.x(), tbw3v.y(), tbw3v.z());
273
274 if (m_magneticField && (abs(pid) != 99666)) {
275 // For monopoles a line segment approximation in the step volume is done,
276 // which is more reasonable, but should be done with a proper catenary FIXME
277 // Cal. distance assuming helix track (still approximation)
278 m_nonUniformField = 1;
279 if (Bfield[0] == 0. && Bfield[1] == 0. &&
280 Bfield[2] != 0.) m_nonUniformField = 0;
281
282 const G4double B_kG[3] = {Bfield[0] / CLHEP::kilogauss,
283 Bfield[1] / CLHEP::kilogauss,
284 Bfield[2] / CLHEP::kilogauss
285 };
286#if defined(CDC_DEBUG)
287 std::cout << "B_kG= " << B_kG[0] << " " << B_kG[1] << " " << B_kG[2] << std::endl;
288 std::cout << "magneticField= " << m_magneticField << std::endl;
289#endif
290
291 const HepPoint3D x(pos[0] / CLHEP::cm, pos[1] / CLHEP::cm, pos[2] / CLHEP::cm);
292 const HepVector3D p(momIn.x() / CLHEP::GeV, momIn.y() / CLHEP::GeV, momIn.z() / CLHEP::GeV);
293 Helix tmp(x, p, charge);
294 tmp.bFieldZ(B_kG[2]);
295 tmp.ignoreErrorMatrix();
296
297 /* // Calculate forward/backward position of current wire
298 const B2Vector3D tfw3v = cdcg.wireForwardPosition(layerId, wires[i]);
299 const B2Vector3D tbw3v = cdcg.wireBackwardPosition(layerId, wires[i]);
300
301 const HepPoint3D fwd(tfw3v.x(), tfw3v.y(), tfw3v.z());
302 const HepPoint3D bck(tbw3v.x(), tbw3v.y(), tbw3v.z());
303 */
304
305 const HepVector3D wire = fwd - bck;
306 HepPoint3D tryp =
307 (x.z() - bck.z()) / wire.z() * wire + bck;
308 tmp.pivot(tryp);
309 tryp = (tmp.x(0.).z() - bck.z()) / wire.z() * wire + bck;
310 tmp.pivot(tryp);
311 tryp = (tmp.x(0.).z() - bck.z()) / wire.z() * wire + bck;
312 tmp.pivot(tryp);
313
314 distance = std::abs(tmp.a()[0]);
315 posW.setX(tryp.x());
316 posW.setY(tryp.y());
317 posW.setZ(tryp.z());
318
319 // HepPoint3D onTrack = tmp.x(0.);
320 onTrack = tmp.x(0.);
321 pOnTrack = tmp.momentum(0.);
322
323 for_Rotat(B_kG);
324 const G4double xwb(bck.x()), ywb(bck.y()), zwb(bck.z());
325 const G4double xwf(fwd.x()), ywf(fwd.y()), zwf(fwd.z());
326 const G4double xp(onTrack.x()), yp(onTrack.y()), zp(onTrack.z());
327 const G4double px(pOnTrack.x()), py(pOnTrack.y()), pz(pOnTrack.z());
328 G4double q2[3] = {0.}, q1[3] = {0.}, q3[3] = {0.};
329 const G4int ntryMax(50); //tentative; too large probably...
330 G4double dist;
331 G4int ntry(999);
332 HELWIR(xwb, ywb, zwb, xwf, ywf, zwf,
333 xp, yp, zp, px, py, pz,
334 B_kG, charge, ntryMax, dist, q2, q1, q3, ntry);
335
336#if defined(CDC_DEBUG)
337 std::cout << "ntry= " << ntry << std::endl;
338 std::cout << "bf distance= " << distance << std::endl;
339 std::cout << "onTrack = " << onTrack << std::endl;
340 std::cout << "posW = " << posW << std::endl;
341#endif
342 if (ntry <= ntryMax) {
343 if (m_wireSag) {
344 G4double ywb_sag, ywf_sag;
345 m_cdcgp->getWireSagEffect(CDCGeometryPar::c_Base, layerIDWithLayerOffset, wires[i], q2[2], ywb_sag, ywf_sag);
346 HELWIR(xwb, ywb_sag, zwb, xwf, ywf_sag, zwf,
347 xp, yp, zp, px, py, pz,
348 B_kG, charge, ntryMax, dist, q2, q1, q3, ntry);
349 }
350 if (ntry <= ntryMax) {
351 distance = dist;
352 onTrack.setX(q1[0]);
353 onTrack.setY(q1[1]);
354 onTrack.setZ(q1[2]);
355 posW.setX(q2[0]);
356 posW.setY(q2[1]);
357 posW.setZ(q2[2]);
358 pOnTrack.setX(q3[0]);
359 pOnTrack.setY(q3[1]);
360 pOnTrack.setZ(q3[2]);
361 }
362#if defined(CDC_DEBUG)
363 std::cout << " " << std::endl;
364 std::cout << "helix distance= " << distance << std::endl;
365 std::cout << "onTrack = " << onTrack << std::endl;
366 std::cout << "posW = " << posW << std::endl;
367 std::cout << "pOnTrack= " << pOnTrack << std::endl;
368 G4ThreeVector bwp(bck.x(), bck.y(), bck.z());
369 G4ThreeVector fwp(fwd.x(), fwd.y(), fwd.z());
370 G4ThreeVector hitPosition, wirePosition;
371 distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
372 hitPosition, wirePosition);
373 if (m_wireSag) {
374 G4double ywb_sag, ywf_sag;
375 m_cdcgp->getWireSagEffect(CDCGeometryPar::c_Base, layerIDWithLayerOffset, wires[i], wirePosition.z(), ywb_sag, ywf_sag);
376 bwp.setY(ywb_sag);
377 fwp.setY(ywf_sag);
378 distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
379 hitPosition, wirePosition);
380 }
381 std::cout << "line distance= " << distance << std::endl;
382 std::cout << "onTrack= " << hitPosition.x() << " " << hitPosition.y() << " " << hitPosition.z() << std::endl;
383 std::cout << "posW = " << wirePosition.x() << " " << wirePosition.y() << " " << wirePosition.z() << std::endl;
384#endif
385 }
386 } else { //no magnetic field case
387 // Cal. distance assuming a line track
388 G4ThreeVector bwp(bck.x(), bck.y(), bck.z());
389 G4ThreeVector fwp(fwd.x(), fwd.y(), fwd.z());
390 G4ThreeVector hitPosition, wirePosition;
391 distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
392 hitPosition, wirePosition);
393 if (m_wireSag) {
394 G4double ywb_sag, ywf_sag;
395 m_cdcgp->getWireSagEffect(CDCGeometryPar::c_Base, layerIDWithLayerOffset, wires[i], wirePosition.z(), ywb_sag, ywf_sag);
396 bwp.setY(ywb_sag);
397 fwp.setY(ywf_sag);
398 distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
399 hitPosition, wirePosition);
400 }
401
402 onTrack.setX(hitPosition.x());
403 onTrack.setY(hitPosition.y());
404 onTrack.setZ(hitPosition.z());
405 posW.setX(wirePosition.x());
406 posW.setY(wirePosition.y());
407 posW.setZ(wirePosition.z());
408 //tentative setting
409 pOnTrack.setX(0.5 * (momIn.x() + momOut.x()) / CLHEP::GeV);
410 pOnTrack.setY(0.5 * (momIn.y() + momOut.y()) / CLHEP::GeV);
411 pOnTrack.setZ(0.5 * (momIn.z() + momOut.z()) / CLHEP::GeV);
412 } //end of magneticfiled on or off
413
414#if defined(CDC_DEBUG)
415 std::cout << "af distance= " << distance << std::endl;
416 std::cout << "onTrack = " << onTrack << std::endl;
417 std::cout << "posW = " << posW << std::endl;
418 std::cout << "pOnTrack = " << pOnTrack << std::endl;
419 if (distance > 2.4) {
420 std::cout << "toolargedriftl" << std::endl;
421 }
422#endif
423 distance *= CLHEP::cm; onTrack *= CLHEP::cm; posW *= CLHEP::cm;
424 pOnTrack *= CLHEP::GeV;
425
426 G4ThreeVector posTrack(onTrack.x(), onTrack.y(), onTrack.z());
427 G4ThreeVector mom(pOnTrack.x(), pOnTrack.y(), pOnTrack.z());
428
429 const B2Vector3D tPosW(posW.x(), posW.y(), posW.z());
430 const B2Vector3D tPosTrack(posTrack.x(), posTrack.y(), posTrack.z());
431 const B2Vector3D tMom(mom.x(), mom.y(), mom.z());
432 G4int lr = m_cdcgp->getOldLeftRight(tPosW, tPosTrack, tMom);
433 G4int newLrRaw = m_cdcgp->getNewLeftRightRaw(tPosW, tPosTrack, tMom);
434 // if(abs(pid) == 11) {
435 // std::cout <<"pid,lr,newLrRaw 4electron= " << pid <<" "<< lr <<" "<< newLrRaw << std::endl;
436 // }
437 G4int newLr = newLrRaw; //to be modified in EndOfEvent
438
439 if (nWires == 1) {
440
441 // saveSimHit(layerId, wires[i], trackID, pid, distance, tofBefore, edep, s_in_layer * cm, momIn, posW, posIn, posOut, posTrack, lr, newLrRaw, newLr, speed);
442 saveSimHit(layerIDWithLayerOffset, wires[i], trackID, pid, distance, tofBefore, edep, s_in_layer * CLHEP::cm, pOnTrack, posW, posIn,
443 posOut,
444 posTrack, lr, newLrRaw, newLr, speed, hitWeight);
445#if defined(CDC_DEBUG)
446 std::cout << "saveSimHit" << std::endl;
447 std::cout << "momIn = " << momIn << std::endl;
448 std::cout << "pOnTrack = " << pOnTrack << std::endl;
449#endif
450
451 } else {
452
453 G4int cel1 = wires[i] + 1;
454 G4int cel2 = cel1;
455 if (i + 1 <= nWires - 1) {
456 cel2 = wires[i + 1] + 1;
457 }
458 const G4double s2 = t.GetTrackLength() / CLHEP::cm; //at post-step
459 G4double s1 = (s2 - s_in_layer); //at pre-step; varied later
460 G4ThreeVector din = momIn;
461 if (din.mag() != 0.) din /= momIn.mag();
462
463 G4double vent[6] = {posIn.x() / CLHEP::cm, posIn.y() / CLHEP::cm, posIn.z() / CLHEP::cm, din.x(), din.y(), din.z()};
464
465 G4ThreeVector dot(momOut.x(), momOut.y(), momOut.z());
466 if (dot.mag() != 0.) {
467 dot /= dot.mag();
468 } else {
469 // Flight-direction is needed to set even when a particle stops
470 dot = din;
471 }
472
473 G4double vext[6] = {posOut.x() / CLHEP::cm, posOut.y() / CLHEP::cm, posOut.z() / CLHEP::cm, dot.x(), dot.y(), dot.z()};
474
475 if (i > 0) {
476 for (int j = 0; j < 6; ++j) vent[j] = xint[j];
477 s1 = sint;
478 }
479
480 // const G4int ic(3); // cubic approximation of the track
481 G4int flag(0);
482 G4double edep_in_cell(0.);
483 G4double eLossInCell(0.);
484
485 if (cel1 != cel2) {
486#if defined(CDC_DEBUG)
487 std::cout << "layerId,cel1,cel2= " << layerId << " " << cel1 << " " << cel2 << std::endl;
488 std::cout << "vent= " << vent[0] << " " << vent[1] << " " << vent[2] << " " << vent[3] << " " << vent[4] << " " << vent[5] <<
489 std::endl;
490 std::cout << "vext= " << vext[0] << " " << vext[1] << " " << vext[2] << " " << vext[3] << " " << vext[4] << " " << vext[5] <<
491 std::endl;
492 std::cout << "s1,s2,ic= " << s1 << " " << s2 << " " << ic << std::endl;
493#endif
494 CellBound(layerIDWithLayerOffset, cel1, cel2, vent, vext, s1, s2, xint, sint, flag);
495#if defined(CDC_DEBUG)
496 std::cout << "flag,sint= " << flag << " " << sint << std::endl;
497 std::cout << "xint= " << xint[0] << " " << xint[1] << " " << xint[2] << " " << xint[3] << " " << xint[4] << " " << xint[5] <<
498 std::endl;
499#endif
500
501 const G4double test = (sint - s1) / s_in_layer;
502 if (test < 0. || test > 1.) {
503 B2WARNING("CDCSensitiveDetector: Strange path length: " << "s1= " << s1 << " sint= " << sint << " s_in_layer= " << s_in_layer <<
504 " test= " << test);
505 }
506 edep_in_cell = edep * std::abs((sint - s1)) / s_in_layer;
507
508 const G4ThreeVector x_In(vent[0]*CLHEP::cm, vent[1]*CLHEP::cm, vent[2]*CLHEP::cm);
509 const G4ThreeVector x_Out(xint[0]*CLHEP::cm, xint[1]*CLHEP::cm, xint[2]*CLHEP::cm);
510 const G4ThreeVector p_In(momBefore * vent[3], momBefore * vent[4], momBefore * vent[5]);
511
512 // saveSimHit(layerId, wires[i], trackID, pid, distance, tofBefore, edep_in_cell, (sint - s1) * cm, p_In, posW, x_In, x_Out, posTrack, lr, newLrRaw, newLr, speed);
513 saveSimHit(layerIDWithLayerOffset, wires[i], trackID, pid, distance, tofBefore, edep_in_cell, std::abs((sint - s1)) * CLHEP::cm,
514 pOnTrack, posW,
515 x_In, x_Out,
516 posTrack, lr, newLrRaw, newLr, speed, hitWeight);
517#if defined(CDC_DEBUG)
518 std::cout << "saveSimHit" << std::endl;
519 std::cout << "p_In = " << p_In << std::endl;
520 std::cout << "pOnTrack= " << pOnTrack << std::endl;
521#endif
522 tofBefore += (sint - s1) / speedInCmPerNs;
523 eLossInCell = eLoss * (sint - s1) / s_in_layer;
524 kinEnergyBefore -= eLossInCell;
525 if (kinEnergyBefore >= 0.) {
526 momBefore = sqrt(kinEnergyBefore * (kinEnergyBefore + 2.*mass));
527 } else {
528 B2WARNING("CDCSensitiveDetector: Kinetic Energy < 0.");
529 momBefore = 0.;
530 }
531
532 } else { //the particle exits
533
534 const G4double test = (s2 - sint) / s_in_layer;
535 if (test < 0. || test > 1.) {
536 B2WARNING("CDCSensitiveDetector: Strange path length: " << "s2= " << s2 << " sint= " << sint << " s_in_layer= " << s_in_layer <<
537 " test= " << test);
538 }
539 edep_in_cell = edep * std::abs((s2 - sint)) / s_in_layer;
540
541 const G4ThreeVector x_In(vent[0]*CLHEP::cm, vent[1]*CLHEP::cm, vent[2]*CLHEP::cm);
542 const G4ThreeVector p_In(momBefore * vent[3], momBefore * vent[4], momBefore * vent[5]);
543
544 // saveSimHit(layerId, wires[i], trackID, pid, distance, tofBefore, edep_in_cell, (s2 - sint) * cm, p_In, posW, x_In, posOut, posTrack, lr, newLrRaw, newLr, speed);
545 saveSimHit(layerIDWithLayerOffset, wires[i], trackID, pid, distance, tofBefore, edep_in_cell, std::abs((s2 - sint)) * CLHEP::cm,
546 pOnTrack, posW,
547 x_In,
548 posOut, posTrack, lr, newLrRaw, newLr, speed, hitWeight);
549#if defined(CDC_DEBUG)
550 std::cout << "saveSimHit" << std::endl;
551 std::cout << "p_In = " << p_In << std::endl;
552 std::cout << "pOnTrack= " << pOnTrack << std::endl;
553#endif
554 }
555 }
556 //setSeenInDetectorFlag(aStep, MCParticle::c_SeenInCDC);
557
561 //StoreArray<Relation> mcPartToSimHits(getRelationCollectionName());
562 //StoreArray<MCParticle> mcPartArray(DEFAULT_MCPARTICLES);
563 //if (saveIndex < 0) {B2FATAL("SimHit wasn't saved despite charge != 0");}
564 //StoreArray<CDCSimHit> m_CDCSimHits(DEFAULT_CDCSIMHITS);
565
566 //new(mcPartToSimHits->AddrAt(saveIndex)) Relation(mcPartArray, m_CDCSimHits, trackID, saveIndex);
567
568 } //end of wire loop
569
570 return true;
571 }
Belle2::CDC::CDCGeometryPar::cellId
unsigned cellId(unsigned layerId, const B2Vector3D &position) const
The method to get cell id based on given layer id and the position.
Definition: CDCGeometryPar.cc:1831
Belle2::CDC::CDCGeometryPar::getWireSagEffect
void getWireSagEffect(EWirePosition set, unsigned layerID, unsigned cellID, double zw, double &ywb_sag, double &ywf_sag) const
Compute effects of the sense wire sag.
Definition: CDCGeometryPar.cc:1960
Belle2::CDC::CDCGeometryPar::getOffsetOfFirstLayer
ushort getOffsetOfFirstLayer() const
Get the offset of the first layer.
Definition: CDCGeometryPar.h:1075
Belle2::CDC::CDCGeometryPar::getOldLeftRight
unsigned short getOldLeftRight(const B2Vector3D &posOnWire, const B2Vector3D &posOnTrack, const B2Vector3D &momentum) const
Returns old left/right.
Definition: CDCGeometryPar.cc:2751
Belle2::CDC::CDCGeometryPar::Instance
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Definition: CDCGeometryPar.cc:32
Belle2::Simulation::UserInfo::getInfo
static Payload getInfo(Carrier &obj)
Static function to just return UserInformation attached to the obj of type Carrier.
Definition: UserInfo.h:100
Belle2::CDC::CDCSensitiveDetector::for_Rotat
void for_Rotat(const G4double bfld[3])
Calculates a rotation matrix.
Definition: CDCSensitiveDetector.cc:971
Belle2::CDC::CDCSensitiveDetector::HELWIR
void HELWIR(const G4double xwb4, const G4double ywb4, const G4double zwb4, const G4double xwf4, const G4double ywf4, const G4double zwf4, const G4double xp, const G4double yp, const G4double zp, const G4double px, const G4double py, const G4double pz, const G4double B_kG[3], const G4double charge, const G4int ntryMax, G4double &distance, G4double q2[3], G4double q1[3], G4double q3[3], G4int &ntry)
Calculate closest points between helix and wire.
Definition: CDCSensitiveDetector.cc:1062
Belle2::CDC::CDCSensitiveDetector::CellBound
void CellBound(const G4int layerId, const G4int ic1, const G4int ic2, const G4double venter[6], const G4double vexit[6], const G4double s1, const G4double s2, G4double xint[6], G4double &sint, G4int &iflag)
Calculate intersection of track with cell boundary.
Definition: CDCSensitiveDetector.cc:689
Belle2::CDC::CDCSensitiveDetector::saveSimHit
void saveSimHit(const G4int layerId, const G4int wireId, const G4int trackID, const G4int pid, const G4double distance, const G4double tof, const G4double edep, const G4double stepLength, const G4ThreeVector &mom, const G4ThreeVector &posW, const G4ThreeVector &posIn, const G4ThreeVector &posOut, const G4ThreeVector &posTrack, const G4int lr, const G4int NewLrRaw, const G4int NewLr, const G4double speed, const G4double hitWeight)
Save CDCSimHit into datastore.
Definition: CDCSensitiveDetector.cc:586
Belle2::CDC::CDCSensitiveDetector::ClosestApproach
G4double ClosestApproach(G4ThreeVector bwp, G4ThreeVector fwp, G4ThreeVector posIn, G4ThreeVector posOut, G4ThreeVector &hitPosition, G4ThreeVector &wirePosition)
Assume line track to calculate distance between track and wire (drift length).
Definition: CDCSensitiveDetector.cc:1412
Belle2::CDC::CDCSensitiveDetector::WireId_in_hit_order
std::vector< int > WireId_in_hit_order(int id0, int id1, int nWires)
Sort wire id.
Definition: CDCSensitiveDetector.cc:1381

◆ WireId_in_hit_order()

std::vector< int > WireId_in_hit_order ( int  id0,
int  id1,
int  nWires 
)
private

Sort wire id.

Definition at line 1381 of file CDCSensitiveDetector.cc.

1382 {
1383 std::vector<int> list;
1384 int i0 = int(id0);
1385 int i1 = int(id1);
1386 if (abs(i0 - i1) * 2 < int(nWires)) {
1387 if (id0 < id1) {
1388 for (int i = id0; i <= id1; ++i)
1389 list.push_back(i);
1390 } else {
1391 for (int i = id0; i >= id1; i--) {
1392 list.push_back(i);
1393 }
1394 }
1395 } else {
1396 if (id0 < id1) {
1397 for (int i = id0; i >= 0; i--)
1398 list.push_back(i);
1399 for (int i = nWires - 1; i >= id1; i--)
1400 list.push_back(i);
1401 } else {
1402 for (int i = id0; i < nWires; ++i)
1403 list.push_back(i);
1404 for (int i = 0; i <= id1; ++i)
1405 list.push_back(i);
1406 }
1407 }
1408
1409 return list;
1410 }