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

Topics

 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  CDCDedxCosLayerAlgorithm
 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  CDCDedxValidationAlgorithm
 A validation algorithm for CDC dE/dx electron. 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  CDCAlphaScaleFactorForAsymmetry
 Database object for scale factors on alpha for CDC hit charge asymmetry. More...
 
class  CDCBadBoards
 Database object for bad boards. 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 1D cell correction 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 cosine 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 1234 of file CDCSensitiveDetector.cc.

1235 {
1236 //require within the same super-layer
1237 if (otherWireId.getISuperLayer() != wireId.getISuperLayer()) return 0;
1238
1239 const signed short iWire = wireId.getIWire();
1240 const signed short iOtherWire = otherWireId.getIWire();
1241 const signed short iCLayer = wireId.getICLayer();
1242 const signed short iOtherCLayer = otherWireId.getICLayer();
1243
1244 //require nearby wire
1245 if (iWire == iOtherWire) {
1246 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1247 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1248 } else {
1249 return 0;
1250 }
1251
1252 signed short iLayerDifference = otherWireId.getILayer() - wireId.getILayer();
1253 if (abs(iLayerDifference) > 1) return 0;
1254
1255 if (iLayerDifference == 0) {
1256 if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer))) return CW_NEIGHBOR;
1257 else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) return CCW_NEIGHBOR;
1258 else return 0;
1259 } else if (iLayerDifference == -1) {
1260 const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1261 m_cdcgp->getShiftInSuperLayer(wireId.getISuperLayer(), wireId.getILayer());
1262 if (iWire == iOtherWire) {
1263 if (deltaShift == CW) return CW_IN_NEIGHBOR;
1264 else if (deltaShift == CCW) return CCW_IN_NEIGHBOR;
1265 else return 0;
1266 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1267 if (deltaShift == CCW) return CW_IN_NEIGHBOR;
1268 else return 0;
1269 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1270 if (deltaShift == CW) return CCW_IN_NEIGHBOR;
1271 else return 0;
1272 } else return 0;
1273 } else if (iLayerDifference == 1) {
1274 const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1275 m_cdcgp->getShiftInSuperLayer(wireId.getISuperLayer(), wireId.getILayer());
1276 if (iWire == iOtherWire) {
1277 if (deltaShift == CW) return CW_OUT_NEIGHBOR;
1278 else if (deltaShift == CCW) return CCW_OUT_NEIGHBOR;
1279 else return 0;
1280 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1281 if (deltaShift == CCW) return CW_OUT_NEIGHBOR;
1282 else return 0;
1283 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1284 if (deltaShift == CW) return CCW_OUT_NEIGHBOR;
1285 else return 0;
1286 } else return 0;
1287 } else return 0;
1288
1289 }

◆ 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 1291 of file CDCSensitiveDetector.cc.

1293 {
1294 //require within the same super-layer
1295 if (otherWireId.getISuperLayer() != iSuperLayer) return 0;
1296
1297 const signed short iOtherWire = otherWireId.getIWire();
1298 const signed short iOtherCLayer = otherWireId.getICLayer();
1299
1300 //require nearby wire
1301 if (iWire == iOtherWire) {
1302 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1303 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1304 } else {
1305 return 0;
1306 }
1307
1308 signed short iLayerDifference = otherWireId.getILayer() - iLayer;
1309 if (abs(iLayerDifference) > 1) return 0;
1310
1311 if (iLayerDifference == 0) {
1312 if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer))) return CW_NEIGHBOR;
1313 else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) return CCW_NEIGHBOR;
1314 else return 0;
1315 } else if (iLayerDifference == -1) {
1316 const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1317 m_cdcgp->getShiftInSuperLayer(iSuperLayer, iLayer);
1318 if (iWire == iOtherWire) {
1319 if (deltaShift == CW) return CW_IN_NEIGHBOR;
1320 else if (deltaShift == CCW) return CCW_IN_NEIGHBOR;
1321 else return 0;
1322 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1323 if (deltaShift == CCW) return CW_IN_NEIGHBOR;
1324 else return 0;
1325 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1326 if (deltaShift == CW) return CCW_IN_NEIGHBOR;
1327 else return 0;
1328 } else return 0;
1329 } else if (iLayerDifference == 1) {
1330 const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1331 m_cdcgp->getShiftInSuperLayer(iSuperLayer, iLayer);
1332 if (iWire == iOtherWire) {
1333 if (deltaShift == CW) return CW_OUT_NEIGHBOR;
1334 else if (deltaShift == CCW) return CCW_OUT_NEIGHBOR;
1335 else return 0;
1336 } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1337 if (deltaShift == CCW) return CW_OUT_NEIGHBOR;
1338 else return 0;
1339 } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1340 if (deltaShift == CW) return CCW_OUT_NEIGHBOR;
1341 else return 0;
1342 } else return 0;
1343 } else return 0;
1344
1345 }

◆ 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(nullptr),
51 m_thresholdEnergyDeposit(thresholdEnergyDeposit),
52 m_thresholdKineticEnergy(thresholdKineticEnergy), m_hitNumber(0)
53 {
54 RelationArray cdcSimHitRel(m_MCParticles, m_CDCSimHits);
55 registerMCParticleRelation(cdcSimHitRel);
56 m_CDCSimHits.registerInDataStore();
57 m_MCParticles.registerRelationTo(m_CDCSimHits);
58
59 const CDCSimControlPar& cntlp = CDCSimControlPar::getInstance();
60
61 m_thresholdEnergyDeposit = cntlp.getThresholdEnergyDeposit();
62 m_thresholdEnergyDeposit *= CLHEP::GeV; //GeV to MeV (=unit in G4)
63 B2DEBUG(150, "CDCSensitiveDetector: Threshold energy (MeV): " << m_thresholdEnergyDeposit);
64 m_thresholdKineticEnergy = 0.0; // Dummy to avoid a warning (tentative).
65
66 m_wireSag = cntlp.getWireSag();
67 B2DEBUG(150, "CDCSensitiveDetector: Sense wire sag on(=1)/off(=0): " << m_wireSag);
68
69 m_modifiedLeftRightFlag = cntlp.getModLeftRightFlag();
70 B2DEBUG(150, "CDCSensitiveDetector: Set left/right flag modified for tracking (=1)/ not set (=0): " << m_modifiedLeftRightFlag);
71
72 m_minTrackLength = cntlp.getMinTrackLength();
73 m_minTrackLength *= CLHEP::cm; //cm to mm (=unit in G4)
74 B2DEBUG(150, "CDCSensitiveDetector: MinTrackLength (mm): " << m_minTrackLength);
75 }

◆ 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 449 of file CDCSensitiveDetector.cc.

457 {
458 //---------------------------------------------------------------------------
459 // (Purpose)
460 // calculate an intersection of track with cell boundary.
461 //
462 // (Relations)
463 // Calls GCUBS
464 //
465 // (Arguments)
466 // input
467 // ic1 serial cell# (start w/ one) of entrance.
468 // ic2 serial cell# (start w/ one) of exit.
469 // venter(6) (x,y,z,px/p,py/p,pz/p) at entrance.
470 // vexit(6) (x,y,z,px/p,py/p,pz/p) at exit.
471 // s1 track length at entrance.
472 // s2 track length at exit.
473 // output
474 // xint(6) (x,y,z,px/p,py/p,pz/p) at intersection of cell boundary.
475 // sint track length at intersection of cell boundary.
476 // iflag return code.
477 //
478 // N.B.(TODO ?) CDC misalignment wrt Belle2 coordinate system is ignored
479 // when calculating the cell-boundary assuming misalign. is small.
480 //--------------------------------------------------------------------------
481
482 G4double div = m_cdcgp->nWiresInLayer(layerId);
483
484 //Check if s1, s2, ic1 and ic2 are ok
485 if (s1 >= s2) {
486 B2ERROR("CDCSensitiveDetector: s1(=" << s1 << ") > s2(=" << s2 << ")");
487 }
488 if (std::abs(ic1 - ic2) != 1) {
489 if (ic1 == 1 && ic2 == div) {
490 } else if (ic1 == div && ic2 == 1) {
491 } else {
492 B2ERROR("CDCSensitiveDetector: |ic1 - ic2| != 1 in CellBound; " << "ic1=" << ic1 << " " << "ic2=" << ic2);
493 }
494 }
495
496 //get wire positions for the entrance cell
497 G4double xwb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).x();
498 G4double ywb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).y();
499 G4double zwb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).z();
500 G4double xwf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).x();
501 G4double ywf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).y();
502 G4double zwf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).z();
503
504 //copy arrays
505 G4double xx1[6], xx2[6];
506 for (int i = 0; i < 6; ++i) {
507 xx1[i] = venter[i];
508 xx2[i] = vexit [i];
509 }
510
511 //calculate the field wire position betw. cell#1 and #2
512 G4double psi = double(ic2 - ic1) * CLHEP::pi / div;
513 if (ic1 == 1 && ic2 == div) {
514 psi = -CLHEP::pi / div;
515 } else if (ic1 == div && ic2 == 1) {
516 psi = CLHEP::pi / div;
517 }
518 G4double cospsi = cos(psi);
519 G4double sinpsi = sin(psi);
520
521 G4double xfwb = cospsi * xwb - sinpsi * ywb;
522 G4double yfwb = sinpsi * xwb + cospsi * ywb;
523 G4double xfwf = cospsi * xwf - sinpsi * ywf;
524 G4double yfwf = sinpsi * xwf + cospsi * ywf;
525 G4double zfwb = zwb;
526 G4double zfwf = zwf;
527
528 //prepare quantities related to the cell-boundary
529 G4double vx = xfwf - xfwb;
530 G4double vy = yfwf - yfwb;
531 G4double vz = zfwf - zfwb;
532 G4double vv = sqrt(vx * vx + vy * vy + vz * vz);
533 vx /= vv; vy /= vv; vz /= vv;
534
535 //translate to make the cubic description easier
536 G4double shiftx = (xx1[0] + xx2[0]) * 0.5;
537 G4double shifty = (xx1[1] + xx2[1]) * 0.5;
538 G4double shiftz = (xx1[2] + xx2[2]) * 0.5;
539 G4double shifts = (s1 + s2) * 0.5;
540 G4double xshft = xx1[0] - shiftx;
541 G4double yshft = xx1[1] - shifty;
542 G4double zshft = xx1[2] - shiftz;
543 G4double sshft = s1 - shifts;
544
545 //approximate the trajectroy by cubic curves
546 G4double pabs1 = sqrt(xx1[3] * xx1[3] + xx1[4] * xx1[4] + xx1[5] * xx1[5]);
547 G4double pabs2 = sqrt(xx2[3] * xx2[3] + xx2[4] * xx2[4] + xx2[5] * xx2[5]);
548
549 G4double a[4] = {0.}, b[4] = {0.}, c[4] = {0.};
550
551 if (m_magneticField) {
552 GCUBS(sshft, xshft, xx1[3] / pabs1, xx2[3] / pabs2, a);
553 GCUBS(sshft, yshft, xx1[4] / pabs1, xx2[4] / pabs2, b);
554 GCUBS(sshft, zshft, xx1[5] / pabs1, xx2[5] / pabs2, c);
555 } else {
556 //n.b. following is really better ?
557 a[1] = xshft / sshft;
558 b[1] = yshft / sshft;
559 c[1] = zshft / sshft;
560 }
561
562 //calculate an int. point betw. the trajectory and the cell-boundary
563 G4double stry(0.), xtry(0.), ytry(0.), ztry(0.);
564 G4double beta(0.), xfw(0.), yfw(0.);
565 G4double sphi(0.), cphi(0.), dphil(0.), dphih(0.);
566 const G4int maxTrials = 100;
567 const G4double eps = 5.e-4;
568 G4double sl = sshft; // negative value
569 G4double sh = -sshft; // positive value
570 G4int i = 0;
571
572 //set initial value (dphil) for the 1st iteration
573 stry = sl;
574 xtry = shiftx + a[0] + stry * (a[1] + stry * (a[2] + stry * a[3]));
575 ytry = shifty + b[0] + stry * (b[1] + stry * (b[2] + stry * b[3]));
576 ztry = shiftz + c[0] + stry * (c[1] + stry * (c[2] + stry * c[3]));
577 beta = (ztry - zfwb) / vz;
578 xfw = xfwb + beta * vx;
579 yfw = yfwb + beta * vy;
580 sphi = (xtry * yfw - ytry * xfw);
581 cphi = (xtry * xfw + ytry * yfw);
582 dphil = atan2(sphi, cphi); //n.b. no need to conv. to dphi...
583
584 iflag = 1;
585
586 while (((sh - sl) > eps) && (i < maxTrials)) {
587 stry = 0.5 * (sl + sh);
588 xtry = shiftx + a[0] + stry * (a[1] + stry * (a[2] + stry * a[3]));
589 ytry = shifty + b[0] + stry * (b[1] + stry * (b[2] + stry * b[3]));
590 ztry = shiftz + c[0] + stry * (c[1] + stry * (c[2] + stry * c[3]));
591 beta = (ztry - zfwb) / vz;
592 xfw = xfwb + beta * vx;
593 yfw = yfwb + beta * vy;
594
595 sphi = (xtry * yfw - ytry * xfw);
596 cphi = (xtry * xfw + ytry * yfw);
597 dphih = atan2(sphi, cphi); //n.b. no need to conv. to dphi...
598
599 if (dphil * dphih > 0.) {
600 sl = stry;
601 } else {
602 sh = stry;
603 }
604 ++i;
605 }
606
607 if (i >= maxTrials - 1) {
608 iflag = 0;
609 B2WARNING("CDCSensitiveDetector: No intersection ?");
610 }
611 sint = stry;
612
613 //get the trajectory at the int. point
614 xint[0] = a[0] + sint * (a[1] + sint * (a[2] + sint * a[3]));
615 xint[1] = b[0] + sint * (b[1] + sint * (b[2] + sint * b[3]));
616 xint[2] = c[0] + sint * (c[1] + sint * (c[2] + sint * c[3]));
617 xint[3] = a[1] + sint * (2. * a[2] + 3. * sint * a[3]);
618 xint[4] = b[1] + sint * (2. * b[2] + 3. * sint * b[3]);
619 xint[5] = c[1] + sint * (2. * c[2] + 3. * sint * c[3]);
620
621 //translate back to the lab. frame
622 xint[0] += shiftx;
623 xint[1] += shifty;
624 xint[2] += shiftz;
625 sint += shifts;
626
627 //re-normalize to one since abs=1 is not guearanteed in the cubic approx.
628 G4double p = sqrt(xint[3] * xint[3] + xint[4] * xint[4] + xint[5] * xint[5]);
629 xint[3] /= p; xint[4] /= p; xint[5] /= p;
630 }
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 1096 of file CDCSensitiveDetector.cc.

1098 {
1099
1100 B2Vector3D tbwp(bwp.x(), bwp.y(), bwp.z());
1101 B2Vector3D tfwp(fwp.x(), fwp.y(), fwp.z());
1102 B2Vector3D tposIn(posIn.x(), posIn.y(), posIn.z());
1103 B2Vector3D tposOut(posOut.x(), posOut.y(), posOut.z());
1104 B2Vector3D thitPosition(0., 0., 0.);
1105 B2Vector3D twirePosition(0., 0., 0.);
1106
1107 G4double distance = CDC::ClosestApproach(tbwp, tfwp, tposIn, tposOut, thitPosition, twirePosition);
1108
1109 hitPosition.setX(thitPosition.x());
1110 hitPosition.setY(thitPosition.y());
1111 hitPosition.setZ(thitPosition.z());
1112
1113 wirePosition.setX(twirePosition.x());
1114 wirePosition.setY(twirePosition.y());
1115 wirePosition.setZ(twirePosition.z());
1116
1117 return distance;
1118 }
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition B2Vector3.h:516

◆ 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 }

◆ 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 378 of file CDCSensitiveDetector.cc.

379 {
380 setModifiedLeftRightFlag();
381 }

◆ 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 667 of file CDCSensitiveDetector.cc.

668 {
669 //Calculates a rotation matrix in advance at a local position in lab.
670 //The rotation is done about the coord. origin; lab.-frame to B-field
671 //frame in which only Bz-comp. is non-zero.
672 //~dead copy of gsim_cdc_for_rotat.F in gsim-cdc for Belle (for tentative use)
673
674 if (m_nonUniformField == 0) return;
675
676 G4double bx, by, bz;
677 bx = bfld[0];
678 by = bfld[1];
679 bz = bfld[2];
680
681 //cal. rotation matrix
682 G4double bxz, bfield;
683 bxz = bx * bx + bz * bz;
684 bfield = bxz + by * by;
685 bxz = sqrt(bxz);
686 bfield = sqrt(bfield);
687
688 m_brot[0][0] = bz / bxz;
689 m_brot[1][0] = 0.;
690 m_brot[2][0] = -bx / bxz;
691 m_brot[0][1] = -by * bx / bxz / bfield;
692 m_brot[1][1] = bxz / bfield;
693 m_brot[2][1] = -by * bz / bxz / bfield;
694 m_brot[0][2] = bx / bfield;
695 m_brot[1][2] = by / bfield;
696 m_brot[2][2] = bz / bfield;
697
698 return;
699
700 }

◆ 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 632 of file CDCSensitiveDetector.cc.

633 {
634 //Original: GCUBS in Geant3
635 // ******************************************************************
636 // * *
637 // * Calculates a cubic through P1,(X,Y),(-X,-Y),P2 *
638 // * Y=A(1)+A(2)*X+A(3)*X**2+A(4)*X**3 *
639 // * The coordinate system is assumed to be the cms system *
640 // * of P1,P2. *
641 // * d1(2): directional cosine at P1(2). *
642 // * *
643 // * ==>Called by : GIPLAN,GICYL *
644 // * Author H.Boerner ********* *
645 // * *
646 // ******************************************************************
647
648 G4double fact(0);
649
650 if (x == 0.) goto L10;
651
652 fact = (d1 - d2) * 0.25;
653 a[0] = - 1. * fact * x;
654 a[2] = fact / x;
655 a[1] = (6. * y - (d1 + d2) * x) / (4. * x);
656 a[3] = ((d1 + d2) * x - 2.*y) / (4.*x * x * x);
657 return;
658
659L10:
660 a[0] = 0.;
661 a[1] = 1.;
662 a[2] = 0.;
663 a[3] = 0.;
664 }

◆ getDedxFromParticle()

CDCDedxTrack const * getDedxFromParticle ( Particle const * particle)

CDC dEdx value from particle.

Definition at line 34 of file DedxVariables.cc.

35 {
36 const Track* track = particle->getTrack();
37 if (!track) {
38 return nullptr;
39 }
40
41 const CDCDedxTrack* dedxTrack = track->getRelatedTo<CDCDedxTrack>();
42 if (!dedxTrack) {
43 return nullptr;
44 }
45
46 return dedxTrack;
47 }
Belle2::Track
Class that bundles various TrackFitResults.
Definition Track.h:25

◆ 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 }

◆ 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 41 of file CDCDedxSigmaPars.cc.

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

◆ 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 756 of file CDCSensitiveDetector.cc.

770 {
771 //~dead copy of gsim_cdc_hit.F in gsim-cdc for Belle (for tentative use)
772 // ---------------------------------------------------------------------
773 // Purpose : Calculate closest points between helix and wire.
774 //
775 // Input
776 // xwb4 : x of wire at backward endplate in lab.
777 // ywb4 : y of wire at backward endplate "
778 // zwb4 : z of wire at backward endplate "
779 // xwf4 : x of wire at forward endplate "
780 // ywf4 : y of wire at forward endplate "
781 // zwf4 : z of wire at forward endplate "
782 //
783 // Output
784 // q2(1) : x of wire at closest point in lab.
785 // q2(2) : y of wire at closest point "
786 // q2(3) : z of wire at closest point "
787 // q1(1) : x of helix at closest point "
788 // q1(2) : y of helix at closest point "
789 // q1(3) : z of helix at closest point "
790 // ntry :
791 // ---------------------------------------------------------------------
792
793 const G4int ndim = 3;
794 const G4double delta = 1.e-5;
795
796
797 G4double xwb, ywb, zwb, xwf, ywf, zwf;
798 G4double xw, yw, zw, xh, yh, zh, pxh, pyh, pzh;
799 G4double fi, fi_corr;
800
801 G4double dr, fi0, cpa, dz, tanl;
802 G4double x0, y0, z0;
803 G4double xc, yc, r;
804 G4double xwm, ywm;
805 G4double sinfi0, cosfi0, sinfi0fi, cosfi0fi;
806
807 G4double vx, vy, vz, vv, cx, cy, cz, tt[3][3];
808 G4double tmp[3];
809
810 G4double xx[3], dxx[3], ddxx[3], pp[3];
811 G4double xxtdxx, dxxtdxx, xxtddxx;
812
813
814 G4double fst = 0.0;
815 G4double f, fderiv, deltafi, fact, eval;
816 G4double dx1, dy1, dx2, dy2, crs, dot;
817
818 G4int iflg;
819
820 //set parameters
821 xwb = xwb4; ywb = ywb4; zwb = zwb4;
822 xwf = xwf4; ywf = ywf4; zwf = zwf4;
823
824 G4double xxx(xp), yyy(yp), zzz(zp);
825 G4double pxx(px), pyy(py), pzz(pz);
826
827 //rotate z-axis to be parallel to B-field in case of non-uniform B
828 Rotat(xwb, ywb, zwb, 1);
829 Rotat(xwf, ywf, zwf, 1);
830 Rotat(xxx, yyy, zzz, 1);
831 Rotat(pxx, pyy, pzz, 1);
832
833 G4double a[8] = {0.};
834 G4double pt = sqrt(pxx * pxx + pyy * pyy);
835 a[1] = atan2(-pxx, pyy);
836 a[2] = charge / pt;
837 a[4] = pzz / pt;
838 a[5] = xxx; a[6] = yyy; a[7] = zzz;
839
840 //calculate unit direction vector of the sense wire
841 vx = xwf - xwb; vy = ywf - ywb; vz = zwf - zwb;
842 vv = sqrt(vx * vx + vy * vy + vz * vz);
843 vx /= vv; vy /= vv; vz /= vv;
844
845 //flag for distinguishing between stereo and axial wire
846 iflg = 0;
847 if (vx == 0. && vy == 0.) iflg = 1;
848
849 //calculate coefficients of f
850 cx = xwb - vx * (vx * xwb + vy * ywb + vz * zwb);
851 cy = ywb - vy * (vx * xwb + vy * ywb + vz * zwb);
852 cz = zwb - vz * (vx * xwb + vy * ywb + vz * zwb);
853
854 //calculate tensor for f
855 tt[0][0] = vx * vx - 1.; tt[1][0] = vx * vy; tt[2][0] = vx * vz;
856 tt[0][1] = vy * vx; tt[1][1] = vy * vy - 1.; tt[2][1] = vy * vz;
857 tt[0][2] = vz * vx; tt[1][2] = vz * vy; tt[2][2] = vz * vz - 1.;
858
859 //set helix parameters
860 dr = a[0]; fi0 = a[1]; cpa = a[2];
861 dz = a[3]; tanl = a[4];
862 x0 = a[5]; y0 = a[6]; z0 = a[7];
863
864 //
865 // set initial value for phi
866 //
867 xwm = xxx;
868 ywm = yyy;
869
870 G4double bfield = sqrt(B_kG[0] * B_kG[0] +
871 B_kG[1] * B_kG[1] +
872 B_kG[2] * B_kG[2]);
873 G4double alpha = 1.e4 / 2.99792458 / bfield;
874 r = alpha / cpa;
875 cosfi0 = cos(fi0);
876 sinfi0 = sin(fi0);
877 xc = x0 + (dr + r) * cosfi0;
878 yc = y0 + (dr + r) * sinfi0;
879 dx1 = x0 - xc;
880 dy1 = y0 - yc;
881 dx2 = xwm - xc;
882 dy2 = ywm - yc;
883 crs = dx1 * dy2 - dy1 * dx2;
884 dot = dx1 * dx2 + dy1 * dy2;
885 fi = atan2(crs, dot);
886
887 //begin iterative procedure for newton 's method '
888 fact = 1.;
889 ntry = 0;
890line1:
891 ntry += 1;
892 cosfi0fi = cos(fi0 + fi);
893 sinfi0fi = sin(fi0 + fi);
894
895 //calculate spatial point Q(x,y,z) along the helix
896 xx[0] = x0 + dr * cosfi0 + r * (cosfi0 - cosfi0fi);
897 xx[1] = y0 + dr * sinfi0 + r * (sinfi0 - sinfi0fi);
898 xx[2] = z0 + dz - r * tanl * fi;
899 pp[0] = -pt * sinfi0fi;
900 pp[1] = pt * cosfi0fi;
901 pp[2] = pt * tanl;
902
903 if (iflg == 1) {
904 q2[0] = xwb; q2[1] = ywb; q2[2] = xx[2];
905 q1[0] = xx[0]; q1[1] = xx[1]; q1[2] = xx[2];
906 q3[0] = pp[0]; q3[1] = pp[1]; q3[2] = pp[2];
907 //inverse rotation to lab. frame in case of non-uniform B
908 Rotat(q1, -1);
909 Rotat(q2, -1);
910 Rotat(q3, -1);
911 distance = sqrt((q2[0] - q1[0]) * (q2[0] - q1[0]) +
912 (q2[1] - q1[1]) * (q2[1] - q1[1]) +
913 (q2[2] - q1[2]) * (q2[2] - q1[2]));
914 return;
915 }
916
917 //calculate direction vector (dx/dphi,dy/dphi,dz/dphi)
918 //on a point along the helix.
919 dxx[0] = r * sinfi0fi; dxx[1] = - r * cosfi0fi; dxx[2] = - r * tanl;
920
921 // In order to derive the closest point between straight line and helix,
922 // we can put following two conditions:
923 // (i) A point H(xh,yh,zh) on the helix given should be on
924 // the plane which is perpendicular to the straight line.
925 // (ii) A line HW from W(xw,yw,zw) which is a point on the straight
926 // line to H(xh,yh,zh) should normal to the direction vector
927 // on the point H.
928 //
929 // Thus, we can make a equation from above conditions.
930 // f(phi) = cx*(dx/dphi) + cy*(dy/dphi) + cz*(dz/dphi)
931 // + (x,y,z)*tt(i,j)*(dx/dphi,dy/dphi,dz/dphi)
932 // = 0,
933 // where
934 // cx = xwb - vx*( vx*xwb + vy*ywb + vz*zwb )
935 // cy = ywb - vy*( vx*xwb + vy*ywb + vz*zwb )
936 // cz = zwb - vz*( vx*xwb + vy*ywb + vz*zwb )
937 //
938 // tt(1,1) = vx*vx - 1 tt(1,2) = vx*vy tt(1,3) = vx*vz
939 // tt(2,1) = vy*vx tt(2,2) = vy*vy - 1 tt(2,3) = vy*vz
940 // tt(3,1) = vz*vx tt(3,2) = vz*vy tt(3,3) = vz*vz - 1
941 //
942 // and the equation of straight line(stereo wire) is written by
943 // (x,y,z) = (xwb,ywb,zwb) + beta*(vx,vy,vz), beta is free parameter.
944
945 //Now calculate f
946 Mvopr(ndim, xx, tt, dxx, tmp, 1);
947 xxtdxx = tmp[0];
948 f = cx * dxx[0] + cy * dxx[1] + cz * dxx[2] + xxtdxx;
949 if (std::abs(f) < delta) goto line100;
950
951 //evaluate fitting result and prepare some factor to multiply to 1/derivative
952 if (ntry > 1) {
953 eval = (1.0 - 0.25 * fact) * std::abs(fst) - std::abs(f);
954 if (eval <= 0.) fact *= 0.5;
955 }
956
957 //calculate derivative of f
958 ddxx[0] = r * cosfi0fi; ddxx[1] = r * sinfi0fi; ddxx[2] = 0.;
959
960 //Now we have derivative of f
961 Mvopr(ndim, dxx, tt, dxx, tmp, 1);
962 dxxtdxx = tmp[0];
963 Mvopr(ndim, xx, tt, ddxx, tmp, 1);
964 xxtddxx = tmp[0];
965 fderiv = cx * ddxx[0] + cy * ddxx[1] + cz * ddxx[2] + dxxtdxx + xxtddxx;
966 // Commented by M. U. June, 2nd, 2013
967 // fist = fi;
968 deltafi = f / fderiv;
969 fi -= fact * deltafi;
970 fst = f;
971
972 if (ntry > ntryMax) {
973 //B2DEBUG(" Exceed max. trials HelWir ");
974 goto line100;
975 }
976 //write(6,*) ntry, fist, deltafi
977 goto line1;
978
979 //check if zh is btw zwb and zwf; if not, set zh=zwb or zh=zwf.
980 //dead regions due to feed-throughs should be considered later.
981line100:
982 zh = z0 + dz - r * tanl * fi;
983 fi_corr = 0.;
984 if (zh < zwb) fi_corr = (zwb - zh) / (-r * tanl);
985 if (zh > zwf) fi_corr = (zwf - zh) / (-r * tanl);
986 fi += fi_corr;
987
988 cosfi0fi = cos(fi0 + fi);
989 sinfi0fi = sin(fi0 + fi);
990
991 xh = x0 + dr * cosfi0 + r * (cosfi0 - cosfi0fi);
992 yh = y0 + dr * sinfi0 + r * (sinfi0 - sinfi0fi);
993 zh = z0 + dz - r * tanl * fi;
994 pxh = -pt * sinfi0fi;
995 pyh = pt * cosfi0fi;
996 pzh = pt * tanl;
997
998 zw = vx * vz * xh + vy * vz * yh + vz * vz * zh + zwb - vz * (vx * xwb + vy * ywb + vz * zwb);
999 xw = xwb + vx * (zw - zwb) / vz;
1000 yw = ywb + vy * (zw - zwb) / vz;
1001
1002 q2[0] = xw; q2[1] = yw; q2[2] = zw;
1003 q1[0] = xh; q1[1] = yh; q1[2] = zh;
1004 q3[0] = pxh; q3[1] = pyh; q3[2] = pzh;
1005
1006 //inverse rotation to lab. frame in case of non-uniform B
1007 Rotat(q1, -1);
1008 Rotat(q2, -1);
1009 Rotat(q3, -1);
1010 distance = sqrt((q2[0] - q1[0]) * (q2[0] - q1[0]) +
1011 (q2[1] - q1[1]) * (q2[1] - q1[1]) +
1012 (q2[2] - q1[2]) * (q2[2] - q1[2]));
1013 return;
1014
1015 }
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

◆ 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 77 of file CDCSensitiveDetector.cc.

78 {
79 // Initialize
80 m_nonUniformField = 0;
81 }

◆ 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 1018 of file CDCSensitiveDetector.cc.

1020 {
1021 //~dead copy of UtilCDC_mvopr in com-cdc for Belle (for tentative use)
1022 //-----------------------------------------------------------------------
1023 // Input
1024 // ndim : dimension
1025 // b(1-ndim) : vector
1026 // m(1-ndim,1-ndim) : matrix
1027 // a(1-ndim) : vector
1028 // c(1-ndim) : vector
1029 // mode : c = m * a for mode=0
1030 // c = b * m * a for mode=1
1031 // Output
1032 // c(1-ndim) : for mode 1, solution is put on c[0]
1033 //-----------------------------------------------------------------------
1034
1035 if (ndim != 3) {
1036 return;
1037 }
1038
1039 for (int i = 0; i < ndim; ++i) c[i] = 0.;
1040 G4double tmp[3];
1041 for (int i = 0; i < ndim; ++i) tmp[i] = 0.;
1042
1043 if (mode == 0) {
1044 for (int i = 0; i < ndim; ++i) {
1045 for (int j = 0; j < ndim; ++j) {
1046 c[i] += m[j][i] * a[j];
1047 }
1048 }
1049 return;
1050 } else if (mode == 1) {
1051 for (int i = 0; i < ndim; ++i) {
1052 for (int j = 0; j < ndim; ++j) {
1053 tmp[i] += m[j][i] * a[j];
1054 }
1055 c[0] += b[i] * tmp[i];
1056 }
1057 } else {
1058 }
1059
1060 return;
1061
1062 }

◆ reAssignLeftRightInfo()

void reAssignLeftRightInfo ( )
private

Re-assign left/right info.

Definition at line 1184 of file CDCSensitiveDetector.cc.

1185 {
1186 CDCSimHit* sHit = nullptr;
1187 WireID sWireId; // = WireID();
1188 B2Vector3D sPos; // = B2Vector3D();
1189
1190 CDCSimHit* pHit = nullptr;
1191 WireID pWireId; // = WireID();
1192
1193 //Find a primary track close to the input 2'ndary hit in question
1194 for (std::vector<CDCSimHit*>::iterator nIt = m_hitWithNegWeight.begin(), nItEnd = m_hitWithNegWeight.end(); nIt != nItEnd; ++nIt) {
1195
1196 sHit = *nIt;
1197 sPos = sHit->getPosTrack();
1198 sWireId = sHit->getWireID();
1199 unsigned short sClayer = sWireId.getICLayer();
1200 unsigned short sSuperLayer = sWireId.getISuperLayer();
1201 unsigned short sLayer = sWireId.getILayer();
1202 unsigned short sWire = sWireId.getIWire();
1203 CDCSimHit* fHit = sHit;
1204
1205 std::multimap<unsigned short, CDCSimHit*>::iterator pItBegin = m_hitWithPosWeight.find(sSuperLayer);
1206 std::multimap<unsigned short, CDCSimHit*>::iterator pItEnd = m_hitWithPosWeight.find(sSuperLayer + 1);
1207
1208 double minDistance2 = DBL_MAX;
1209
1210 for (std::multimap<unsigned short, CDCSimHit*>::iterator pIt = pItBegin; pIt != pItEnd; ++pIt) {
1211
1212 // scan hits in the same/neighboring cells
1213 pHit = pIt->second;
1214 pWireId = pHit->getWireID();
1215 unsigned short neighb = areNeighbors(sClayer, sSuperLayer, sLayer, sWire, pWireId);
1216 if (neighb != 0 || pWireId == sWireId) {
1217 double distance2 = (pHit->getPosTrack() - sPos).Mag2();
1218 if (distance2 < minDistance2) {
1219 fHit = pHit;
1220 minDistance2 = distance2;
1221 }
1222 }
1223 }
1224
1225 //reassign LR using the momentum-direction of the primary particle found
1226 unsigned short lR = m_cdcgp->getNewLeftRightRaw(sHit->getPosWire(),
1227 sHit->getPosTrack(),
1228 fHit->getMomentum());
1229 sHit->setLeftRightPassage(lR);
1230 }
1231 }

◆ 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 703 of file CDCSensitiveDetector.cc.

705 {
706 //Translates (x,y,z) in lab. to (x,y,z) in B-field frame (mode=1), or reverse
707 // translation (mode=-1).
708 //~dead copy (for tentative use) of gsim_cdc_rotat/irotat.F in gsim-cdc
709 //for Belle
710
711 if (m_nonUniformField == 0) return;
712
713 G4double x0(x), y0(y), z0(z);
714
715 if (mode == 1) {
716 x = m_brot[0][0] * x0 + m_brot[1][0] * y0 + m_brot[2][0] * z0;
717 y = m_brot[0][1] * x0 + m_brot[1][1] * y0 + m_brot[2][1] * z0;
718 z = m_brot[0][2] * x0 + m_brot[1][2] * y0 + m_brot[2][2] * z0;
719 } else if (mode == -1) {
720 x = m_brot[0][0] * x0 + m_brot[0][1] * y0 + m_brot[0][2] * z0;
721 y = m_brot[1][0] * x0 + m_brot[1][1] * y0 + m_brot[1][2] * z0;
722 z = m_brot[2][0] * x0 + m_brot[2][1] * y0 + m_brot[2][2] * z0;
723 } else {
724 }
725 return;
726
727 }

◆ Rotat() [2/2]

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

Overloaded translation method.

Definition at line 730 of file CDCSensitiveDetector.cc.

731 {
732 //Translates (x,y,z) in lab. to (x,y,z) in B-field frame (mode=1), or reverse
733 // translation (mode=-1).
734 //~dead copy (for tentative use) of gsim_cdc_rotat/irotat.F in gsim-cdc
735 //for Belle
736
737 if (m_nonUniformField == 0) return;
738
739 G4double x0(x[0]), y0(x[1]), z0(x[2]);
740
741 if (mode == 1) {
742 x[0] = m_brot[0][0] * x0 + m_brot[1][0] * y0 + m_brot[2][0] * z0;
743 x[1] = m_brot[0][1] * x0 + m_brot[1][1] * y0 + m_brot[2][1] * z0;
744 x[2] = m_brot[0][2] * x0 + m_brot[1][2] * y0 + m_brot[2][2] * z0;
745 } else if (mode == -1) {
746 x[0] = m_brot[0][0] * x0 + m_brot[0][1] * y0 + m_brot[0][2] * z0;
747 x[1] = m_brot[1][0] * x0 + m_brot[1][1] * y0 + m_brot[1][2] * z0;
748 x[2] = m_brot[2][0] * x0 + m_brot[2][1] * y0 + m_brot[2][2] * z0;
749 } else {
750 }
751 return;
752
753 }

◆ 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 384 of file CDCSensitiveDetector.cc.

402 {
403
404 // Discard the hit below Edep_th
405 if (edep <= m_thresholdEnergyDeposit) return;
406
407 //compute tof at the closest point; linear approx.
408 const G4double sign = (posTrack - posIn).dot(mom) < 0. ? -1. : 1.;
409 const G4double CorrectTof = tof + sign * (posTrack - posIn).mag() / speed;
410
411 RelationArray cdcSimHitRel(m_MCParticles, m_CDCSimHits);
412
413 m_hitNumber = m_CDCSimHits.getEntries();
414
415 CDCSimHit* simHit = m_CDCSimHits.appendNew();
416
417 simHit->setWireID(layerId, wireId);
418 simHit->setTrackId(trackID);
419 simHit->setPDGCode(pid);
420 simHit->setDriftLength(distance / CLHEP::cm);
421 simHit->setFlightTime(CorrectTof / CLHEP::ns);
422 simHit->setGlobalTime(CorrectTof / CLHEP::ns);
423 simHit->setEnergyDep(edep / CLHEP::GeV);
424 simHit->setStepLength(stepLength / CLHEP::cm);
425 B2Vector3D momentum(mom.getX() / CLHEP::GeV, mom.getY() / CLHEP::GeV, mom.getZ() / CLHEP::GeV);
426 simHit->setMomentum(momentum);
427 B2Vector3D posWire(posW.getX() / CLHEP::cm, posW.getY() / CLHEP::cm, posW.getZ() / CLHEP::cm);
428 simHit->setPosWire(posWire);
429 B2Vector3D positionIn(posIn.getX() / CLHEP::cm, posIn.getY() / CLHEP::cm, posIn.getZ() / CLHEP::cm);
430 simHit->setPosIn(positionIn);
431 B2Vector3D positionOut(posOut.getX() / CLHEP::cm, posOut.getY() / CLHEP::cm, posOut.getZ() / CLHEP::cm);
432 simHit->setPosOut(positionOut);
433 B2Vector3D positionTrack(posTrack.getX() / CLHEP::cm, posTrack.getY() / CLHEP::cm, posTrack.getZ() / CLHEP::cm);
434 simHit->setPosTrack(positionTrack);
435 simHit->setPosFlag(lr);
436 simHit->setLeftRightPassageRaw(newLrRaw);
437 simHit->setLeftRightPassage(newLr);
438
439 B2DEBUG(150, "HitNumber: " << m_hitNumber);
440 if (m_modifiedLeftRightFlag) {
441 //N.B. Negative hitWeight is allowed intentionally here; all weights are to be reset to positive in EndOfEvent
442 cdcSimHitRel.add(trackID, m_hitNumber, hitWeight);
443 } else {
444 cdcSimHitRel.add(trackID, m_hitNumber);
445 }
446 }

◆ setModifiedLeftRightFlag()

void setModifiedLeftRightFlag ( )
private

set left/right flag modified for tracking

Definition at line 1121 of file CDCSensitiveDetector.cc.

1122 {
1123 if (!m_modifiedLeftRightFlag) return;
1124
1125 // Get SimHit array and relation betw. MC and SimHit
1126 // N.B. MCParticle is incomplete at this stage; the relation betw it and
1127 // simHit is Okay.
1128 // MCParticle will be completed after all sub-detectors' EndOfEvent calls.
1129 RelationArray mcPartToSimHits(m_MCParticles, m_CDCSimHits);
1130 int nRelationsMinusOne = mcPartToSimHits.getEntries() - 1;
1131
1132 if (nRelationsMinusOne == -1) return;
1133
1134 //reset some of negative weights to positive; this is needed for the hits
1135 //created by secondary particles whose track-lengths get larger than the
1136 //threshold (set by the user) during G4 swimming (i.e. the weights are
1137 //first set to negative as far as the track-lengths are shorther than the
1138 //threshold; set to positive when the track-lengths exceed the threshold).
1139
1140 size_t iRelation = 0;
1141 int trackIdOld = INT_MAX;
1142 m_hitWithPosWeight.clear();
1143 m_hitWithNegWeight.clear();
1144
1145 for (int it = nRelationsMinusOne; it >= 0; --it) {
1146 RelationElement& mcPartToSimHit = const_cast<RelationElement&>(mcPartToSimHits[it]);
1147 size_t nRelatedHits = mcPartToSimHit.getSize();
1148 if (nRelatedHits > 1) B2FATAL("CDCSensitiveDetector::EndOfEvent: MCParticle<-> CDCSimHit relation is not one-to-one !");
1149
1150 unsigned short trackId = mcPartToSimHit.getFromIndex();
1151 RelationElement::weight_type weight = mcPartToSimHit.getWeight(iRelation);
1152 if (weight > 0.) {
1153 trackIdOld = trackId;
1154 } else if (weight <= 0. && trackId == trackIdOld) {
1155 weight *= -1.;
1156 mcPartToSimHit.setToIndex(mcPartToSimHit.getToIndex(iRelation), weight);
1157 trackIdOld = trackId;
1158 }
1159
1160 CDCSimHit* sHit = m_CDCSimHits[mcPartToSimHit.getToIndex(iRelation)];
1161
1162 if (weight > 0.) {
1163 m_hitWithPosWeight.insert(std::pair<unsigned short, CDCSimHit*>(sHit->getWireID().getISuperLayer(), sHit));
1164 } else {
1165 m_hitWithNegWeight.push_back(sHit);
1166 }
1167 }
1168
1169 //reassign L/R flag
1170 reAssignLeftRightInfo();
1171
1172 //reset all weights positive; this is required for completing MCParticle object at the EndOfEvent action of FullSim
1173 // is this part really needed ??? check again !
1174 for (int it = 0; it <= nRelationsMinusOne; ++it) {
1175 RelationElement& mcPartToSimHit = const_cast<RelationElement&>(mcPartToSimHits[it]);
1176 RelationElement::weight_type weight = mcPartToSimHit.getWeight(iRelation);
1177 if (weight < 0.) {
1178 mcPartToSimHit.setToIndex(mcPartToSimHit.getToIndex(iRelation), -1.*weight);
1179 }
1180 }
1181
1182 }

◆ sigmaCurve()

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

parameterized resolution for predictions

Definition at line 20 of file CDCDedxSigmaPars.cc.

21 {
22 // calculate the predicted mean value as a function of beta-gamma (bg)
23 // this is done with a different function depending dE/dx, nhit, and sin(theta)
24 double f = 0;
25
26 if (version == 0) {
27 if (par[0] == 1) { // return dedx parameterization
28 f = par[1] + par[2] * x;
29 } else if (par[0] == 2) { // return nhit or sin(theta) parameterization
30 f = par[1] * x * x * x * x + par[2] * x * x * x + par[3] * x * x + par[4] * x + par[5];
31 } else if (par[0] == 3) { // return cos(theta) parameterization
32 f = par[1] * std::exp(-0.5 * square(((x - par[2]) / par[3]))) +
33 par[4] * (pow5(x) * x) + par[5] * pow5(x) + par[6] * pow4(x) +
34 par[7] * x * x * x + par[8] * x * x + par[9] * x + par[10];
35 }
36 }
37
38 return f;
39 }
Belle2::square
constexpr T square(const T &x)
Calculate the square of the input.
Definition MathHelpers.h:21
Belle2::pow5
constexpr T pow5(const T &x)
Calculate the fifth power of the input.
Definition MathHelpers.h:46
Belle2::pow4
constexpr T pow4(const T &x)
Calculate the fourth power of the input.
Definition MathHelpers.h:37

◆ step()

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

Process each step and calculate variables defined in CDCB4VHit.

Implements SensitiveDetectorBase.

Definition at line 86 of file CDCSensitiveDetector.cc.

87 {
88 m_cdcgp = &CDCGeometryPar::Instance();
89 m_nonUniformField = 0;
90
91 // Get deposited energy
92 const G4double edep = aStep->GetTotalEnergyDeposit();
93
94 // Discard the hit below Edep_th
95 if (edep <= m_thresholdEnergyDeposit) return false;
96
97 // Get step length
98 const G4double stepLength = aStep->GetStepLength();
99 if (stepLength == 0.) return false;
100
101 // Get step information
102 const G4Track& t = * aStep->GetTrack();
103
104 G4double hitWeight = Simulation::TrackInfo::getInfo(t).getIgnore() ? -1 : 1;
105 // save in MCParticle if track-length is enough long
106 if (t.GetTrackLength() > m_minTrackLength) {
107 Simulation::TrackInfo::getInfo(t).setIgnore(false);
108 hitWeight = 1.;
109 }
110
111 const G4int pid = t.GetDefinition()->GetPDGEncoding();
112 const G4double charge = t.GetDefinition()->GetPDGCharge();
113 const G4int trackID = t.GetTrackID();
114 const G4VPhysicalVolume& v = * t.GetVolume();
115 const G4StepPoint& in = * aStep->GetPreStepPoint();
116 const G4StepPoint& out = * aStep->GetPostStepPoint();
117 const G4ThreeVector& posIn = in.GetPosition();
118 const G4ThreeVector& posOut = out.GetPosition();
119 const G4ThreeVector momIn(in.GetMomentum().x(), in.GetMomentum().y(),
120 in.GetMomentum().z());
121
122 // Get layer ID
123 const unsigned layerId = v.GetCopyNo();
124 const unsigned layerIDWithLayerOffset = layerId + m_cdcgp->getOffsetOfFirstLayer();
125 B2DEBUG(150, "LayerID in continuous counting method: " << layerId);
126
127 // If neutral particles, ignore them, unless monopoles.
128 if ((charge == 0.) && (abs(pid) != 99666)) return false;
129
130 // Calculate cell ID
131 B2Vector3D tposIn(posIn.x() / CLHEP::cm, posIn.y() / CLHEP::cm, posIn.z() / CLHEP::cm);
132 B2Vector3D tposOut(posOut.x() / CLHEP::cm, posOut.y() / CLHEP::cm, posOut.z() / CLHEP::cm);
133 const unsigned idIn = m_cdcgp->cellId(layerIDWithLayerOffset, tposIn);
134 const unsigned idOut = m_cdcgp->cellId(layerIDWithLayerOffset, tposOut);
135
136 // Calculate drift length
137 std::vector<int> wires = WireId_in_hit_order(idIn, idOut, m_cdcgp->nWiresInLayer(layerIDWithLayerOffset));
138 G4double sint(0.);
139 const G4double s_in_layer = stepLength / CLHEP::cm;
140 G4double xint[6] = {0};
141
142 const G4ThreeVector momOut(out.GetMomentum().x(), out.GetMomentum().y(),
143 out.GetMomentum().z());
144 const G4double speedIn = in.GetVelocity();
145 const G4double speedOut = out.GetVelocity();
146 const G4double speed = 0.5 * (speedIn + speedOut);
147 const G4double speedInCmPerNs = speed / CLHEP::cm;
148
149 const unsigned int nWires = wires.size();
150 G4double tofBefore = in.GetGlobalTime();
151 G4double kinEnergyBefore = in.GetKineticEnergy();
152 G4double momBefore = momIn.mag();
153 const G4double eLoss = kinEnergyBefore - out.GetKineticEnergy(); //n.b. not always equal to edep
154 const G4double mass = t.GetDefinition()->GetPDGMass();
155
156 const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
157
158 for (unsigned i = 0; i < nWires; ++i) {
159
160 const G4double pos[3] = {posIn.x(), posIn.y(), posIn.z()};
161 G4double Bfield[3];
162 field->GetFieldValue(pos, Bfield);
163 m_magneticField = (Bfield[0] == 0. && Bfield[1] == 0. &&
164 Bfield[2] == 0.) ? false : true;
165
166 double distance = 0;
167 G4ThreeVector posW(0, 0, 0);
168 HepPoint3D onTrack;
169 HepPoint3D pOnTrack;
170
171 // Calculate forward/backward position of current wire
172 const B2Vector3D tfw3v = m_cdcgp->wireForwardPosition(layerIDWithLayerOffset, wires[i]);
173 const B2Vector3D tbw3v = m_cdcgp->wireBackwardPosition(layerIDWithLayerOffset, wires[i]);
174
175 const HepPoint3D fwd(tfw3v.x(), tfw3v.y(), tfw3v.z());
176 const HepPoint3D bck(tbw3v.x(), tbw3v.y(), tbw3v.z());
177
178 if (m_magneticField && (abs(pid) != 99666)) {
179 // For monopoles a line segment approximation in the step volume is done,
180 // which is more reasonable, but should be done with a proper catenary FIXME
181 // Cal. distance assuming helix track (still approximation)
182 m_nonUniformField = 1;
183 if (Bfield[0] == 0. && Bfield[1] == 0. &&
184 Bfield[2] != 0.) m_nonUniformField = 0;
185
186 const G4double B_kG[3] = {Bfield[0] / CLHEP::kilogauss,
187 Bfield[1] / CLHEP::kilogauss,
188 Bfield[2] / CLHEP::kilogauss
189 };
190
191 const HepPoint3D x(pos[0] / CLHEP::cm, pos[1] / CLHEP::cm, pos[2] / CLHEP::cm);
192 const HepVector3D p(momIn.x() / CLHEP::GeV, momIn.y() / CLHEP::GeV, momIn.z() / CLHEP::GeV);
193 Helix tmp(x, p, charge);
194 tmp.bFieldZ(B_kG[2]);
195 tmp.ignoreErrorMatrix();
196
197 const HepVector3D wire = fwd - bck;
198 HepPoint3D tryp =
199 (x.z() - bck.z()) / wire.z() * wire + bck;
200 tmp.pivot(tryp);
201 tryp = (tmp.x(0.).z() - bck.z()) / wire.z() * wire + bck;
202 tmp.pivot(tryp);
203 tryp = (tmp.x(0.).z() - bck.z()) / wire.z() * wire + bck;
204 tmp.pivot(tryp);
205
206 distance = std::abs(tmp.a()[0]);
207 posW.setX(tryp.x());
208 posW.setY(tryp.y());
209 posW.setZ(tryp.z());
210
211 onTrack = tmp.x(0.);
212 pOnTrack = tmp.momentum(0.);
213
214 for_Rotat(B_kG);
215 const G4double xwb(bck.x()), ywb(bck.y()), zwb(bck.z());
216 const G4double xwf(fwd.x()), ywf(fwd.y()), zwf(fwd.z());
217 const G4double xp(onTrack.x()), yp(onTrack.y()), zp(onTrack.z());
218 const G4double px(pOnTrack.x()), py(pOnTrack.y()), pz(pOnTrack.z());
219 G4double q2[3] = {0.}, q1[3] = {0.}, q3[3] = {0.};
220 const G4int ntryMax(50); //tentative; too large probably...
221 G4double dist;
222 G4int ntry(999);
223 HELWIR(xwb, ywb, zwb, xwf, ywf, zwf,
224 xp, yp, zp, px, py, pz,
225 B_kG, charge, ntryMax, dist, q2, q1, q3, ntry);
226
227 if (ntry <= ntryMax) {
228 if (m_wireSag) {
229 G4double ywb_sag, ywf_sag;
230 m_cdcgp->getWireSagEffect(CDCGeometryPar::c_Base, layerIDWithLayerOffset, wires[i], q2[2], ywb_sag, ywf_sag);
231 HELWIR(xwb, ywb_sag, zwb, xwf, ywf_sag, zwf,
232 xp, yp, zp, px, py, pz,
233 B_kG, charge, ntryMax, dist, q2, q1, q3, ntry);
234 }
235 if (ntry <= ntryMax) {
236 distance = dist;
237 onTrack.setX(q1[0]);
238 onTrack.setY(q1[1]);
239 onTrack.setZ(q1[2]);
240 posW.setX(q2[0]);
241 posW.setY(q2[1]);
242 posW.setZ(q2[2]);
243 pOnTrack.setX(q3[0]);
244 pOnTrack.setY(q3[1]);
245 pOnTrack.setZ(q3[2]);
246 }
247 }
248 } else { //no magnetic field case
249 // Cal. distance assuming a line track
250 G4ThreeVector bwp(bck.x(), bck.y(), bck.z());
251 G4ThreeVector fwp(fwd.x(), fwd.y(), fwd.z());
252 G4ThreeVector hitPosition, wirePosition;
253 distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
254 hitPosition, wirePosition);
255 if (m_wireSag) {
256 G4double ywb_sag, ywf_sag;
257 m_cdcgp->getWireSagEffect(CDCGeometryPar::c_Base, layerIDWithLayerOffset, wires[i], wirePosition.z(), ywb_sag, ywf_sag);
258 bwp.setY(ywb_sag);
259 fwp.setY(ywf_sag);
260 distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
261 hitPosition, wirePosition);
262 }
263
264 onTrack.setX(hitPosition.x());
265 onTrack.setY(hitPosition.y());
266 onTrack.setZ(hitPosition.z());
267 posW.setX(wirePosition.x());
268 posW.setY(wirePosition.y());
269 posW.setZ(wirePosition.z());
270 //tentative setting
271 pOnTrack.setX(0.5 * (momIn.x() + momOut.x()) / CLHEP::GeV);
272 pOnTrack.setY(0.5 * (momIn.y() + momOut.y()) / CLHEP::GeV);
273 pOnTrack.setZ(0.5 * (momIn.z() + momOut.z()) / CLHEP::GeV);
274 } //end of magneticfiled on or off
275
276 distance *= CLHEP::cm; onTrack *= CLHEP::cm; posW *= CLHEP::cm;
277 pOnTrack *= CLHEP::GeV;
278
279 G4ThreeVector posTrack(onTrack.x(), onTrack.y(), onTrack.z());
280 G4ThreeVector mom(pOnTrack.x(), pOnTrack.y(), pOnTrack.z());
281
282 const B2Vector3D tPosW(posW.x(), posW.y(), posW.z());
283 const B2Vector3D tPosTrack(posTrack.x(), posTrack.y(), posTrack.z());
284 const B2Vector3D tMom(mom.x(), mom.y(), mom.z());
285 G4int lr = m_cdcgp->getOldLeftRight(tPosW, tPosTrack, tMom);
286 G4int newLrRaw = m_cdcgp->getNewLeftRightRaw(tPosW, tPosTrack, tMom);
287 G4int newLr = newLrRaw; //to be modified in EndOfEvent
288
289 if (nWires == 1) {
290 saveSimHit(layerIDWithLayerOffset, wires[i], trackID, pid, distance, tofBefore, edep, s_in_layer * CLHEP::cm, pOnTrack, posW, posIn,
291 posOut,
292 posTrack, lr, newLrRaw, newLr, speed, hitWeight);
293
294 } else {
295
296 G4int cel1 = wires[i] + 1;
297 G4int cel2 = cel1;
298 if (i + 1 <= nWires - 1) {
299 cel2 = wires[i + 1] + 1;
300 }
301 const G4double s2 = t.GetTrackLength() / CLHEP::cm; //at post-step
302 G4double s1 = (s2 - s_in_layer); //at pre-step; varied later
303 G4ThreeVector din = momIn;
304 if (din.mag() != 0.) din /= momIn.mag();
305
306 G4double vent[6] = {posIn.x() / CLHEP::cm, posIn.y() / CLHEP::cm, posIn.z() / CLHEP::cm, din.x(), din.y(), din.z()};
307
308 G4ThreeVector dot(momOut.x(), momOut.y(), momOut.z());
309 if (dot.mag() != 0.) {
310 dot /= dot.mag();
311 } else {
312 // Flight-direction is needed to set even when a particle stops
313 dot = din;
314 }
315
316 G4double vext[6] = {posOut.x() / CLHEP::cm, posOut.y() / CLHEP::cm, posOut.z() / CLHEP::cm, dot.x(), dot.y(), dot.z()};
317
318 if (i > 0) {
319 for (int j = 0; j < 6; ++j) vent[j] = xint[j];
320 s1 = sint;
321 }
322
323 G4int flag(0);
324 G4double edep_in_cell(0.);
325 G4double eLossInCell(0.);
326
327 if (cel1 != cel2) {
328 CellBound(layerIDWithLayerOffset, cel1, cel2, vent, vext, s1, s2, xint, sint, flag);
329 const G4double test = (sint - s1) / s_in_layer;
330 if (test < 0. || test > 1.) {
331 B2WARNING("CDCSensitiveDetector: Strange path length: " << "s1= " << s1 << " sint= " << sint << " s_in_layer= " << s_in_layer <<
332 " test= " << test);
333 }
334 edep_in_cell = edep * std::abs((sint - s1)) / s_in_layer;
335
336 const G4ThreeVector x_In(vent[0]*CLHEP::cm, vent[1]*CLHEP::cm, vent[2]*CLHEP::cm);
337 const G4ThreeVector x_Out(xint[0]*CLHEP::cm, xint[1]*CLHEP::cm, xint[2]*CLHEP::cm);
338 const G4ThreeVector p_In(momBefore * vent[3], momBefore * vent[4], momBefore * vent[5]);
339
340 saveSimHit(layerIDWithLayerOffset, wires[i], trackID, pid, distance, tofBefore, edep_in_cell, std::abs((sint - s1)) * CLHEP::cm,
341 pOnTrack, posW,
342 x_In, x_Out,
343 posTrack, lr, newLrRaw, newLr, speed, hitWeight);
344 tofBefore += (sint - s1) / speedInCmPerNs;
345 eLossInCell = eLoss * (sint - s1) / s_in_layer;
346 kinEnergyBefore -= eLossInCell;
347 if (kinEnergyBefore >= 0.) {
348 momBefore = sqrt(kinEnergyBefore * (kinEnergyBefore + 2.*mass));
349 } else {
350 B2WARNING("CDCSensitiveDetector: Kinetic Energy < 0.");
351 momBefore = 0.;
352 }
353
354 } else { //the particle exits
355
356 const G4double test = (s2 - sint) / s_in_layer;
357 if (test < 0. || test > 1.) {
358 B2WARNING("CDCSensitiveDetector: Strange path length: " << "s2= " << s2 << " sint= " << sint << " s_in_layer= " << s_in_layer <<
359 " test= " << test);
360 }
361 edep_in_cell = edep * std::abs((s2 - sint)) / s_in_layer;
362
363 const G4ThreeVector x_In(vent[0]*CLHEP::cm, vent[1]*CLHEP::cm, vent[2]*CLHEP::cm);
364 const G4ThreeVector p_In(momBefore * vent[3], momBefore * vent[4], momBefore * vent[5]);
365
366 saveSimHit(layerIDWithLayerOffset, wires[i], trackID, pid, distance, tofBefore, edep_in_cell, std::abs((s2 - sint)) * CLHEP::cm,
367 pOnTrack, posW,
368 x_In,
369 posOut, posTrack, lr, newLrRaw, newLr, speed, hitWeight);
370 }
371 }
372
373 } //end of wire loop
374
375 return true;
376 }
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

◆ WireId_in_hit_order()

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

Sort wire id.

Definition at line 1065 of file CDCSensitiveDetector.cc.

1066 {
1067 std::vector<int> list;
1068 int i0 = int(id0);
1069 int i1 = int(id1);
1070 if (abs(i0 - i1) * 2 < int(nWires)) {
1071 if (id0 < id1) {
1072 for (int i = id0; i <= id1; ++i)
1073 list.push_back(i);
1074 } else {
1075 for (int i = id0; i >= id1; i--) {
1076 list.push_back(i);
1077 }
1078 }
1079 } else {
1080 if (id0 < id1) {
1081 for (int i = id0; i >= 0; i--)
1082 list.push_back(i);
1083 for (int i = nWires - 1; i >= id1; i--)
1084 list.push_back(i);
1085 } else {
1086 for (int i = id0; i < nWires; ++i)
1087 list.push_back(i);
1088 for (int i = 0; i <= id1; ++i)
1089 list.push_back(i);
1090 }
1091 }
1092
1093 return list;
1094 }