Belle II Software  release-08-01-10
CDCSensitiveDetector.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 #include <cdc/simulation/CDCSensitiveDetector.h>
10 
11 #include <cdc/simulation/CDCSimControlPar.h>
12 #include <cdc/simulation/Helix.h>
13 #include <cdc/geometry/CDCGeometryPar.h>
14 #include <cdc/utilities/ClosestApproach.h>
15 #include <framework/logging/Logger.h>
16 #include <framework/datastore/RelationArray.h>
17 #include <framework/geometry/B2Vector3.h>
18 
19 #include "G4Step.hh"
20 #include "G4TransportationManager.hh"
21 #include "G4Field.hh"
22 #include "G4FieldManager.hh"
23 
24 #include "CLHEP/Geometry/Vector3D.h"
25 #include "CLHEP/Geometry/Point3D.h"
26 
27 
28 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
29 typedef HepGeom::Point3D<double> HepPoint3D;
30 #endif
31 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
32 typedef HepGeom::Vector3D<double> HepVector3D;
33 #endif
34 
35 
36 
37 namespace Belle2 {
42  using namespace CDC;
43 
44  //N.B.#0: Do not put CDCGeometryPar::Instance() in the initializing list of the constructor,
45  //because it is called before CDCGeometryPar(geom) is called in case of UseDB=true of Geometry module.
46  //N.B.#1: Do not call AddNewDetector(), because it'll cause a job crash currently.
47 
48  CDCSensitiveDetector::CDCSensitiveDetector(G4String name, G4double thresholdEnergyDeposit, G4double thresholdKineticEnergy):
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  }
88 
89  void CDCSensitiveDetector::Initialize(G4HCofThisEvent*)
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  }
111 
112  //-----------------------------------------------------
113  // Method invoked for every step in sensitive detector
114  //-----------------------------------------------------
115  bool CDCSensitiveDetector::step(G4Step* aStep, G4TouchableHistory*)
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  }
572 
573  /*
574  void CDCSensitiveDetector::BeginOfEvent(G4HCofThisEvent*)
575  {
576  std::cout <<"CDCSensitiveDetector::BeginOfEvent callded." << std::endl;
577  }
578  */
579 
580  void CDCSensitiveDetector::EndOfEvent(G4HCofThisEvent*)
581  {
582  setModifiedLeftRightFlag();
583  }
584 
585  void
586  CDCSensitiveDetector::saveSimHit(const G4int layerId,
587  const G4int wireId,
588  const G4int trackID,
589  const G4int pid,
590  const G4double distance,
591  const G4double tof,
592  const G4double edep,
593  const G4double stepLength,
594  const G4ThreeVector& mom,
595  const G4ThreeVector& posW,
596  const G4ThreeVector& posIn,
597  const G4ThreeVector& posOut,
598  const G4ThreeVector& posTrack,
599  const G4int lr,
600  const G4int newLrRaw,
601  const G4int newLr,
602  const G4double speed,
603  const G4double hitWeight)
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  }
671 
672 
673  /*
674  void SensitiveDetector::AddbgOne(bool doit) {
675  Belle::Datcdc_olhit_Manager& olhitmgr=Belle::Datcdc_olhit_Manager::get_manager();
676  if(doit) {
677  for( int i=0;i<olhitmgr.count();i++ ){
678  Belle::Datcdc_olhit& h = olhitmgr[i];
679  }
680  dout(Debugout::B2INFO,"SensitiveDetector")
681  << "AddbgOne " << olhitmgr.size()
682  << std::endl;
683  }
684  olhitmgr.remove();
685  }
686  */
687 
688  void
689  CDCSensitiveDetector::CellBound(const G4int layerId,
690  const G4int ic1,
691  const G4int ic2,
692  const G4double venter[6],
693  const G4double vexit[6],
694  const G4double s1, const G4double s2,
695  G4double xint[6],
696  G4double& sint, G4int& iflag)
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  }
935 
936  void CDCSensitiveDetector::GCUBS(const G4double x, const G4double y, const G4double d1, const G4double d2, G4double a[4])
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 
963 L10:
964  a[0] = 0.;
965  a[1] = 1.;
966  a[2] = 0.;
967  a[3] = 0.;
968  }
969 
970  void
971  CDCSensitiveDetector::for_Rotat(const G4double bfld[3])
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 tentaive 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  }
1005 
1006  void
1007  CDCSensitiveDetector::Rotat(G4double& x, G4double& y, G4double& z,
1008  const int mode)
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 tentaive 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 << "specifed");
1029  }
1030  return;
1031 
1032  }
1033 
1034  void
1035  CDCSensitiveDetector::Rotat(G4double x[3], const int mode)
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 tentaive 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 << "specifed");
1056  }
1057  return;
1058 
1059  }
1060 
1061  void
1062  CDCSensitiveDetector::HELWIR(const G4double xwb4, const G4double ywb4,
1063  const G4double zwb4,
1064  const G4double xwf4, const G4double ywf4,
1065  const G4double zwf4,
1066  const G4double xp, const G4double yp,
1067  const G4double zp,
1068  const G4double px, const G4double py,
1069  const G4double pz,
1070  const G4double B_kG[3],
1071  const G4double charge, const G4int ntryMax,
1072  G4double& distance,
1073  G4double q2[3], G4double q1[3],
1074  G4double q3[3],
1075  G4int& ntry)
1076  {
1077  //~dead copy of gsim_cdc_hit.F in gsim-cdc for Belle (for tentaive 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 distingushing 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;
1201 line1:
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 pont 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.
1292 line100:
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  }
1330 
1331  void
1332  CDCSensitiveDetector::Mvopr(const G4int ndim, const G4double b[3], const G4double m[3][3],
1333  const G4double a[3], G4double c[3], const G4int mode)
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  }
1379 
1380  std::vector<int>
1381  CDCSensitiveDetector::WireId_in_hit_order(int id0, int id1, int nWires)
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  }
1411 
1412  G4double CDCSensitiveDetector::ClosestApproach(const G4ThreeVector bwp, const G4ThreeVector fwp, const G4ThreeVector posIn,
1413  const G4ThreeVector posOut, G4ThreeVector& hitPosition, G4ThreeVector& wirePosition)//,G4double& transferT)
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  }
1436 
1437 
1438  //The following-to-end is for setting of left/right flag modified for tracking
1439  void CDCSensitiveDetector::setModifiedLeftRightFlag()
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  }
1521 
1522 
1523  void CDCSensitiveDetector::reAssignLeftRightInfo()
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  }
1607 
1608 
1609  unsigned short CDCSensitiveDetector::areNeighbors(const WireID& wireId, const WireID& otherWireId) const
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  }
1670 
1671  unsigned short CDCSensitiveDetector::areNeighbors(unsigned short iCLayer, unsigned short iSuperLayer, unsigned short iLayer,
1672  unsigned short iWire, const WireID& otherWireId) const
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  }
1731 
1733 } // namespace Belle2
Belle2::B2Vector3::y
DataType y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:427
Belle2::B2Vector3::z
DataType z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:429
Belle2::B2Vector3::x
DataType x() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:425
Belle2::CDCSimHit
Example Detector.
Definition: CDCSimHit.h:21
Belle2::CDCSimHit::setLeftRightPassageRaw
void setLeftRightPassageRaw(int minusOneOrZeroOrOne)
The method to set new left/right info. for digitization.
Definition: CDCSimHit.h:156
Belle2::CDCSimHit::setPosIn
void setPosIn(const B2Vector3D &posIn)
The method to set position of pre-step.
Definition: CDCSimHit.h:122
Belle2::CDCSimHit::setPosFlag
void setPosFlag(int zeroOrOne)
The method to set position flag.
Definition: CDCSimHit.h:149
Belle2::CDCSimHit::getPosWire
B2Vector3D getPosWire() const
The method to get position on wire.
Definition: CDCSimHit.h:199
Belle2::CDCSimHit::getWireID
WireID getWireID() const
Getter for WireID object.
Definition: CDCSimHit.h:172
Belle2::CDCSimHit::setWireID
void setWireID(int iCLayerID, int iWireID)
Setter for Wire ID.
Definition: CDCSimHit.h:79
Belle2::CDCSimHit::setGlobalTime
void setGlobalTime(double globalTime)
The method to set global time.
Definition: CDCSimHit.h:97
Belle2::CDCSimHit::setTrackId
void setTrackId(int trackId)
The method to set track id.
Definition: CDCSimHit.h:85
Belle2::CDCSimHit::setPDGCode
void setPDGCode(int pdg)
The method to set PDG code.
Definition: CDCSimHit.h:88
Belle2::CDCSimHit::setMomentum
void setMomentum(const B2Vector3D &momentum)
The method to set momentum.
Definition: CDCSimHit.h:106
Belle2::CDCSimHit::setPosWire
void setPosWire(const B2Vector3D &posWire)
The method to set position on wire.
Definition: CDCSimHit.h:114
Belle2::CDCSimHit::setEnergyDep
void setEnergyDep(double edep)
The method to set deposited energy.
Definition: CDCSimHit.h:100
Belle2::CDCSimHit::setPosOut
void setPosOut(const B2Vector3D &posOut)
The method to set position of post-step.
Definition: CDCSimHit.h:130
Belle2::CDCSimHit::setDriftLength
void setDriftLength(double driftLength)
The method to set drift length.
Definition: CDCSimHit.h:91
Belle2::CDCSimHit::setStepLength
void setStepLength(double stepLength)
The method to set step length.
Definition: CDCSimHit.h:103
Belle2::CDCSimHit::getPosTrack
B2Vector3D getPosTrack() const
The method to get position on the track.
Definition: CDCSimHit.h:217
Belle2::CDCSimHit::getMomentum
B2Vector3D getMomentum() const
The method to get momentum.
Definition: CDCSimHit.h:193
Belle2::CDCSimHit::setFlightTime
void setFlightTime(double flightTime)
The method to set flight time.
Definition: CDCSimHit.h:94
Belle2::CDCSimHit::setPosTrack
void setPosTrack(const B2Vector3D &posTrack)
The method to set position on the track.
Definition: CDCSimHit.h:138
Belle2::CDCSimHit::setLeftRightPassage
void setLeftRightPassage(int minusOneOrZeroOrOne)
The method to set new left/right info. for tracking.
Definition: CDCSimHit.h:164
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:1838
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:1967
Belle2::CDC::CDCGeometryPar::getOffsetOfFirstLayer
ushort getOffsetOfFirstLayer() const
Get the offset of the first layer.
Definition: CDCGeometryPar.h:1071
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:1706
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:1755
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1250
Belle2::CDC::CDCGeometryPar::getShiftInSuperLayer
signed short getShiftInSuperLayer(unsigned short iSuperLayer, unsigned short iLayer) const
Returns shift in the super-layer.
Definition: CDCGeometryPar.cc:3031
Belle2::CDC::CDCGeometryPar::getOldLeftRight
unsigned short getOldLeftRight(const B2Vector3D &posOnWire, const B2Vector3D &posOnTrack, const B2Vector3D &momentum) const
Returns old left/right.
Definition: CDCGeometryPar.cc:2758
Belle2::CDC::CDCGeometryPar::Instance
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Definition: CDCGeometryPar.cc:39
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:2793
Belle2::CDC::CDCSensitiveDetector::CW_NEIGHBOR
const signed short CW_NEIGHBOR
Constant for clockwise.
Definition: CDCSensitiveDetector.h:293
Belle2::CDC::CDCSensitiveDetector::m_thresholdEnergyDeposit
G4double m_thresholdEnergyDeposit
Threshold energy deposit to be stored.
Definition: CDCSensitiveDetector.h:264
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::m_hitWithPosWeight
std::multimap< unsigned short, CDCSimHit * > m_hitWithPosWeight
Map containing hits with positive weight.
Definition: CDCSensitiveDetector.h:280
Belle2::CDC::CDCSensitiveDetector::CCW_OUT_NEIGHBOR
const signed short CCW_OUT_NEIGHBOR
Constant for counterclockwise outwards.
Definition: CDCSensitiveDetector.h:297
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::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::m_modifiedLeftRightFlag
G4bool m_modifiedLeftRightFlag
Switch for left/right flag modified for tracking.
Definition: CDCSensitiveDetector.h:273
Belle2::CDC::CDCSensitiveDetector::m_magneticField
G4int m_magneticField
Magnetic field is on or off.
Definition: CDCSensitiveDetector.h:248
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_brot
G4double m_brot[3][3]
a rotation matrix.
Definition: CDCSensitiveDetector.h:256
Belle2::CDC::CDCSensitiveDetector::m_hitWithNegWeight
std::vector< CDCSimHit * > m_hitWithNegWeight
Vector containing hits with negative weight.
Definition: CDCSensitiveDetector.h:282
Belle2::CDC::CDCSensitiveDetector::m_nonUniformField
G4int m_nonUniformField
Magnetic field is uniform or non-uniform.
Definition: CDCSensitiveDetector.h:254
Belle2::CDC::CDCSensitiveDetector::m_MCParticles
StoreArray< MCParticle > m_MCParticles
MC particles.
Definition: CDCSensitiveDetector.h:240
Belle2::CDC::CDCSensitiveDetector::CW
const signed short CW
Constant for clockwise orientation.
Definition: CDCSensitiveDetector.h:291
Belle2::CDC::CDCSimControlPar
The Class for CDC Simulation Control Parameters.
Definition: CDCSimControlPar.h:20
Belle2::CDC::CDCSimControlPar::getModLeftRightFlag
bool getModLeftRightFlag() const
Get modified left/right flag.
Definition: CDCSimControlPar.h:100
Belle2::CDC::CDCSimControlPar::getMinTrackLength
double getMinTrackLength() const
Get minimum track length.
Definition: CDCSimControlPar.h:124
Belle2::CDC::CDCSimControlPar::getThresholdEnergyDeposit
double getThresholdEnergyDeposit() const
Get threshold for Energy Deposit;.
Definition: CDCSimControlPar.h:116
Belle2::CDC::CDCSimControlPar::getInstance
static CDCSimControlPar & getInstance()
Static method to get a reference to the CDCSimControlPar instance.
Definition: CDCSimControlPar.cc:16
Belle2::CDC::CDCSimControlPar::getWireSag
bool getWireSag() const
Get wiresag flag.
Definition: CDCSimControlPar.h:92
Belle2::CDC::Helix
Helix parameter class.
Definition: Helix.h:48
Belle2::Const
This class provides a set of constants for the framework.
Definition: Const.h:34
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
Belle2::RelationArray::getEntries
int getEntries() const
Get the number of elements.
Definition: RelationArray.h:242
Belle2::RelationArray::add
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Definition: RelationArray.h:281
Belle2::RelationElement
Class to store a single element of a relation.
Definition: RelationElement.h:25
Belle2::RelationElement::getFromIndex
index_type getFromIndex() const
Get index we point from.
Definition: RelationElement.h:70
Belle2::RelationElement::getToIndex
index_type getToIndex(size_t n=0) const
Get nth index we point to.
Definition: RelationElement.h:79
Belle2::RelationElement::setToIndex
void setToIndex(index_type to, weight_type weight=1.0)
Set index we point to, converts relation to 1:1 and discards all existing to-indices.
Definition: RelationElement.h:94
Belle2::RelationElement::weight_type
float weight_type
type used for weights.
Definition: RelationElement.h:32
Belle2::RelationElement::getSize
size_t getSize() const
Get number of indices we points to.
Definition: RelationElement.h:73
Belle2::RelationElement::getWeight
weight_type getWeight(size_t n=0) const
Get nth weight we point to.
Definition: RelationElement.h:82
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::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::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
Belle2::WireID
Class to identify a wire inside the CDC.
Definition: WireID.h:34
Belle2::WireID::getICLayer
unsigned short getICLayer() const
Getter for continuous layer numbering.
Definition: WireID.cc:24
Belle2::WireID::getIWire
unsigned short getIWire() const
Getter for wire within the layer.
Definition: WireID.h:145
Belle2::WireID::getISuperLayer
unsigned short getISuperLayer() const
Getter for Super-Layer.
Definition: WireID.h:130
Belle2::WireID::getILayer
unsigned short getILayer() const
Getter for layer within the Super-Layer.
Definition: WireID.h:136
Belle2::CDC::CDCSensitiveDetector::step
bool step(G4Step *aStep, G4TouchableHistory *history) override
Process each step and calculate variables defined in CDCB4VHit.
Definition: CDCSensitiveDetector.cc:115
Belle2::CDC::CDCSensitiveDetector::setModifiedLeftRightFlag
void setModifiedLeftRightFlag()
set left/right flag modified for tracking
Definition: CDCSensitiveDetector.cc:1439
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
Belle2::CDC::CDCSensitiveDetector::EndOfEvent
void EndOfEvent(G4HCofThisEvent *) override
Do what you want to do at the beginning of each event (why this is not called ?)
Definition: CDCSensitiveDetector.cc:580
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::CDCSensitiveDetector
CDCSensitiveDetector(G4String name, G4double thresholdEnergyDeposit, G4double thresholdKineticEnergy)
Constructor.
Definition: CDCSensitiveDetector.cc:48
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::Rotat
void Rotat(G4double &x, G4double &y, G4double &z, const int mode)
Translation method.
Definition: CDCSensitiveDetector.cc:1007
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::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::CDC::CDCSensitiveDetector::Initialize
void Initialize(G4HCofThisEvent *) override
Register CDC hits collection into G4HCofThisEvent.
Definition: CDCSensitiveDetector.cc:89
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::GCUBS
void GCUBS(const G4double x, const G4double y, const G4double d1, const G4double d2, G4double a[4])
Definition: CDCSensitiveDetector.cc:936
Belle2::CDC::CDCSensitiveDetector::reAssignLeftRightInfo
void reAssignLeftRightInfo()
Re-assign left/right info.
Definition: CDCSensitiveDetector.cc:1523
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
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition: beamHelpers.h:28
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::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
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:17