Belle II Software  release-05-02-19
CDCSensitiveDetector.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Guofu Cao *
7  * 2012.03.05 SensitiveDetector -> CDCSensitiveDetector by M. Uchida *
8  * *
9  * This software is provided "as is" without any warranty. *
10  **************************************************************************/
11 
12 #include <cdc/simulation/CDCSensitiveDetector.h>
13 
14 #include <cdc/simulation/CDCSimControlPar.h>
15 #include <cdc/simulation/Helix.h>
16 #include <cdc/geometry/CDCGeometryPar.h>
17 #include <cdc/utilities/ClosestApproach.h>
18 #include <framework/logging/Logger.h>
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/datastore/RelationArray.h>
21 #include <cdc/dataobjects/CDCSimHit.h>
22 
23 #include "G4Step.hh"
24 #include "G4TransportationManager.hh"
25 #include "G4Field.hh"
26 #include "G4FieldManager.hh"
27 
28 #include "CLHEP/Geometry/Vector3D.h"
29 #include "CLHEP/Geometry/Point3D.h"
30 
31 #include "TVector3.h"
32 
33 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
34 typedef HepGeom::Point3D<double> HepPoint3D;
35 #endif
36 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
37 typedef HepGeom::Vector3D<double> HepVector3D;
38 #endif
39 
40 
41 
42 namespace Belle2 {
47  using namespace CDC;
48 
49  //N.B.#0: Do not put CDCGeometryPar::Instance() in the initializing list of the constructor,
50  //because it is called before CDCGeometryPar(geom) is called in case of UseDB=true of Geometry module.
51  //N.B.#1: Do not call AddNewDetector(), because it'll cause a job crash currently.
52 
53  CDCSensitiveDetector::CDCSensitiveDetector(G4String name, G4double thresholdEnergyDeposit, G4double thresholdKineticEnergy):
54  SensitiveDetectorBase(name, Const::CDC),
55  // m_cdcgp(CDCGeometryPar::Instance()),
56  m_cdcgp(nullptr),
57  m_thresholdEnergyDeposit(thresholdEnergyDeposit),
58  m_thresholdKineticEnergy(thresholdKineticEnergy), m_hitNumber(0)
59  {
60  StoreArray<MCParticle> mcParticles;
61  StoreArray<CDCSimHit> cdcSimHits;
62  RelationArray cdcSimHitRel(mcParticles, cdcSimHits);
63  registerMCParticleRelation(cdcSimHitRel);
64  // registerMCParticleRelation(cdcSimHitRel, RelationArray::c_doNothing);
65  // registerMCParticleRelation(cdcSimHitRel, RelationArray::c_negativeWeight);
66  // registerMCParticleRelation(cdcSimHitRel, RelationArray::c_deleteElement);
67  // registerMCParticleRelation(cdcSimHitRel, RelationArray::c_zeroWeight);
68  cdcSimHits.registerInDataStore();
69  mcParticles.registerRelationTo(cdcSimHits);
70 
71  CDCSimControlPar& cntlp = CDCSimControlPar::getInstance();
72 
73  m_thresholdEnergyDeposit = cntlp.getThresholdEnergyDeposit();
74  m_thresholdEnergyDeposit *= CLHEP::GeV; //GeV to MeV (=unit in G4)
75  B2DEBUG(150, "CDCSensitiveDetector: Threshold energy (MeV): " << m_thresholdEnergyDeposit);
76  m_thresholdKineticEnergy = 0.0; // Dummy to avoid a warning (tentative).
77 
78  m_wireSag = cntlp.getWireSag();
79  // m_wireSag = false;
80  B2DEBUG(150, "CDCSensitiveDetector: Sense wire sag on(=1)/off(=0): " << m_wireSag);
81 
82  m_modifiedLeftRightFlag = cntlp.getModLeftRightFlag();
83  B2DEBUG(150, "CDCSensitiveDetector: Set left/right flag modified for tracking (=1)/ not set (=0): " << m_modifiedLeftRightFlag);
84 
85  m_minTrackLength = cntlp.getMinTrackLength();
86  m_minTrackLength *= CLHEP::cm; //cm to mm (=unit in G4)
87  B2DEBUG(150, "CDCSensitiveDetector: MinTrackLength (mm): " << m_minTrackLength);
88 
89  //For activating Initialize and EndOfEvent functions
90  //but not work --> commented out for a while...
91  // if (m_modifiedLeftRightFlag) {
92  // G4SDManager::GetSDMpointer()->AddNewDetector(this);
93  // }
94  }
95 
96  void CDCSensitiveDetector::Initialize(G4HCofThisEvent*)
97  {
98  /*
99  m_cdcgp = &CDCGeometryPar::Instance();
100 
101  m_thresholdEnergyDeposit = m_cdcgp->getThresholdEnergyDeposit();
102  m_thresholdEnergyDeposit *= CLHEP::GeV; //GeV to MeV (=unit in G4)
103  B2INFO("CDCSensitiveDetector: Threshold energy (MeV): " << m_thresholdEnergyDeposit);
104  m_modifiedLeftRightFlag = m_cdcgp->isModifiedLeftRightFlagOn();
105  B2INFO("CDCSensitiveDetector: Set left/right flag modified for tracking (=1)/ not set (=0): " << m_modifiedLeftRightFlag);
106 
107  m_minTrackLength = m_cdcgp->getMinTrackLength();
108  m_minTrackLength *= CLHEP::cm; //cm to mm (=unit in G4)
109  B2INFO("CDCSensitiveDetector: MinTrackLength (mm): " << m_minTrackLength);
110  */
111 
112  // Initialize
113  m_nonUniformField = 0;
114  // std::cout << "Initialize called" << std::endl;
115  // m_nPosHits = 0;
116  // m_nNegHits = 0;
117  }
118 
119  //-----------------------------------------------------
120  // Method invoked for every step in sensitive detector
121  //-----------------------------------------------------
122  bool CDCSensitiveDetector::step(G4Step* aStep, G4TouchableHistory*)
123  {
124  // static bool firstCall = true;
125  // if (firstCall) {
126  // firstCall = false;
127  m_cdcgp = &CDCGeometryPar::Instance();
128  // CDCSimControlPar & m_cntlp = CDCSimControlPar::getInstance();
129 
130  // // m_thresholdEnergyDeposit = m_cdcgp->getThresholdEnergyDeposit();
131  // m_thresholdEnergyDeposit = m_cntlp.getThresholdEnergyDeposit();
132  // m_thresholdEnergyDeposit *= CLHEP::GeV; //GeV to MeV (=unit in G4)
133  // B2INFO("CDCSensitiveDetector: Threshold energy (MeV): " << m_thresholdEnergyDeposit);
134 
135  // // m_modifiedLeftRightFlag = m_cdcgp->isModifiedLeftRightFlagOn();
136  // m_modifiedLeftRightFlag = m_cntlp.getModLeftRightFlag();
137  // B2INFO("CDCSensitiveDetector: Set left/right flag modified for tracking (=1)/ not set (=0): " << m_modifiedLeftRightFlag);
138 
139  // // m_minTrackLength = m_cdcgp->getMinTrackLength();
140  // m_minTrackLength = m_cntlp.getMinTrackLength();
141  // m_minTrackLength *= CLHEP::cm; //cm to mm (=unit in G4)
142  // B2INFO("CDCSensitiveDetector: MinTrackLength (mm): " << m_minTrackLength);
143 
144  // // m_wireSag = m_cdcgp->isWireSagOn();
145  // m_wireSag = m_cntlp.getWireSag();
146 
147  m_nonUniformField = 0;
148  // }
149 
150 #if defined(CDC_DEBUG)
151  std::cout << " " << std::endl;
152  std::cout << "********* step in ********" << std::endl;
153 #endif
154  // Get deposited energy
155  const G4double edep = aStep->GetTotalEnergyDeposit();
156 
157 
158  // Discard the hit below Edep_th
159  if (edep <= m_thresholdEnergyDeposit) return false;
160 
161  // Get step length
162  const G4double stepLength = aStep->GetStepLength();
163  if (stepLength == 0.) return false;
164 
165  // Get step information
166  const G4Track& t = * aStep->GetTrack();
167 
168  G4double hitWeight = Simulation::TrackInfo::getInfo(t).getIgnore() ? -1 : 1;
169  // save in MCParticle if track-length is enough long
170  if (t.GetTrackLength() > m_minTrackLength) {
171  Simulation::TrackInfo::getInfo(t).setIgnore(false);
172  hitWeight = 1.;
173  // std::cout <<"setignore=false for track= "<< t.GetTrackID() << std::endl;
174  // } else {
175  // std::cout <<"setignore=true for track= "<< t.GetTrackID() << std::endl;
176  }
177 
178  // const G4double tof = t.GetGlobalTime(); //tof at post step point
179  // if (isnan(tof)) {
180  // B2ERROR("SensitiveDetector: global time is nan");
181  // return false;
182  // }
183 
184  const G4int pid = t.GetDefinition()->GetPDGEncoding();
185  const G4double charge = t.GetDefinition()->GetPDGCharge();
186  const G4int trackID = t.GetTrackID();
187  // std::cout << "pid,stepl,trackID,trackl,weight= " << pid <<" "<< stepLength <<" "<< trackID <<" "<< t.GetTrackLength() <<" "<< hitWeight << std::endl;
188 
189  const G4VPhysicalVolume& v = * t.GetVolume();
190  const G4StepPoint& in = * aStep->GetPreStepPoint();
191  const G4StepPoint& out = * aStep->GetPostStepPoint();
192  const G4ThreeVector& posIn = in.GetPosition();
193  const G4ThreeVector& posOut = out.GetPosition();
194  const G4ThreeVector momIn(in.GetMomentum().x(), in.GetMomentum().y(),
195  in.GetMomentum().z());
196 #if defined(CDC_DEBUG)
197  std::cout << "pid = " << pid << std::endl;
198  std::cout << "mass = " << t.GetDefinition()->GetPDGMass() << std::endl;
199  std::cout << "posIn = " << posIn << std::endl;
200  std::cout << "posOut= " << posOut << std::endl;
201  std::cout << "tof at post-step = " << out.GetGlobalTime() << std::endl;
202  std::cout << "stepl = " << stepLength << std::endl;
203 #endif
204 
205  // Get layer ID
206  const unsigned layerId = v.GetCopyNo();
207  B2DEBUG(150, "LayerID in continuous counting method: " << layerId);
208 
209  // If neutral particles, ignore them, unless monopoles.
210 
211  if ((charge == 0.) && (abs(pid) != 99666)) return false;
212 
213  // Calculate cell ID
214  TVector3 tposIn(posIn.x() / CLHEP::cm, posIn.y() / CLHEP::cm, posIn.z() / CLHEP::cm);
215  TVector3 tposOut(posOut.x() / CLHEP::cm, posOut.y() / CLHEP::cm, posOut.z() / CLHEP::cm);
216  const unsigned idIn = m_cdcgp->cellId(layerId, tposIn);
217  const unsigned idOut = m_cdcgp->cellId(layerId, tposOut);
218 #if defined(CDC_DEBUG)
219  std::cout << "edep= " << edep << std::endl;
220  std::cout << "idIn,idOut= " << idIn << " " << idOut << std::endl;
221 #endif
222 
223  // Calculate drift length
224  std::vector<int> wires = WireId_in_hit_order(idIn, idOut, m_cdcgp->nWiresInLayer(layerId));
225  G4double sint(0.);
226  const G4double s_in_layer = stepLength / CLHEP::cm;
227  G4double xint[6] = {0};
228 
229  const G4ThreeVector momOut(out.GetMomentum().x(), out.GetMomentum().y(),
230  out.GetMomentum().z());
231  const G4double speedIn = in.GetVelocity();
232  const G4double speedOut = out.GetVelocity();
233  const G4double speed = 0.5 * (speedIn + speedOut);
234  const G4double speedInCmPerNs = speed / CLHEP::cm;
235 
236  const unsigned int nWires = wires.size();
237  G4double tofBefore = in.GetGlobalTime();
238  G4double kinEnergyBefore = in.GetKineticEnergy();
239  G4double momBefore = momIn.mag();
240  const G4double eLoss = kinEnergyBefore - out.GetKineticEnergy(); //n.b. not always equal to edep
241  const G4double mass = t.GetDefinition()->GetPDGMass();
242 #if defined(CDC_DEBUG)
243  std::cout << "momBefore = " << momBefore << std::endl;
244  std::cout << "momIn = " << momIn.x() << " " << momIn.y() << " " << momIn.z() << std::endl;
245  std::cout << "momOut= " << momOut.x() << " " << momOut.y() << " " << momOut.z() << std::endl;
246  std::cout << "speedIn,speedOut= " << speedIn << " " << speedOut << std::endl;
247  std::cout << " speedInCmPerNs= " << speedInCmPerNs << std::endl;
248  std::cout << "tofBefore= " << tofBefore << std::endl;
249 #endif
250 
251  const G4Field* field = G4TransportationManager::GetTransportationManager()->GetFieldManager()->GetDetectorField();
252 
253  for (unsigned i = 0; i < nWires; ++i) {
254 #if defined(CDC_DEBUG)
255  std::cout << "============ i,wires[i]= " << i << " " << wires[i] << std::endl;
256 #endif
257 
258  const G4double pos[3] = {posIn.x(), posIn.y(), posIn.z()};
259  G4double Bfield[3];
260  field->GetFieldValue(pos, Bfield);
261  m_magneticField = (Bfield[0] == 0. && Bfield[1] == 0. &&
262  Bfield[2] == 0.) ? false : true;
263 #if defined(CDC_DEBUG)
264  std::cout << "Bfield= " << Bfield[0] << " " << Bfield[1] << " " << Bfield[2] << std::endl;
265  std::cout << "magneticField= " << m_magneticField << std::endl;
266 #endif
267 
268  double distance = 0;
269  G4ThreeVector posW(0, 0, 0);
270  HepPoint3D onTrack;
271  HepPoint3D pOnTrack;
272 
273  // Calculate forward/backward position of current wire
274  const TVector3 tfw3v = m_cdcgp->wireForwardPosition(layerId, wires[i]);
275  const TVector3 tbw3v = m_cdcgp->wireBackwardPosition(layerId, wires[i]);
276 
277  const HepPoint3D fwd(tfw3v.x(), tfw3v.y(), tfw3v.z());
278  const HepPoint3D bck(tbw3v.x(), tbw3v.y(), tbw3v.z());
279 
280  if (m_magneticField && (abs(pid) != 99666)) {
281  // For monopoles a line segment approximation in the step volume is done,
282  // which is more reasonable, but should be done with a proper catenary FIXME
283  // Cal. distance assuming helix track (still approximation)
284  m_nonUniformField = 1;
285  if (Bfield[0] == 0. && Bfield[1] == 0. &&
286  Bfield[2] != 0.) m_nonUniformField = 0;
287 
288  const G4double B_kG[3] = {Bfield[0] / CLHEP::kilogauss,
289  Bfield[1] / CLHEP::kilogauss,
290  Bfield[2] / CLHEP::kilogauss
291  };
292 #if defined(CDC_DEBUG)
293  std::cout << "B_kG= " << B_kG[0] << " " << B_kG[1] << " " << B_kG[2] << std::endl;
294  std::cout << "magneticField= " << m_magneticField << std::endl;
295 #endif
296 
297  const HepPoint3D x(pos[0] / CLHEP::cm, pos[1] / CLHEP::cm, pos[2] / CLHEP::cm);
298  const HepVector3D p(momIn.x() / CLHEP::GeV, momIn.y() / CLHEP::GeV, momIn.z() / CLHEP::GeV);
299  Helix tmp(x, p, charge);
300  tmp.bFieldZ(B_kG[2]);
301  tmp.ignoreErrorMatrix();
302 
303  /* // Calculate forward/backward position of current wire
304  const TVector3 tfw3v = cdcg.wireForwardPosition(layerId, wires[i]);
305  const TVector3 tbw3v = cdcg.wireBackwardPosition(layerId, wires[i]);
306 
307  const HepPoint3D fwd(tfw3v.x(), tfw3v.y(), tfw3v.z());
308  const HepPoint3D bck(tbw3v.x(), tbw3v.y(), tbw3v.z());
309  */
310 
311  const HepVector3D wire = fwd - bck;
312  HepPoint3D tryp =
313  (x.z() - bck.z()) / wire.z() * wire + bck;
314  tmp.pivot(tryp);
315  tryp = (tmp.x(0.).z() - bck.z()) / wire.z() * wire + bck;
316  tmp.pivot(tryp);
317  tryp = (tmp.x(0.).z() - bck.z()) / wire.z() * wire + bck;
318  tmp.pivot(tryp);
319 
320  distance = std::abs(tmp.a()[0]);
321  posW.setX(tryp.x());
322  posW.setY(tryp.y());
323  posW.setZ(tryp.z());
324 
325  // HepPoint3D onTrack = tmp.x(0.);
326  onTrack = tmp.x(0.);
327  pOnTrack = tmp.momentum(0.);
328 
329  for_Rotat(B_kG);
330  const G4double xwb(bck.x()), ywb(bck.y()), zwb(bck.z());
331  const G4double xwf(fwd.x()), ywf(fwd.y()), zwf(fwd.z());
332  const G4double xp(onTrack.x()), yp(onTrack.y()), zp(onTrack.z());
333  const G4double px(pOnTrack.x()), py(pOnTrack.y()), pz(pOnTrack.z());
334  G4double q2[3] = {0.}, q1[3] = {0.}, q3[3] = {0.};
335  const G4int ntryMax(50); //tentative; too large probably...
336  G4double dist;
337  G4int ntry(999);
338  HELWIR(xwb, ywb, zwb, xwf, ywf, zwf,
339  xp, yp, zp, px, py, pz,
340  B_kG, charge, ntryMax, dist, q2, q1, q3, ntry);
341 
342 #if defined(CDC_DEBUG)
343  std::cout << "ntry= " << ntry << std::endl;
344  std::cout << "bf distance= " << distance << std::endl;
345  std::cout << "onTrack = " << onTrack << std::endl;
346  std::cout << "posW = " << posW << std::endl;
347 #endif
348  if (ntry <= ntryMax) {
349  if (m_wireSag) {
350  G4double ywb_sag, ywf_sag;
351  m_cdcgp->getWireSagEffect(CDCGeometryPar::c_Base, layerId, wires[i], q2[2], ywb_sag, ywf_sag);
352  HELWIR(xwb, ywb_sag, zwb, xwf, ywf_sag, zwf,
353  xp, yp, zp, px, py, pz,
354  B_kG, charge, ntryMax, dist, q2, q1, q3, ntry);
355  }
356  if (ntry <= ntryMax) {
357  distance = dist;
358  onTrack.setX(q1[0]);
359  onTrack.setY(q1[1]);
360  onTrack.setZ(q1[2]);
361  posW.setX(q2[0]);
362  posW.setY(q2[1]);
363  posW.setZ(q2[2]);
364  pOnTrack.setX(q3[0]);
365  pOnTrack.setY(q3[1]);
366  pOnTrack.setZ(q3[2]);
367  }
368 #if defined(CDC_DEBUG)
369  std::cout << " " << std::endl;
370  std::cout << "helix distance= " << distance << std::endl;
371  std::cout << "onTrack = " << onTrack << std::endl;
372  std::cout << "posW = " << posW << std::endl;
373  std::cout << "pOnTrack= " << pOnTrack << std::endl;
374  G4ThreeVector bwp(bck.x(), bck.y(), bck.z());
375  G4ThreeVector fwp(fwd.x(), fwd.y(), fwd.z());
376  G4ThreeVector hitPosition, wirePosition;
377  distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
378  hitPosition, wirePosition);
379  if (m_wireSag) {
380  G4double ywb_sag, ywf_sag;
381  m_cdcgp->getWireSagEffect(CDCGeometryPar::c_Base, layerId, wires[i], wirePosition.z(), ywb_sag, ywf_sag);
382  bwp.setY(ywb_sag);
383  fwp.setY(ywf_sag);
384  distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
385  hitPosition, wirePosition);
386  }
387  std::cout << "line distance= " << distance << std::endl;
388  std::cout << "onTrack= " << hitPosition.x() << " " << hitPosition.y() << " " << hitPosition.z() << std::endl;
389  std::cout << "posW = " << wirePosition.x() << " " << wirePosition.y() << " " << wirePosition.z() << std::endl;
390 #endif
391  }
392  } else { //no magnetic field case
393  // Cal. distance assuming a line track
394  G4ThreeVector bwp(bck.x(), bck.y(), bck.z());
395  G4ThreeVector fwp(fwd.x(), fwd.y(), fwd.z());
396  G4ThreeVector hitPosition, wirePosition;
397  distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
398  hitPosition, wirePosition);
399  if (m_wireSag) {
400  G4double ywb_sag, ywf_sag;
401  m_cdcgp->getWireSagEffect(CDCGeometryPar::c_Base, layerId, wires[i], wirePosition.z(), ywb_sag, ywf_sag);
402  bwp.setY(ywb_sag);
403  fwp.setY(ywf_sag);
404  distance = ClosestApproach(bwp, fwp, posIn / CLHEP::cm, posOut / CLHEP::cm,
405  hitPosition, wirePosition);
406  }
407 
408  onTrack.setX(hitPosition.x());
409  onTrack.setY(hitPosition.y());
410  onTrack.setZ(hitPosition.z());
411  posW.setX(wirePosition.x());
412  posW.setY(wirePosition.y());
413  posW.setZ(wirePosition.z());
414  //tentative setting
415  pOnTrack.setX(0.5 * (momIn.x() + momOut.x()) / CLHEP::GeV);
416  pOnTrack.setY(0.5 * (momIn.y() + momOut.y()) / CLHEP::GeV);
417  pOnTrack.setZ(0.5 * (momIn.z() + momOut.z()) / CLHEP::GeV);
418  } //end of magneticfiled on or off
419 
420 #if defined(CDC_DEBUG)
421  std::cout << "af distance= " << distance << std::endl;
422  std::cout << "onTrack = " << onTrack << std::endl;
423  std::cout << "posW = " << posW << std::endl;
424  std::cout << "pOnTrack = " << pOnTrack << std::endl;
425  if (distance > 2.4) {
426  std::cout << "toolargedriftl" << std::endl;
427  }
428 #endif
429  distance *= CLHEP::cm; onTrack *= CLHEP::cm; posW *= CLHEP::cm;
430  pOnTrack *= CLHEP::GeV;
431 
432  G4ThreeVector posTrack(onTrack.x(), onTrack.y(), onTrack.z());
433  G4ThreeVector mom(pOnTrack.x(), pOnTrack.y(), pOnTrack.z());
434 
435  const TVector3 tPosW(posW.x(), posW.y(), posW.z());
436  const TVector3 tPosTrack(posTrack.x(), posTrack.y(), posTrack.z());
437  const TVector3 tMom(mom.x(), mom.y(), mom.z());
438  G4int lr = m_cdcgp->getOldLeftRight(tPosW, tPosTrack, tMom);
439  G4int newLrRaw = m_cdcgp->getNewLeftRightRaw(tPosW, tPosTrack, tMom);
440  // if(abs(pid) == 11) {
441  // std::cout <<"pid,lr,newLrRaw 4electron= " << pid <<" "<< lr <<" "<< newLrRaw << std::endl;
442  // }
443  G4int newLr = newLrRaw; //to be modified in EndOfEvent
444 
445  if (nWires == 1) {
446 
447  // saveSimHit(layerId, wires[i], trackID, pid, distance, tofBefore, edep, s_in_layer * cm, momIn, posW, posIn, posOut, posTrack, lr, newLrRaw, newLr, speed);
448  saveSimHit(layerId, wires[i], trackID, pid, distance, tofBefore, edep, s_in_layer * CLHEP::cm, pOnTrack, posW, posIn, posOut,
449  posTrack, lr, newLrRaw, newLr, speed, hitWeight);
450 #if defined(CDC_DEBUG)
451  std::cout << "saveSimHit" << std::endl;
452  std::cout << "momIn = " << momIn << std::endl;
453  std::cout << "pOnTrack = " << pOnTrack << std::endl;
454 #endif
455 
456  } else {
457 
458  G4int cel1 = wires[i] + 1;
459  G4int cel2 = cel1;
460  if (i + 1 <= nWires - 1) {
461  cel2 = wires[i + 1] + 1;
462  }
463  const G4double s2 = t.GetTrackLength() / CLHEP::cm; //at post-step
464  G4double s1 = (s2 - s_in_layer); //at pre-step; varied later
465  G4ThreeVector din = momIn;
466  if (din.mag() != 0.) din /= momIn.mag();
467 
468  G4double vent[6] = {posIn.x() / CLHEP::cm, posIn.y() / CLHEP::cm, posIn.z() / CLHEP::cm, din.x(), din.y(), din.z()};
469 
470  G4ThreeVector dot(momOut.x(), momOut.y(), momOut.z());
471  if (dot.mag() != 0.) {
472  dot /= dot.mag();
473  } else {
474  // Flight-direction is needed to set even when a particle stops
475  dot = din;
476  }
477 
478  G4double vext[6] = {posOut.x() / CLHEP::cm, posOut.y() / CLHEP::cm, posOut.z() / CLHEP::cm, dot.x(), dot.y(), dot.z()};
479 
480  if (i > 0) {
481  for (int j = 0; j < 6; ++j) vent[j] = xint[j];
482  s1 = sint;
483  }
484 
485  // const G4int ic(3); // cubic approximation of the track
486  G4int flag(0);
487  G4double edep_in_cell(0.);
488  G4double eLossInCell(0.);
489 
490  if (cel1 != cel2) {
491 #if defined(CDC_DEBUG)
492  std::cout << "layerId,cel1,cel2= " << layerId << " " << cel1 << " " << cel2 << std::endl;
493  std::cout << "vent= " << vent[0] << " " << vent[1] << " " << vent[2] << " " << vent[3] << " " << vent[4] << " " << vent[5] <<
494  std::endl;
495  std::cout << "vext= " << vext[0] << " " << vext[1] << " " << vext[2] << " " << vext[3] << " " << vext[4] << " " << vext[5] <<
496  std::endl;
497  std::cout << "s1,s2,ic= " << s1 << " " << s2 << " " << ic << std::endl;
498 #endif
499  CellBound(layerId, cel1, cel2, vent, vext, s1, s2, xint, sint, flag);
500 #if defined(CDC_DEBUG)
501  std::cout << "flag,sint= " << flag << " " << sint << std::endl;
502  std::cout << "xint= " << xint[0] << " " << xint[1] << " " << xint[2] << " " << xint[3] << " " << xint[4] << " " << xint[5] <<
503  std::endl;
504 #endif
505 
506  const G4double test = (sint - s1) / s_in_layer;
507  if (test < 0. || test > 1.) {
508  B2WARNING("CDCSensitiveDetector: Strange path length: " << "s1= " << s1 << " sint= " << sint << " s_in_layer= " << s_in_layer <<
509  " test= " << test);
510  }
511  edep_in_cell = edep * std::abs((sint - s1)) / s_in_layer;
512 
513  const G4ThreeVector x_In(vent[0]*CLHEP::cm, vent[1]*CLHEP::cm, vent[2]*CLHEP::cm);
514  const G4ThreeVector x_Out(xint[0]*CLHEP::cm, xint[1]*CLHEP::cm, xint[2]*CLHEP::cm);
515  const G4ThreeVector p_In(momBefore * vent[3], momBefore * vent[4], momBefore * vent[5]);
516 
517  // 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);
518  saveSimHit(layerId, wires[i], trackID, pid, distance, tofBefore, edep_in_cell, std::abs((sint - s1)) * CLHEP::cm, pOnTrack, posW,
519  x_In, x_Out,
520  posTrack, lr, newLrRaw, newLr, speed, hitWeight);
521 #if defined(CDC_DEBUG)
522  std::cout << "saveSimHit" << std::endl;
523  std::cout << "p_In = " << p_In << std::endl;
524  std::cout << "pOnTrack= " << pOnTrack << std::endl;
525 #endif
526  tofBefore += (sint - s1) / speedInCmPerNs;
527  eLossInCell = eLoss * (sint - s1) / s_in_layer;
528  kinEnergyBefore -= eLossInCell;
529  if (kinEnergyBefore >= 0.) {
530  momBefore = sqrt(kinEnergyBefore * (kinEnergyBefore + 2.*mass));
531  } else {
532  B2WARNING("CDCSensitiveDetector: Kinetic Energy < 0.");
533  momBefore = 0.;
534  }
535 
536  } else { //the particle exits
537 
538  const G4double test = (s2 - sint) / s_in_layer;
539  if (test < 0. || test > 1.) {
540  B2WARNING("CDCSensitiveDetector: Strange path length: " << "s2= " << s2 << " sint= " << sint << " s_in_layer= " << s_in_layer <<
541  " test= " << test);
542  }
543  edep_in_cell = edep * std::abs((s2 - sint)) / s_in_layer;
544 
545  const G4ThreeVector x_In(vent[0]*CLHEP::cm, vent[1]*CLHEP::cm, vent[2]*CLHEP::cm);
546  const G4ThreeVector p_In(momBefore * vent[3], momBefore * vent[4], momBefore * vent[5]);
547 
548  // 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);
549  saveSimHit(layerId, wires[i], trackID, pid, distance, tofBefore, edep_in_cell, std::abs((s2 - sint)) * CLHEP::cm, pOnTrack, posW,
550  x_In,
551  posOut, posTrack, lr, newLrRaw, newLr, speed, hitWeight);
552 #if defined(CDC_DEBUG)
553  std::cout << "saveSimHit" << std::endl;
554  std::cout << "p_In = " << p_In << std::endl;
555  std::cout << "pOnTrack= " << pOnTrack << std::endl;
556 #endif
557  }
558  }
559  //setSeenInDetectorFlag(aStep, MCParticle::c_SeenInCDC);
560 
564  //StoreArray<Relation> mcPartToSimHits(getRelationCollectionName());
565  //StoreArray<MCParticle> mcPartArray(DEFAULT_MCPARTICLES);
566  //if (saveIndex < 0) {B2FATAL("SimHit wasn't saved despite charge != 0");}
567  //StoreArray<CDCSimHit> cdcArray(DEFAULT_CDCSIMHITS);
568 
569  //new(mcPartToSimHits->AddrAt(saveIndex)) Relation(mcPartArray, cdcArray, trackID, saveIndex);
570 
571  } //end of wire loop
572 
573  return true;
574  }
575 
576  /*
577  void CDCSensitiveDetector::BeginOfEvent(G4HCofThisEvent*)
578  {
579  std::cout <<"CDCSensitiveDetector::BeginOfEvent callded." << std::endl;
580  }
581  */
582 
583  void CDCSensitiveDetector::EndOfEvent(G4HCofThisEvent*)
584  {
585  setModifiedLeftRightFlag();
586  }
587 
588  void
589  CDCSensitiveDetector::saveSimHit(const G4int layerId,
590  const G4int wireId,
591  const G4int trackID,
592  const G4int pid,
593  const G4double distance,
594  const G4double tof,
595  const G4double edep,
596  const G4double stepLength,
597  const G4ThreeVector& mom,
598  const G4ThreeVector& posW,
599  const G4ThreeVector& posIn,
600  const G4ThreeVector& posOut,
601  const G4ThreeVector& posTrack,
602  const G4int lr,
603  const G4int newLrRaw,
604  const G4int newLr,
605  const G4double speed,
606  const G4double hitWeight)
607  {
608 
609  // Discard the hit below Edep_th
610  // if (edep <= m_thresholdEnergyDeposit) return 0;
611  if (edep <= m_thresholdEnergyDeposit) return;
612 
613  //compute tof at the closest point; linear approx.
614  const G4double sign = (posTrack - posIn).dot(mom) < 0. ? -1. : 1.;
615  const G4double CorrectTof = tof + sign * (posTrack - posIn).mag() / speed;
616  // if (sign < 0.) std::cout <<"deltatof= "<< sign * (posTrack - posIn).mag() / speed << std::endl;
617 #if defined(CDC_DEBUG)
618  std::cout << "posIn= " << posIn.x() << " " << posIn.y() << " " << posIn.z() << std::endl;
619  std::cout << "posOut= " << posOut.x() << " " << posOut.y() << " " << posOut.z() << std::endl;
620  std::cout << "posTrack= " << posTrack.x() << " " << posTrack.y() << " " << posTrack.z() << std::endl;
621  std::cout << "posW= " << posW.x() << " " << posW.y() << " " << posW.z() << std::endl;
622  std::cout << "tof = " << tof << std::endl;
623  std::cout << "deltaTof = " << (posTrack - posIn).mag() / speed << std::endl;
624  std::cout << "CorrectTof= " << CorrectTof << std::endl;
625  if (CorrectTof > 95) {
626  std::cout << "toolargecorrecttof" << std::endl;
627  }
628 #endif
629 
630  StoreArray<MCParticle> mcParticles;
631 
632  StoreArray<CDCSimHit> cdcArray;
633 
634  RelationArray cdcSimHitRel(mcParticles, cdcArray);
635 
636  m_hitNumber = cdcArray.getEntries();
637 
638  CDCSimHit* simHit = cdcArray.appendNew();
639 
640  simHit->setWireID(layerId, wireId);
641  simHit->setTrackId(trackID);
642  simHit->setPDGCode(pid);
643  simHit->setDriftLength(distance / CLHEP::cm);
644  simHit->setFlightTime(CorrectTof / CLHEP::ns);
645  simHit->setGlobalTime(CorrectTof / CLHEP::ns);
646  simHit->setEnergyDep(edep / CLHEP::GeV);
647  simHit->setStepLength(stepLength / CLHEP::cm);
648  TVector3 momentum(mom.getX() / CLHEP::GeV, mom.getY() / CLHEP::GeV, mom.getZ() / CLHEP::GeV);
649  simHit->setMomentum(momentum);
650  TVector3 posWire(posW.getX() / CLHEP::cm, posW.getY() / CLHEP::cm, posW.getZ() / CLHEP::cm);
651  simHit->setPosWire(posWire);
652  TVector3 positionIn(posIn.getX() / CLHEP::cm, posIn.getY() / CLHEP::cm, posIn.getZ() / CLHEP::cm);
653  simHit->setPosIn(positionIn);
654  TVector3 positionOut(posOut.getX() / CLHEP::cm, posOut.getY() / CLHEP::cm, posOut.getZ() / CLHEP::cm);
655  simHit->setPosOut(positionOut);
656  TVector3 positionTrack(posTrack.getX() / CLHEP::cm, posTrack.getY() / CLHEP::cm, posTrack.getZ() / CLHEP::cm);
657  simHit->setPosTrack(positionTrack);
658  simHit->setPosFlag(lr);
659  simHit->setLeftRightPassageRaw(newLrRaw);
660  simHit->setLeftRightPassage(newLr);
661 #if defined(CDC_DEBUG)
662  std::cout << "sensitived,oldlr,newlrRaw,newlr= " << lr << " " << newLrRaw << " " << newLr << std::endl;
663 #endif
664 
665  B2DEBUG(150, "HitNumber: " << m_hitNumber);
666  if (m_modifiedLeftRightFlag) {
667  //N.B. Negative hitWeight is allowed intentionally here; all weights are to be reset to positive in EndOfEvent
668  cdcSimHitRel.add(trackID, m_hitNumber, hitWeight);
669  } else {
670  cdcSimHitRel.add(trackID, m_hitNumber);
671  }
672 
673  // if (hitWeight > 0) m_nPosHits++;
674  // if (hitWeight < 0) m_nNegHits++;
675  // std::cout <<"trackID,HitNumber,weight,driftL,edep= "<< trackID <<" "<< m_hitNumber <<" "<< hitWeight <<" "<< distance <<" "<< edep << std::endl;
676  // return (m_hitNumber);
677  }
678 
679 
680  /*
681  void SensitiveDetector::AddbgOne(bool doit) {
682  Belle::Datcdc_olhit_Manager& olhitmgr=Belle::Datcdc_olhit_Manager::get_manager();
683  if(doit) {
684  for( int i=0;i<olhitmgr.count();i++ ){
685  Belle::Datcdc_olhit& h = olhitmgr[i];
686  }
687  dout(Debugout::B2INFO,"SensitiveDetector")
688  << "AddbgOne " << olhitmgr.size()
689  << std::endl;
690  }
691  olhitmgr.remove();
692  }
693  */
694 
695  void
696  CDCSensitiveDetector::CellBound(const G4int layerId,
697  const G4int ic1,
698  const G4int ic2,
699  const G4double venter[6],
700  const G4double vexit[6],
701  const G4double s1, const G4double s2,
702  G4double xint[6],
703  G4double& sint, G4int& iflag)
704  {
705  //---------------------------------------------------------------------------
706  // (Purpose)
707  // calculate an intersection of track with cell boundary.
708  //
709  // (Relations)
710  // Calls GCUBS
711  //
712  // (Arguments)
713  // input
714  // ic1 serial cell# (start w/ one) of entrance.
715  // ic2 serial cell# (start w/ one) of exit.
716  // venter(6) (x,y,z,px/p,py/p,pz/p) at entrance.
717  // vexit(6) (x,y,z,px/p,py/p,pz/p) at exit.
718  // s1 track length at entrance.
719  // s2 track length at exit.
720  // output
721  // xint(6) (x,y,z,px/p,py/p,pz/p) at intersection of cell boundary.
722  // sint track length at intersection of cell boundary.
723  // iflag return code.
724  //
725  // N.B.(TODO ?) CDC misalignment wrt Belle2 coordinate system is ignored
726  // when calculating the cell-boundary assuming misalign. is small.
727  //--------------------------------------------------------------------------
728 
729  G4double div = m_cdcgp->nWiresInLayer(layerId);
730 
731  //Check if s1, s2, ic1 and ic2 are ok
732  if (s1 >= s2) {
733  B2ERROR("CDCSensitiveDetector: s1(=" << s1 << ") > s2(=" << s2 << ")");
734  }
735  if (std::abs(ic1 - ic2) != 1) {
736  if (ic1 == 1 && ic2 == div) {
737  } else if (ic1 == div && ic2 == 1) {
738  } else {
739  B2ERROR("CDCSensitiveDetector: |ic1 - ic2| != 1 in CellBound; " << "ic1=" << ic1 << " " << "ic2=" << ic2);
740  }
741  }
742 
743  //get wire positions for the entrance cell
744  G4double xwb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).x();
745  G4double ywb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).y();
746  G4double zwb = (m_cdcgp->wireBackwardPosition(layerId, ic1 - 1)).z();
747  G4double xwf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).x();
748  G4double ywf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).y();
749  G4double zwf = (m_cdcgp->wireForwardPosition(layerId, ic1 - 1)).z();
750 
751  /*
752  G4double pathl = sqrt((vexit[0] - venter[0]) * (vexit[0] - venter[0])
753  + (vexit[1] - venter[1]) * (vexit[1] - venter[1])
754  + (vexit[2] - venter[2]) * (vexit[2] - venter[2]));
755  std::cout << "app pathl= " << pathl << std::endl;
756  G4double dot = venter[3] * vexit[3] + venter[4] * vexit[4];
757  dot /= sqrt(venter[3] * venter[3] + venter[4] * venter[4]);
758  dot /= sqrt( vexit[3] * vexit[3] + vexit[4] * vexit[4]);
759  if (dot < 0.) std::cout <<"negativedot= " << dot << std::endl;
760  */
761 
762  //copy arrays
763  G4double xx1[6], xx2[6];
764  for (int i = 0; i < 6; ++i) {
765  xx1[i] = venter[i];
766  xx2[i] = vexit [i];
767  }
768 
769  //calculate the field wire position betw. cell#1 and #2
770  G4double psi = double(ic2 - ic1) * CLHEP::pi / div;
771  if (ic1 == 1 && ic2 == div) {
772  psi = -CLHEP::pi / div;
773  } else if (ic1 == div && ic2 == 1) {
774  psi = CLHEP::pi / div;
775  }
776  G4double cospsi = cos(psi);
777  G4double sinpsi = sin(psi);
778 
779  G4double xfwb = cospsi * xwb - sinpsi * ywb;
780  G4double yfwb = sinpsi * xwb + cospsi * ywb;
781  G4double xfwf = cospsi * xwf - sinpsi * ywf;
782  G4double yfwf = sinpsi * xwf + cospsi * ywf;
783  G4double zfwb = zwb;
784  G4double zfwf = zwf;
785 
786  //prepare quantities related to the cell-boundary
787  G4double vx = xfwf - xfwb;
788  G4double vy = yfwf - yfwb;
789  G4double vz = zfwf - zfwb;
790  G4double vv = sqrt(vx * vx + vy * vy + vz * vz);
791  vx /= vv; vy /= vv; vz /= vv;
792 
793  //translate to make the cubic description easier
794  G4double shiftx = (xx1[0] + xx2[0]) * 0.5;
795  G4double shifty = (xx1[1] + xx2[1]) * 0.5;
796  G4double shiftz = (xx1[2] + xx2[2]) * 0.5;
797  G4double shifts = (s1 + s2) * 0.5;
798  G4double xshft = xx1[0] - shiftx;
799  G4double yshft = xx1[1] - shifty;
800  G4double zshft = xx1[2] - shiftz;
801  G4double sshft = s1 - shifts;
802 
803  //approximate the trajectroy by cubic curves
804  G4double pabs1 = sqrt(xx1[3] * xx1[3] + xx1[4] * xx1[4] + xx1[5] * xx1[5]);
805  G4double pabs2 = sqrt(xx2[3] * xx2[3] + xx2[4] * xx2[4] + xx2[5] * xx2[5]);
806  // std::cout << "pabs1,2= " << pabs1 <<" "<< pabs2 << std::endl;
807 
808  G4double a[4] = {0.}, b[4] = {0.}, c[4] = {0.};
809 
810  if (m_magneticField) {
811  GCUBS(sshft, xshft, xx1[3] / pabs1, xx2[3] / pabs2, a);
812  GCUBS(sshft, yshft, xx1[4] / pabs1, xx2[4] / pabs2, b);
813  GCUBS(sshft, zshft, xx1[5] / pabs1, xx2[5] / pabs2, c);
814  // std::cout <<"a= " << a[0] <<" "<< a[1] <<" "<< a[2] <<" "<< a[3] << std::endl;
815  // std::cout <<"b= " << b[0] <<" "<< b[1] <<" "<< b[2] <<" "<< b[3] << std::endl;
816  // std::cout <<"c= " << c[0] <<" "<< c[1] <<" "<< c[2] <<" "<< c[3] << std::endl;
817  } else {
818  //n.b. following is really better ?
819  a[1] = xshft / sshft;
820  b[1] = yshft / sshft;
821  c[1] = zshft / sshft;
822  }
823 
824  //calculate an int. point betw. the trajectory and the cell-boundary
825  G4double stry(0.), xtry(0.), ytry(0.), ztry(0.);
826  G4double beta(0.), xfw(0.), yfw(0.);
827  G4double sphi(0.), cphi(0.), dphil(0.), dphih(0.);
828  const G4int maxTrials = 100;
829  const G4double eps = 5.e-4;
830  G4double sl = sshft; // negative value
831  G4double sh = -sshft; // positive value
832  G4int i = 0;
833 
834  //set initial value (dphil) for the 1st iteration
835  stry = sl;
836  xtry = shiftx + a[0] + stry * (a[1] + stry * (a[2] + stry * a[3]));
837  ytry = shifty + b[0] + stry * (b[1] + stry * (b[2] + stry * b[3]));
838  ztry = shiftz + c[0] + stry * (c[1] + stry * (c[2] + stry * c[3]));
839  beta = (ztry - zfwb) / vz;
840  xfw = xfwb + beta * vx;
841  yfw = yfwb + beta * vy;
842  sphi = (xtry * yfw - ytry * xfw);
843  cphi = (xtry * xfw + ytry * yfw);
844  dphil = atan2(sphi, cphi); //n.b. no need to conv. to dphi...
845 
846  iflag = 1;
847 
848  while (((sh - sl) > eps) && (i < maxTrials)) {
849  stry = 0.5 * (sl + sh);
850  xtry = shiftx + a[0] + stry * (a[1] + stry * (a[2] + stry * a[3]));
851  ytry = shifty + b[0] + stry * (b[1] + stry * (b[2] + stry * b[3]));
852  ztry = shiftz + c[0] + stry * (c[1] + stry * (c[2] + stry * c[3]));
853  beta = (ztry - zfwb) / vz;
854  xfw = xfwb + beta * vx;
855  yfw = yfwb + beta * vy;
856 
857  sphi = (xtry * yfw - ytry * xfw);
858  cphi = (xtry * xfw + ytry * yfw);
859  dphih = atan2(sphi, cphi); //n.b. no need to conv. to dphi...
860 
861  if (dphil * dphih > 0.) {
862  sl = stry;
863  } else {
864  sh = stry;
865  }
866  ++i;
867  }
868 
869  // std::cout << "itry= " << i << std::endl;
870  if (i >= maxTrials - 1) {
871  iflag = 0;
872  B2WARNING("CDCSensitiveDetector: No intersection ?");
873 
874  /* G4double ds = 1.e-4;
875  G4int imax = (s2 - s1) / ds + 1;
876  G4double rdphimin = DBL_MAX;
877 
878  for (i=0; i <= imax; ++i) {
879  stry = sshft + i * ds;
880  xtry = shiftx + a[0] + stry * (a[1] + stry * (a[2] + stry * a[3]));
881  ytry = shifty + b[0] + stry * (b[1] + stry * (b[2] + stry * b[3]));
882  ztry = shiftz + c[0] + stry * (c[1] + stry * (c[2] + stry * c[3]));
883  beta = (ztry - zfwb) / vz;
884  xfw = xfwb + beta * vx;
885  yfw = yfwb + beta * vy;
886 
887  sphi = (xtry * yfw - ytry * xfw);
888  cphi = (xtry * xfw + ytry * yfw);
889  dphi = atan2(sphi, cphi);
890  rdphi = sqrt(xfw * xfw + yfw * yfw) * dphi;
891 
892  if ( std::abs(rdphi) < std::abs(rdphimin)) {
893  rdphimin = rdphi;
894  imin = i;
895  }
896  }
897  */
898  }
899  // sint = sshft + imin * ds;
900  sint = stry;
901 
902  // std::cout <<"i,dphil,dphih,sint= " << i <<" "<< dphil <<" "<< dphih <<" "<< sint << std::endl;
903  //get the trajectory at the int. point
904  xint[0] = a[0] + sint * (a[1] + sint * (a[2] + sint * a[3]));
905  xint[1] = b[0] + sint * (b[1] + sint * (b[2] + sint * b[3]));
906  xint[2] = c[0] + sint * (c[1] + sint * (c[2] + sint * c[3]));
907  xint[3] = a[1] + sint * (2. * a[2] + 3. * sint * a[3]);
908  xint[4] = b[1] + sint * (2. * b[2] + 3. * sint * b[3]);
909  xint[5] = c[1] + sint * (2. * c[2] + 3. * sint * c[3]);
910 
911  //translate back to the lab. frame
912  xint[0] += shiftx;
913  xint[1] += shifty;
914  xint[2] += shiftz;
915  sint += shifts;
916  /*
917  std::cout <<"s1,s2,sint= " << s1 <<" "<< s2 <<" "<< sint << std::endl;
918  std::cout <<" xx1= " << xx1[0] <<" "<< xx1[1] <<" "<< xx1[2] << std::endl;
919  std::cout <<" xx2= " << xx2[0] <<" "<< xx2[1] <<" "<< xx2[2] << std::endl;
920  std::cout <<"xint= " << xint[0] <<" "<< xint[1] <<" "<< xint[2] << std::endl;
921  */
922 
923  /* if (((xx1[0] <= xint[0] && xint[0] <= xx2[0]) ||
924  (xx2[0] <= xint[0] && xint[0] <= xx1[0])) &&
925  ((xx1[1] <= xint[1] && xint[1] <= xx2[1]) ||
926  (xx2[1] <= xint[1] && xint[1] <= xx1[1])) &&
927  ((xx1[2] <= xint[2] && xint[2] <= xx2[2]) ||
928  (xx2[2] <= xint[2] && xint[2] <= xx1[2])) &&
929  (s1 <= sint && sint <= s2)) {
930  } else {
931  std::cout << "strangeinttersection" << std::endl;
932  }
933  */
934  //re-normalize to one since abs=1 is not guearanteed in the cubic approx.
935  G4double p = sqrt(xint[3] * xint[3] + xint[4] * xint[4] + xint[5] * xint[5]);
936  xint[3] /= p; xint[4] /= p; xint[5] /= p;
937  // std::cout << "norm= " << p << std::endl;
938  // std::cout <<"s1,s2,sint= " << s1 <<" "<< s2 <<" "<< sint << std::endl;
939  // std::cout <<"xint= " << xint[0] <<" "<< xint[1] <<" "<< xint[2] << std::endl;
940  // std::cout <<"xint= " << xint[3] <<" "<< xint[4] <<" "<< xint[5] << std::endl;
941  }
942 
943  void CDCSensitiveDetector::GCUBS(const G4double x, const G4double y, const G4double d1, const G4double d2, G4double a[4])
944  {
945  //Original: GCUBS in Geant3
946  // ******************************************************************
947  // * *
948  // * Calculates a cubic through P1,(X,Y),(-X,-Y),P2 *
949  // * Y=A(1)+A(2)*X+A(3)*X**2+A(4)*X**3 *
950  // * The coordinate system is assumed to be the cms system *
951  // * of P1,P2. *
952  // * d1(2): directional cosine at P1(2). *
953  // * *
954  // * ==>Called by : GIPLAN,GICYL *
955  // * Author H.Boerner ********* *
956  // * *
957  // ******************************************************************
958 
959  G4double fact(0);
960 
961  if (x == 0.) goto L10;
962 
963  fact = (d1 - d2) * 0.25;
964  a[0] = - 1. * fact * x;
965  a[2] = fact / x;
966  a[1] = (6. * y - (d1 + d2) * x) / (4. * x);
967  a[3] = ((d1 + d2) * x - 2.*y) / (4.*x * x * x);
968  return;
969 
970 L10:
971  a[0] = 0.;
972  a[1] = 1.;
973  a[2] = 0.;
974  a[3] = 0.;
975  }
976 
977  void
978  CDCSensitiveDetector::for_Rotat(const G4double bfld[3])
979  {
980  //Calculates a rotation matrix in advance at a local position in lab.
981  //The rotation is done about the coord. origin; lab.-frame to B-field
982  //frame in which only Bz-comp. is non-zero.
983  //~dead copy of gsim_cdc_for_rotat.F in gsim-cdc for Belle (for tentaive use)
984 
985  if (m_nonUniformField == 0) return;
986 
987  G4double bx, by, bz;
988  bx = bfld[0];
989  by = bfld[1];
990  bz = bfld[2];
991 
992  //cal. rotation matrix
993  G4double bxz, bfield;
994  bxz = bx * bx + bz * bz;
995  bfield = bxz + by * by;
996  bxz = sqrt(bxz);
997  bfield = sqrt(bfield);
998 
999  m_brot[0][0] = bz / bxz;
1000  m_brot[1][0] = 0.;
1001  m_brot[2][0] = -bx / bxz;
1002  m_brot[0][1] = -by * bx / bxz / bfield;
1003  m_brot[1][1] = bxz / bfield;
1004  m_brot[2][1] = -by * bz / bxz / bfield;
1005  m_brot[0][2] = bx / bfield;
1006  m_brot[1][2] = by / bfield;
1007  m_brot[2][2] = bz / bfield;
1008 
1009  return;
1010 
1011  }
1012 
1013  void
1014  CDCSensitiveDetector::Rotat(G4double& x, G4double& y, G4double& z,
1015  const int mode)
1016  {
1017  //Translates (x,y,z) in lab. to (x,y,z) in B-field frame (mode=1), or reverse
1018  // translation (mode=-1).
1019  //~dead copy (for tentaive use) of gsim_cdc_rotat/irotat.F in gsim-cdc
1020  //for Belle
1021 
1022  if (m_nonUniformField == 0) return;
1023 
1024  G4double x0(x), y0(y), z0(z);
1025 
1026  if (mode == 1) {
1027  x = m_brot[0][0] * x0 + m_brot[1][0] * y0 + m_brot[2][0] * z0;
1028  y = m_brot[0][1] * x0 + m_brot[1][1] * y0 + m_brot[2][1] * z0;
1029  z = m_brot[0][2] * x0 + m_brot[1][2] * y0 + m_brot[2][2] * z0;
1030  } else if (mode == -1) {
1031  x = m_brot[0][0] * x0 + m_brot[0][1] * y0 + m_brot[0][2] * z0;
1032  y = m_brot[1][0] * x0 + m_brot[1][1] * y0 + m_brot[1][2] * z0;
1033  z = m_brot[2][0] * x0 + m_brot[2][1] * y0 + m_brot[2][2] * z0;
1034  } else {
1035  //B2ERROR("SensitiveDetector " <<"invalid mode " << mode << "specifed");
1036  }
1037  return;
1038 
1039  }
1040 
1041  void
1042  CDCSensitiveDetector::Rotat(G4double x[3], const int mode)
1043  {
1044  //Translates (x,y,z) in lab. to (x,y,z) in B-field frame (mode=1), or reverse
1045  // translation (mode=-1).
1046  //~dead copy (for tentaive use) of gsim_cdc_rotat/irotat.F in gsim-cdc
1047  //for Belle
1048 
1049  if (m_nonUniformField == 0) return;
1050 
1051  G4double x0(x[0]), y0(x[1]), z0(x[2]);
1052 
1053  if (mode == 1) {
1054  x[0] = m_brot[0][0] * x0 + m_brot[1][0] * y0 + m_brot[2][0] * z0;
1055  x[1] = m_brot[0][1] * x0 + m_brot[1][1] * y0 + m_brot[2][1] * z0;
1056  x[2] = m_brot[0][2] * x0 + m_brot[1][2] * y0 + m_brot[2][2] * z0;
1057  } else if (mode == -1) {
1058  x[0] = m_brot[0][0] * x0 + m_brot[0][1] * y0 + m_brot[0][2] * z0;
1059  x[1] = m_brot[1][0] * x0 + m_brot[1][1] * y0 + m_brot[1][2] * z0;
1060  x[2] = m_brot[2][0] * x0 + m_brot[2][1] * y0 + m_brot[2][2] * z0;
1061  } else {
1062  //B2ERROR("SensitiveDetector " <<"invalid mode " << mode << "specifed");
1063  }
1064  return;
1065 
1066  }
1067 
1068  void
1069  CDCSensitiveDetector::HELWIR(const G4double xwb4, const G4double ywb4,
1070  const G4double zwb4,
1071  const G4double xwf4, const G4double ywf4,
1072  const G4double zwf4,
1073  const G4double xp, const G4double yp,
1074  const G4double zp,
1075  const G4double px, const G4double py,
1076  const G4double pz,
1077  const G4double B_kG[3],
1078  const G4double charge, const G4int ntryMax,
1079  G4double& distance,
1080  G4double q2[3], G4double q1[3],
1081  G4double q3[3],
1082  G4int& ntry)
1083  {
1084  //~dead copy of gsim_cdc_hit.F in gsim-cdc for Belle (for tentaive use)
1085  // ---------------------------------------------------------------------
1086  // Purpose : Calculate closest points between helix and wire.
1087  //
1088  // Input
1089  // xwb4 : x of wire at backward endplate in lab.
1090  // ywb4 : y of wire at backward endplate "
1091  // zwb4 : z of wire at backward endplate "
1092  // xwf4 : x of wire at forward endplate "
1093  // ywf4 : y of wire at forward endplate "
1094  // zwf4 : z of wire at forward endplate "
1095  //
1096  // Output
1097  // q2(1) : x of wire at closest point in lab.
1098  // q2(2) : y of wire at closest point "
1099  // q2(3) : z of wire at closest point "
1100  // q1(1) : x of helix at closest point "
1101  // q1(2) : y of helix at closest point "
1102  // q1(3) : z of helix at closest point "
1103  // ntry :
1104  // ---------------------------------------------------------------------
1105 
1106  const G4int ndim = 3;
1107  const G4double delta = 1.e-5;
1108 
1109 
1110  G4double xwb, ywb, zwb, xwf, ywf, zwf;
1111  G4double xw, yw, zw, xh, yh, zh, pxh, pyh, pzh;
1112  G4double fi, fi_corr;
1113 
1114  G4double dr, fi0, cpa, dz, tanl;
1115  G4double x0, y0, z0;
1116  // "chrg" removed by M. U. June, 2nd, 2013
1117  // G4double xc, yc, r, chrg;
1118  G4double xc, yc, r;
1119  G4double xwm, ywm;
1120  G4double sinfi0, cosfi0, sinfi0fi, cosfi0fi;
1121 
1122  G4double vx, vy, vz, vv, cx, cy, cz, tt[3][3];
1123  G4double tmp[3];
1124 
1125  G4double xx[3], dxx[3], ddxx[3], pp[3];
1126  G4double xxtdxx, dxxtdxx, xxtddxx;
1127 
1128 
1129  G4double fst = 0.0;
1130  G4double f, fderiv, deltafi, fact, eval;
1131  G4double dx1, dy1, dx2, dy2, crs, dot;
1132 
1133  G4int iflg;
1134 
1135  //set parameters
1136  xwb = xwb4; ywb = ywb4; zwb = zwb4;
1137  xwf = xwf4; ywf = ywf4; zwf = zwf4;
1138 
1139  G4double xxx(xp), yyy(yp), zzz(zp);
1140  G4double pxx(px), pyy(py), pzz(pz);
1141 
1142  //rotate z-axis to be parallel to B-field in case of non-uniform B
1143  Rotat(xwb, ywb, zwb, 1);
1144  Rotat(xwf, ywf, zwf, 1);
1145  Rotat(xxx, yyy, zzz, 1);
1146  Rotat(pxx, pyy, pzz, 1);
1147 
1148  G4double a[8] = {0.};
1149  G4double pt = sqrt(pxx * pxx + pyy * pyy);
1150  a[1] = atan2(-pxx, pyy);
1151  a[2] = charge / pt;
1152  a[4] = pzz / pt;
1153  a[5] = xxx; a[6] = yyy; a[7] = zzz;
1154 
1155  //calculate unit direction vector of the sense wire
1156  vx = xwf - xwb; vy = ywf - ywb; vz = zwf - zwb;
1157  vv = sqrt(vx * vx + vy * vy + vz * vz);
1158  vx /= vv; vy /= vv; vz /= vv;
1159 
1160  //flag for distingushing between stereo and axial wire
1161  iflg = 0;
1162  if (vx == 0. && vy == 0.) iflg = 1;
1163  // std::cout << "iflg= " << iflg << std::endl;
1164  //write(6,*) ' hlx2wir ', xwb, ywb, zwb, vx, vy, vz
1165 
1166  //calculate coefficients of f
1167  cx = xwb - vx * (vx * xwb + vy * ywb + vz * zwb);
1168  cy = ywb - vy * (vx * xwb + vy * ywb + vz * zwb);
1169  cz = zwb - vz * (vx * xwb + vy * ywb + vz * zwb);
1170 
1171  //calculate tensor for f
1172  tt[0][0] = vx * vx - 1.; tt[1][0] = vx * vy; tt[2][0] = vx * vz;
1173  tt[0][1] = vy * vx; tt[1][1] = vy * vy - 1.; tt[2][1] = vy * vz;
1174  tt[0][2] = vz * vx; tt[1][2] = vz * vy; tt[2][2] = vz * vz - 1.;
1175 
1176  //set helix parameters
1177  dr = a[0]; fi0 = a[1]; cpa = a[2];
1178  dz = a[3]; tanl = a[4];
1179  x0 = a[5]; y0 = a[6]; z0 = a[7];
1180 
1181  //
1182  // set initial value for phi
1183  //
1184 
1185  xwm = xxx;
1186  ywm = yyy;
1187  //r(cm) = alpha/cpa = alpha * pt(GeV); bfield(kG)
1188  G4double bfield = sqrt(B_kG[0] * B_kG[0] +
1189  B_kG[1] * B_kG[1] +
1190  B_kG[2] * B_kG[2]);
1191  G4double alpha = 1.e4 / 2.99792458 / bfield;
1192  r = alpha / cpa;
1193  cosfi0 = cos(fi0);
1194  sinfi0 = sin(fi0);
1195  xc = x0 + (dr + r) * cosfi0;
1196  yc = y0 + (dr + r) * sinfi0;
1197  dx1 = x0 - xc;
1198  dy1 = y0 - yc;
1199  dx2 = xwm - xc;
1200  dy2 = ywm - yc;
1201  crs = dx1 * dy2 - dy1 * dx2;
1202  dot = dx1 * dx2 + dy1 * dy2;
1203  fi = atan2(crs, dot);
1204 
1205  //begin iterative procedure for newton 's method '
1206  fact = 1.;
1207  ntry = 0;
1208 line1:
1209  ntry += 1;
1210  cosfi0fi = cos(fi0 + fi);
1211  sinfi0fi = sin(fi0 + fi);
1212 
1213  //calculate spatial point Q(x,y,z) along the helix
1214  xx[0] = x0 + dr * cosfi0 + r * (cosfi0 - cosfi0fi);
1215  xx[1] = y0 + dr * sinfi0 + r * (sinfi0 - sinfi0fi);
1216  xx[2] = z0 + dz - r * tanl * fi;
1217  pp[0] = -pt * sinfi0fi;
1218  pp[1] = pt * cosfi0fi;
1219  pp[2] = pt * tanl;
1220 
1221  if (iflg == 1) {
1222  q2[0] = xwb; q2[1] = ywb; q2[2] = xx[2];
1223  q1[0] = xx[0]; q1[1] = xx[1]; q1[2] = xx[2];
1224  q3[0] = pp[0]; q3[1] = pp[1]; q3[2] = pp[2];
1225  //inverse rotation to lab. frame in case of non-uniform B
1226  Rotat(q1, -1);
1227  Rotat(q2, -1);
1228  Rotat(q3, -1);
1229  distance = sqrt((q2[0] - q1[0]) * (q2[0] - q1[0]) +
1230  (q2[1] - q1[1]) * (q2[1] - q1[1]) +
1231  (q2[2] - q1[2]) * (q2[2] - q1[2]));
1232  return;
1233  }
1234 
1235  //calculate direction vector (dx/dphi,dy/dphi,dz/dphi)
1236  //on a point along the helix.
1237  dxx[0] = r * sinfi0fi; dxx[1] = - r * cosfi0fi; dxx[2] = - r * tanl;
1238 
1239  // In order to derive the closest pont between straight line and helix,
1240  // we can put following two conditions:
1241  // (i) A point H(xh,yh,zh) on the helix given should be on
1242  // the plane which is perpendicular to the straight line.
1243  // (ii) A line HW from W(xw,yw,zw) which is a point on the straight
1244  // line to H(xh,yh,zh) should normal to the direction vector
1245  // on the point H.
1246  //
1247  // Thus, we can make a equation from above conditions.
1248  // f(phi) = cx*(dx/dphi) + cy*(dy/dphi) + cz*(dz/dphi)
1249  // + (x,y,z)*tt(i,j)*(dx/dphi,dy/dphi,dz/dphi)
1250  // = 0,
1251  // where
1252  // cx = xwb - vx*( vx*xwb + vy*ywb + vz*zwb )
1253  // cy = ywb - vy*( vx*xwb + vy*ywb + vz*zwb )
1254  // cz = zwb - vz*( vx*xwb + vy*ywb + vz*zwb )
1255  //
1256  // tt(1,1) = vx*vx - 1 tt(1,2) = vx*vy tt(1,3) = vx*vz
1257  // tt(2,1) = vy*vx tt(2,2) = vy*vy - 1 tt(2,3) = vy*vz
1258  // tt(3,1) = vz*vx tt(3,2) = vz*vy tt(3,3) = vz*vz - 1
1259  //
1260  // and the equation of straight line(stereo wire) is written by
1261  // (x,y,z) = (xwb,ywb,zwb) + beta*(vx,vy,vz), beta is free parameter.
1262 
1263  //Now calculate f
1264  Mvopr(ndim, xx, tt, dxx, tmp, 1);
1265  xxtdxx = tmp[0];
1266  f = cx * dxx[0] + cy * dxx[1] + cz * dxx[2] + xxtdxx;
1267  if (std::abs(f) < delta) goto line100;
1268 
1269  //evaluate fitting result and prepare some factor to multiply to 1/derivative
1270  if (ntry > 1) {
1271  eval = (1.0 - 0.25 * fact) * std::abs(fst) - std::abs(f);
1272  if (eval <= 0.) fact *= 0.5;
1273  }
1274 
1275  //calculate derivative of f
1276  ddxx[0] = r * cosfi0fi; ddxx[1] = r * sinfi0fi; ddxx[2] = 0.;
1277 
1278  //Now we have derivative of f
1279  Mvopr(ndim, dxx, tt, dxx, tmp, 1);
1280  dxxtdxx = tmp[0];
1281  Mvopr(ndim, xx, tt, ddxx, tmp, 1);
1282  xxtddxx = tmp[0];
1283  fderiv = cx * ddxx[0] + cy * ddxx[1] + cz * ddxx[2] + dxxtdxx + xxtddxx;
1284  // Commented by M. U. June, 2nd, 2013
1285  // fist = fi;
1286  deltafi = f / fderiv;
1287  fi -= fact * deltafi;
1288  fst = f;
1289 
1290  if (ntry > ntryMax) {
1291  //B2DEBUG(" Exceed max. trials HelWir ");
1292  goto line100;
1293  }
1294  //write(6,*) ntry, fist, deltafi
1295  goto line1;
1296 
1297  //check if zh is btw zwb and zwf; if not, set zh=zwb or zh=zwf.
1298  //dead regions due to feed-throughs should be considered later.
1299 line100:
1300  zh = z0 + dz - r * tanl * fi;
1301  fi_corr = 0.;
1302  if (zh < zwb) fi_corr = (zwb - zh) / (-r * tanl);
1303  if (zh > zwf) fi_corr = (zwf - zh) / (-r * tanl);
1304  fi += fi_corr;
1305 
1306  cosfi0fi = cos(fi0 + fi);
1307  sinfi0fi = sin(fi0 + fi);
1308 
1309  xh = x0 + dr * cosfi0 + r * (cosfi0 - cosfi0fi);
1310  yh = y0 + dr * sinfi0 + r * (sinfi0 - sinfi0fi);
1311  zh = z0 + dz - r * tanl * fi;
1312  pxh = -pt * sinfi0fi;
1313  pyh = pt * cosfi0fi;
1314  pzh = pt * tanl;
1315 
1316  //write(6,*) 'fi_corr, zh, zwb, zwf=', fi_corr, zh, zwb, zwf
1317  //write(6,*) 'zh = ', z0, dz, r, tanl, fi
1318 
1319  zw = vx * vz * xh + vy * vz * yh + vz * vz * zh + zwb - vz * (vx * xwb + vy * ywb + vz * zwb);
1320  xw = xwb + vx * (zw - zwb) / vz;
1321  yw = ywb + vy * (zw - zwb) / vz;
1322 
1323  q2[0] = xw; q2[1] = yw; q2[2] = zw;
1324  q1[0] = xh; q1[1] = yh; q1[2] = zh;
1325  q3[0] = pxh; q3[1] = pyh; q3[2] = pzh;
1326 
1327  //inverse rotation to lab. frame in case of non-uniform B
1328  Rotat(q1, -1);
1329  Rotat(q2, -1);
1330  Rotat(q3, -1);
1331  distance = sqrt((q2[0] - q1[0]) * (q2[0] - q1[0]) +
1332  (q2[1] - q1[1]) * (q2[1] - q1[1]) +
1333  (q2[2] - q1[2]) * (q2[2] - q1[2]));
1334  return;
1335 
1336  }
1337 
1338  void
1339  CDCSensitiveDetector::Mvopr(const G4int ndim, const G4double b[3], const G4double m[3][3],
1340  const G4double a[3], G4double c[3], const G4int mode)
1341  {
1342  //~dead copy of UtilCDC_mvopr in com-cdc for Belle (for tentative use)
1343  //-----------------------------------------------------------------------
1344  // Input
1345  // ndim : dimension
1346  // b(1-ndim) : vector
1347  // m(1-ndim,1-ndim) : matrix
1348  // a(1-ndim) : vector
1349  // c(1-ndim) : vector
1350  // mode : c = m * a for mode=0
1351  // c = b * m * a for mode=1
1352  // Output
1353  // c(1-ndim) : for mode 1, solution is put on c[0]
1354  //-----------------------------------------------------------------------
1355 
1356  if (ndim != 3) {
1357  //B2ERROR("invalid ndim " << ndim << " specified");
1358  return;
1359  }
1360 
1361  for (int i = 0; i < ndim; ++i) c[i] = 0.;
1362  G4double tmp[3];
1363  for (int i = 0; i < ndim; ++i) tmp[i] = 0.;
1364 
1365  if (mode == 0) {
1366  for (int i = 0; i < ndim; ++i) {
1367  for (int j = 0; j < ndim; ++j) {
1368  c[i] += m[j][i] * a[j];
1369  }
1370  }
1371  return;
1372  } else if (mode == 1) {
1373  for (int i = 0; i < ndim; ++i) {
1374  for (int j = 0; j < ndim; ++j) {
1375  tmp[i] += m[j][i] * a[j];
1376  }
1377  c[0] += b[i] * tmp[i];
1378  }
1379  } else {
1380  //B2ERROR("Error, you specified invalid mode= " << mode);
1381  }
1382 
1383  return;
1384 
1385  }
1386 
1387  std::vector<int>
1388  CDCSensitiveDetector::WireId_in_hit_order(int id0, int id1, int nWires)
1389  {
1390  std::vector<int> list;
1391  int i0 = int(id0);
1392  int i1 = int(id1);
1393  if (abs(i0 - i1) * 2 < int(nWires)) {
1394  if (id0 < id1) {
1395  for (int i = id0; i <= id1; ++i)
1396  list.push_back(i);
1397  } else {
1398  for (int i = id0; i >= id1; i--) {
1399  list.push_back(i);
1400  }
1401  }
1402  } else {
1403  if (id0 < id1) {
1404  for (int i = id0; i >= 0; i--)
1405  list.push_back(i);
1406  for (int i = nWires - 1; i >= id1; i--)
1407  list.push_back(i);
1408  } else {
1409  for (int i = id0; i < nWires; ++i)
1410  list.push_back(i);
1411  for (int i = 0; i <= id1; ++i)
1412  list.push_back(i);
1413  }
1414  }
1415 
1416  return list;
1417  }
1418 
1419  G4double CDCSensitiveDetector::ClosestApproach(const G4ThreeVector bwp, const G4ThreeVector fwp, const G4ThreeVector posIn,
1420  const G4ThreeVector posOut, G4ThreeVector& hitPosition, G4ThreeVector& wirePosition)//,G4double& transferT)
1421  {
1422 
1423  TVector3 tbwp(bwp.x(), bwp.y(), bwp.z());
1424  TVector3 tfwp(fwp.x(), fwp.y(), fwp.z());
1425  TVector3 tposIn(posIn.x(), posIn.y(), posIn.z());
1426  TVector3 tposOut(posOut.x(), posOut.y(), posOut.z());
1427  TVector3 thitPosition(0., 0., 0.);
1428  TVector3 twirePosition(0., 0., 0.);
1429 
1430  // G4double distance = m_cdcgp.ClosestApproach(tbwp, tfwp, tposIn, tposOut, thitPosition, twirePosition);
1431  G4double distance = CDC::ClosestApproach(tbwp, tfwp, tposIn, tposOut, thitPosition, twirePosition);
1432 
1433  hitPosition.setX(thitPosition.x());
1434  hitPosition.setY(thitPosition.y());
1435  hitPosition.setZ(thitPosition.z());
1436 
1437  wirePosition.setX(twirePosition.x());
1438  wirePosition.setY(twirePosition.y());
1439  wirePosition.setZ(twirePosition.z());
1440 
1441  return distance;
1442  }
1443 
1444 
1445  //The following-to-end is for setting of left/right flag modified for tracking
1446  void CDCSensitiveDetector::setModifiedLeftRightFlag()
1447  {
1448  if (!m_modifiedLeftRightFlag) return;
1449 
1450  // std::cout <<"#posHits,#negHits= " << m_nPosHits <<" "<< m_nNegHits << std::endl;
1451 
1452  // Get SimHit array and relation betw. MC and SimHit
1453  // N.B. MCParticle is incomplete at this stage; the relation betw it and
1454  // simHit is Okay.
1455  // MCParticle will be completed after all sub-detectors' EndOfEvent calls.
1456  StoreArray<CDCSimHit> simHits;
1457  StoreArray<MCParticle> mcParticles;
1458  RelationArray mcPartToSimHits(mcParticles, simHits);
1459  int nRelationsMinusOne = mcPartToSimHits.getEntries() - 1;
1460 
1461  if (nRelationsMinusOne == -1) return;
1462 
1463  // std::cout <<"#simHits= " << simHits.getEntries() << std::endl;
1464  // std::cout <<"#mcParticles= " << mcParticles.getEntries() << std::endl;
1465  // std::cout <<"#mcPartToSimHits= " << mcPartToSimHits.getEntries() << std::endl;
1466 
1467  //reset some of negative weights to positive; this is needed for the hits
1468  //created by secondary particles whose track-lengths get larger than the
1469  //threshold (set by the user) during G4 swimming (i.e. the weights are
1470  //first set to negative as far as the track-lengths are shorther than the
1471  //threshold; set to positive when the track-lengths exceed the threshold).
1472 
1473  size_t iRelation = 0;
1474  int trackIdOld = INT_MAX;
1475  // std::cout << "INT_MAX= " << INT_MAX << std::endl;
1476  m_hitWithPosWeight.clear();
1477  m_hitWithNegWeight.clear();
1478 
1479  for (int it = nRelationsMinusOne; it >= 0; --it) {
1480  RelationElement& mcPartToSimHit = const_cast<RelationElement&>(mcPartToSimHits[it]);
1481  size_t nRelatedHits = mcPartToSimHit.getSize();
1482  if (nRelatedHits > 1) B2FATAL("CDCSensitiveDetector::EndOfEvent: MCParticle<-> CDCSimHit relation is not one-to-one !");
1483 
1484  unsigned short trackId = mcPartToSimHit.getFromIndex();
1485  RelationElement::weight_type weight = mcPartToSimHit.getWeight(iRelation);
1486  if (weight > 0.) {
1487  trackIdOld = trackId;
1488  } else if (weight <= 0. && trackId == trackIdOld) {
1489  // RelationElement::index_type iSimHit = mcPartToSimHit.getToIndex(iRelation);
1490  weight *= -1.;
1491  mcPartToSimHit.setToIndex(mcPartToSimHit.getToIndex(iRelation), weight);
1492  trackIdOld = trackId;
1493  // std::cout <<"trackId,,iSimHit,wgtafterreset= "<< trackId <<" "<< iSimHit <<" "<< mcPartToSimHit.getWeight(iRelation) << std::endl;
1494  }
1495 
1496  CDCSimHit* sHit = simHits[mcPartToSimHit.getToIndex(iRelation)];
1497 
1498  if (weight > 0.) {
1499  m_hitWithPosWeight.insert(std::pair<unsigned short, CDCSimHit*>(sHit->getWireID().getISuperLayer(), sHit));
1500  } else {
1501  m_hitWithNegWeight.push_back(sHit);
1502  }
1503  }
1504 
1505  /*
1506  // std::cout <<"m_hitWithPosWeight.size= " << m_hitWithPosWeight.size() << std::endl;
1507  for(int i=0; i<9; ++i) {
1508  // if (m_hitWithPosWeight.find(i) != m_hitWithPosWeight.end()) {
1509  // std::cout << i << " found" << std::endl;
1510  // }
1511  m_posWeightMapItBegin.push_back(m_hitWithPosWeight.find(i));
1512  m_posWeightMapItEnd.push_back(m_hitWithPosWeight.find(i+1));
1513  }
1514  */
1515 
1516  //reassign L/R flag
1517  reAssignLeftRightInfo();
1518 
1519  //reset all weights positive; this is required for completing MCParticle object at the EndOfEvent action of FullSim
1520  // is this part really needed ??? check again !
1521  for (int it = 0; it <= nRelationsMinusOne; ++it) {
1522  RelationElement& mcPartToSimHit = const_cast<RelationElement&>(mcPartToSimHits[it]);
1523  RelationElement::weight_type weight = mcPartToSimHit.getWeight(iRelation);
1524  if (weight < 0.) {
1525  mcPartToSimHit.setToIndex(mcPartToSimHit.getToIndex(iRelation), -1.*weight);
1526  }
1527  }
1528 
1529  }
1530 
1531 
1532  void CDCSensitiveDetector::reAssignLeftRightInfo()
1533  {
1534  CDCSimHit* sHit = nullptr;
1535  WireID sWireId; // = WireID();
1536  TVector3 sPos; // = TVector3();
1537 
1538  CDCSimHit* pHit = nullptr;
1539  WireID pWireId; // = WireID();
1540  // double minDistance2 = DBL_MAX;
1541  // double distance2 = DBL_MAX;
1542  // unsigned short bestNeighb = 0;
1543  // unsigned short neighb = 0;
1544 
1545  // std::multimap<unsigned short, CDCSimHit*>::iterator pItBegin = m_hitWithPosWeight.begin();
1546  // std::multimap<unsigned short, CDCSimHit*>::iterator pItEnd = m_hitWithPosWeight.end();
1547 
1548  // unsigned short sClayer = 0;
1549  // unsigned short sSuperLayer = 0;
1550  // unsigned short sLayer = 0;
1551  // unsigned short sWire = 0;
1552  // CDCSimHit* fHit = nullptr;
1553 
1554  //Find a primary track close to the input 2'ndary hit in question
1555  for (std::vector<CDCSimHit*>::iterator nIt = m_hitWithNegWeight.begin(), nItEnd = m_hitWithNegWeight.end(); nIt != nItEnd; ++nIt) {
1556 
1557  sHit = *nIt;
1558  sPos = sHit->getPosTrack();
1559  sWireId = sHit->getWireID();
1560  // sClayer = sWireId.getICLayer();
1561  // sSuperLayer = sWireId.getISuperLayer();
1562  // sLayer = sWireId.getILayer();
1563  // sWire = sWireId.getIWire();
1564  // fHit = sHit;
1565  unsigned short sClayer = sWireId.getICLayer();
1566  unsigned short sSuperLayer = sWireId.getISuperLayer();
1567  unsigned short sLayer = sWireId.getILayer();
1568  unsigned short sWire = sWireId.getIWire();
1569  CDCSimHit* fHit = sHit;
1570 
1571  std::multimap<unsigned short, CDCSimHit*>::iterator pItBegin = m_hitWithPosWeight.find(sSuperLayer);
1572  std::multimap<unsigned short, CDCSimHit*>::iterator pItEnd = m_hitWithPosWeight.find(sSuperLayer + 1);
1573  /*
1574  if (sSuperLayer <= 8) {
1575  pItBegin = m_posWeightMapItBegin.at(sSuperLayer);
1576  pItEnd = m_posWeightMapItEnd.at(sSuperLayer);
1577  } else {
1578  B2FATAL("CDCSensitiveDetector::EndOfEvent: invalid super-layer id ! " << sSuperLayer);
1579  }
1580  */
1581 
1582  double minDistance2 = DBL_MAX;
1583  // bestNeighb = 0;
1584 
1585  /* for (std::multimap<unsigned short, CDCSimHit*>::iterator pIt = m_hitWithPosWeight.begin(); pIt != m_hitWithPosWeight.end(); ++pIt) {
1586  std::cout <<"superLyr#= " << pIt->first << std::endl;
1587  }
1588  */
1589 
1590  for (std::multimap<unsigned short, CDCSimHit*>::iterator pIt = pItBegin; pIt != pItEnd; ++pIt) {
1591 
1592  //scan hits in the same/neighboring cells
1593  pHit = pIt->second;
1594  pWireId = pHit->getWireID();
1595  // neigh = areNeighbors(sWireId, pWireId);
1596  unsigned short neighb = areNeighbors(sClayer, sSuperLayer, sLayer, sWire, pWireId);
1597  if (neighb != 0 || pWireId == sWireId) {
1598  double distance2 = (pHit->getPosTrack() - sPos).Mag2();
1599  if (distance2 < minDistance2) {
1600  fHit = pHit;
1601  minDistance2 = distance2;
1602  // bestNeighb = neighb;
1603  }
1604  }
1605  }
1606 
1607  //reassign LR using the momentum-direction of the primary particle found
1608  unsigned short lR = m_cdcgp->getNewLeftRightRaw(sHit->getPosWire(),
1609  sHit->getPosTrack(),
1610  fHit->getMomentum());
1611  // unsigned short bflr = sHit->getLeftRightPassage();
1612  sHit->setLeftRightPassage(lR);
1613  // std::cout <<"neighb, bfaf lrs, minDistance= " << bestNeighb <<" "<<" "<< bflr <<" "<< sHit->getLeftRightPassage() <<" "<< std::scientific << sqrt(minDistance2) << std::endl;
1614  }
1615  }
1616 
1617 
1618  unsigned short CDCSensitiveDetector::areNeighbors(const WireID& wireId, const WireID& otherWireId) const
1619  {
1620  //require within the same super-layer
1621  if (otherWireId.getISuperLayer() != wireId.getISuperLayer()) return 0;
1622 
1623  const signed short iWire = wireId.getIWire();
1624  const signed short iOtherWire = otherWireId.getIWire();
1625  const signed short iCLayer = wireId.getICLayer();
1626  const signed short iOtherCLayer = otherWireId.getICLayer();
1627 
1628  //require nearby wire
1629  if (iWire == iOtherWire) {
1630  } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1631  } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1632  } else {
1633  return 0;
1634  }
1635  // std::cout <<"iCLayer,iLayer,nShifts= " << iCLayer <<" "<< iLayer <<" "<< nShifts(iCLayer) << std::endl;
1636 
1637  signed short iLayerDifference = otherWireId.getILayer() - wireId.getILayer();
1638  if (abs(iLayerDifference) > 1) return 0;
1639 
1640  if (iLayerDifference == 0) {
1641  if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer))) return CW_NEIGHBOR;
1642  else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) return CCW_NEIGHBOR;
1643  else return 0;
1644  } else if (iLayerDifference == -1) {
1645  // const CCWInfo deltaShift = otherLayer.getShift() - layer.getShift();
1646  const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1647  m_cdcgp->getShiftInSuperLayer(wireId.getISuperLayer(), wireId.getILayer());
1648  // std::cout <<"in deltaShift,iOtherWire,iWire= " << deltaShift <<" "<< iOtherWire <<" "<< iWire << std::endl;
1649  if (iWire == iOtherWire) {
1650  if (deltaShift == CW) return CW_IN_NEIGHBOR;
1651  else if (deltaShift == CCW) return CCW_IN_NEIGHBOR;
1652  else return 0;
1653  } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1654  if (deltaShift == CCW) return CW_IN_NEIGHBOR;
1655  else return 0;
1656  } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1657  if (deltaShift == CW) return CCW_IN_NEIGHBOR;
1658  else return 0;
1659  } else return 0;
1660  } else if (iLayerDifference == 1) {
1661  // const CCWInfo deltaShift = otherLayer.getShift() - layer.getShift();
1662  const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1663  m_cdcgp->getShiftInSuperLayer(wireId.getISuperLayer(), wireId.getILayer());
1664  // std::cout <<"out deltaShift,iOtherWire,iWire= " << deltaShift <<" "<< iOtherWire <<" "<< iWire << std::endl;
1665  if (iWire == iOtherWire) {
1666  if (deltaShift == CW) return CW_OUT_NEIGHBOR;
1667  else if (deltaShift == CCW) return CCW_OUT_NEIGHBOR;
1668  else return 0;
1669  } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1670  if (deltaShift == CCW) return CW_OUT_NEIGHBOR;
1671  else return 0;
1672  } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1673  if (deltaShift == CW) return CCW_OUT_NEIGHBOR;
1674  else return 0;
1675  } else return 0;
1676  } else return 0;
1677 
1678  }
1679 
1680  unsigned short CDCSensitiveDetector::areNeighbors(unsigned short iCLayer, unsigned short iSuperLayer, unsigned short iLayer,
1681  unsigned short iWire, const WireID& otherWireId) const
1682  {
1683  //require within the same super-layer
1684  if (otherWireId.getISuperLayer() != iSuperLayer) return 0;
1685 
1686  const signed short iOtherWire = otherWireId.getIWire();
1687  const signed short iOtherCLayer = otherWireId.getICLayer();
1688 
1689  //require nearby wire
1690  if (iWire == iOtherWire) {
1691  } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1692  } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1693  } else {
1694  return 0;
1695  }
1696 
1697  // std::cout <<"iCLayer,iLayer,nShifts= " << iCLayer <<" "<< iLayer <<" "<< nShifts(iCLayer) << std::endl;
1698  signed short iLayerDifference = otherWireId.getILayer() - iLayer;
1699  if (abs(iLayerDifference) > 1) return 0;
1700 
1701  if (iLayerDifference == 0) {
1702  if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer))) return CW_NEIGHBOR;
1703  else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) return CCW_NEIGHBOR;
1704  else return 0;
1705  } else if (iLayerDifference == -1) {
1706  // const CCWInfo deltaShift = otherLayer.getShift() - layer.getShift();
1707  const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1708  m_cdcgp->getShiftInSuperLayer(iSuperLayer, iLayer);
1709  // std::cout <<"in deltaShift,iOtherWire,iWire= " << deltaShift <<" "<< iOtherWire <<" "<< iWire << std::endl;
1710  if (iWire == iOtherWire) {
1711  if (deltaShift == CW) return CW_IN_NEIGHBOR;
1712  else if (deltaShift == CCW) return CCW_IN_NEIGHBOR;
1713  else return 0;
1714  } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1715  if (deltaShift == CCW) return CW_IN_NEIGHBOR;
1716  else return 0;
1717  } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1718  if (deltaShift == CW) return CCW_IN_NEIGHBOR;
1719  else return 0;
1720  } else return 0;
1721  } else if (iLayerDifference == 1) {
1722  // const CCWInfo deltaShift = otherLayer.getShift() - layer.getShift();
1723  const signed short deltaShift = m_cdcgp->getShiftInSuperLayer(otherWireId.getISuperLayer(), otherWireId.getILayer()) -
1724  m_cdcgp->getShiftInSuperLayer(iSuperLayer, iLayer);
1725  // std::cout <<"out deltaShift,iOtherWire,iWire= " << deltaShift <<" "<< iOtherWire <<" "<< iWire << std::endl;
1726  if (iWire == iOtherWire) {
1727  if (deltaShift == CW) return CW_OUT_NEIGHBOR;
1728  else if (deltaShift == CCW) return CCW_OUT_NEIGHBOR;
1729  else return 0;
1730  } else if (iWire == (iOtherWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iOtherCLayer))) {
1731  if (deltaShift == CCW) return CW_OUT_NEIGHBOR;
1732  else return 0;
1733  } else if ((iWire + 1) % static_cast<signed short>(m_cdcgp->nWiresInLayer(iCLayer)) == iOtherWire) {
1734  if (deltaShift == CW) return CCW_OUT_NEIGHBOR;
1735  else return 0;
1736  } else return 0;
1737  } else return 0;
1738 
1739  }
1740 
1742 } // namespace Belle2
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::WireID
Class to identify a wire inside the CDC.
Definition: WireID.h:44
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:1339
Belle2::CDCSimHit::setLeftRightPassage
void setLeftRightPassage(int minusOneOrZeroOrOne)
The method to set new left/right info. for tracking.
Definition: CDCSimHit.h:181
Belle2::CDC::CDCGeometryPar::wireForwardPosition
const TVector3 wireForwardPosition(int layerId, int cellId, EWirePosition set=c_Base) const
Returns the forward position of the input sense wire.
Definition: CDCGeometryPar.cc:1625
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:150
Belle2::CDC::CDCSensitiveDetector::CW_IN_NEIGHBOR
const signed short CW_IN_NEIGHBOR
Constant for clockwise inwards.
Definition: CDCSensitiveDetector.h:303
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:1388
Belle2::WireID::getILayer
unsigned short getILayer() const
Getter for layer within the Super-Layer.
Definition: WireID.h:146
Belle2::CDC::CDCSensitiveDetector::CW
const signed short CW
Constant for clockwise orientation.
Definition: CDCSensitiveDetector.h:300
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:110
Belle2::CDCSimHit::getPosTrack
TVector3 getPosTrack() const
The method to get position on the track.
Definition: CDCSimHit.h:234
Belle2::CDCSimHit::setPosTrack
void setPosTrack(TVector3 posTrack)
The method to set position on the track.
Definition: CDCSimHit.h:155
Belle2::CDC::CDCSensitiveDetector::m_thresholdKineticEnergy
G4double m_thresholdKineticEnergy
Threshold kinetic energy to be stored.
Definition: CDCSensitiveDetector.h:278
Belle2::CDC::CDCSensitiveDetector::m_magneticField
G4int m_magneticField
Magnetic field is on or off.
Definition: CDCSensitiveDetector.h:257
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:696
Belle2::CDC::CDCSensitiveDetector::CCW
const signed short CCW
Constant for counterclockwise orientation.
Definition: CDCSensitiveDetector.h:299
Belle2::WireID::getISuperLayer
unsigned short getISuperLayer() const
Getter for Super-Layer.
Definition: WireID.h:140
Belle2::CDC::CDCSensitiveDetector::GCUBS
void GCUBS(const G4double x, const G4double y, const G4double d1, const G4double d2, G4double a[4])
Definition: CDCSensitiveDetector.cc:943
Belle2::CDC::CDCSensitiveDetector::reAssignLeftRightInfo
void reAssignLeftRightInfo()
Re-assign left/right info.
Definition: CDCSensitiveDetector.cc:1532
Belle2::CDC::CDCSensitiveDetector::CDCSensitiveDetector
CDCSensitiveDetector(G4String name, G4double thresholdEnergyDeposit, G4double thresholdKineticEnergy)
Constructor.
Definition: CDCSensitiveDetector.cc:53
Belle2::CDC::CDCGeometryPar::getNewLeftRightRaw
unsigned short getNewLeftRightRaw(const TVector3 &posOnWire, const TVector3 &posOnTrack, const TVector3 &momentum) const
Returns new left/right_raw.
Definition: CDCGeometryPar.cc:2652
Belle2::CDCSimHit::setDriftLength
void setDriftLength(double driftLength)
The method to set drift length.
Definition: CDCSimHit.h:108
Belle2::CDCSimHit::setStepLength
void setStepLength(double stepLength)
The method to set step length.
Definition: CDCSimHit.h:120
Belle2::CDC::CDCSensitiveDetector::Rotat
void Rotat(G4double &x, G4double &y, G4double &z, const int mode)
Translation method.
Definition: CDCSensitiveDetector.cc:1014
Belle2::CDCSimHit
Example Detector.
Definition: CDCSimHit.h:33
Belle2::CDC::ClosestApproach
double ClosestApproach(const TVector3 &bwp, const TVector3 &fwp, const TVector3 &posIn, const TVector3 &posOut, TVector3 &hitPosition, TVector3 &wirePosition)
Returns a closest distance between a track and a wire.
Definition: ClosestApproach.cc:27
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:24
Belle2::CDC::CDCSensitiveDetector::step
bool step(G4Step *aStep, G4TouchableHistory *history) override
Process each step and calculate variables defined in CDCB4VHit.
Definition: CDCSensitiveDetector.cc:122
Belle2::CDCSimHit::setWireID
void setWireID(int iCLayerID, int iWireID)
Setter for Wire ID.
Definition: CDCSimHit.h:96
Belle2::CDCSimHit::setPosOut
void setPosOut(TVector3 posOut)
The method to set position of post-step.
Definition: CDCSimHit.h:147
Belle2::CDC::CDCSensitiveDetector::CW_OUT_NEIGHBOR
const signed short CW_OUT_NEIGHBOR
Constant for clockwise outwards.
Definition: CDCSensitiveDetector.h:301
Belle2::CDCSimHit::setEnergyDep
void setEnergyDep(double edep)
The method to set deposited energy.
Definition: CDCSimHit.h:117
Belle2::CDC::CDCSimControlPar::getInstance
static CDCSimControlPar & getInstance()
Static method to get a reference to the CDCSimControlPar instance.
Definition: CDCSimControlPar.cc:18
Belle2::RelationElement
Class to store a single element of a relation.
Definition: RelationElement.h:33
Belle2::CDCSimHit::setPDGCode
void setPDGCode(int pdg)
The method to set PDG code.
Definition: CDCSimHit.h:105
Belle2::CDC::CDCSensitiveDetector::m_hitWithNegWeight
std::vector< CDCSimHit * > m_hitWithNegWeight
Vector containing hits with negative weight.
Definition: CDCSensitiveDetector.h:291
Belle2::CDC::CDCSimControlPar::getMinTrackLength
double getMinTrackLength() const
Get minimum track length.
Definition: CDCSimControlPar.h:135
Belle2::CDC::CDCGeometryPar::wireBackwardPosition
const TVector3 wireBackwardPosition(int layerId, int cellId, EWirePosition set=c_Base) const
Returns the backward position of the input sense wire.
Definition: CDCGeometryPar.cc:1662
Belle2::CDC::CDCGeometryPar::Instance
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Definition: CDCGeometryPar.cc:41
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::CDC::CDCSimControlPar
The Class for CDC Simulation Control Parameters.
Definition: CDCSimControlPar.h:31
Belle2::CDC::CDCSensitiveDetector::m_hitNumber
int m_hitNumber
The current number of created hits in an event.
Definition: CDCSensitiveDetector.h:286
Belle2::CDCSimHit::setLeftRightPassageRaw
void setLeftRightPassageRaw(int minusOneOrZeroOrOne)
The method to set new left/right info. for digitization.
Definition: CDCSimHit.h:173
Belle2::CDC::CDCSensitiveDetector::m_cdcgp
CDCGeometryPar * m_cdcgp
Pointer to CDCGeometryPar object.
Definition: CDCSensitiveDetector.h:268
Belle2::CDC::CDCSensitiveDetector::for_Rotat
void for_Rotat(const G4double bfld[3])
Calculates a rotation matrix.
Definition: CDCSensitiveDetector.cc:978
Belle2::CDC::CDCGeometryPar::getOldLeftRight
unsigned short getOldLeftRight(const TVector3 &posOnWire, const TVector3 &posOnTrack, const TVector3 &momentum) const
Returns old left/right.
Definition: CDCGeometryPar.cc:2617
Belle2::CDCSimHit::setTrackId
void setTrackId(int trackId)
The method to set track id.
Definition: CDCSimHit.h:102
Belle2::CDC::CDCGeometryPar::getShiftInSuperLayer
signed short getShiftInSuperLayer(unsigned short iSuperLayer, unsigned short iLayer) const
Returns shift in the super-layer.
Definition: CDCGeometryPar.cc:2889
Belle2::CDCSimHit::setPosFlag
void setPosFlag(int zeroOrOne)
The method to set position flag.
Definition: CDCSimHit.h:166
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:291
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:589
Belle2::CDC::CDCGeometryPar::cellId
unsigned cellId(unsigned layerId, const TVector3 &position) const
The method to get cell id based on given layer id and the position.
Definition: CDCGeometryPar.cc:1733
Belle2::CDC::CDCSensitiveDetector::CCW_OUT_NEIGHBOR
const signed short CCW_OUT_NEIGHBOR
Constant for counterclockwise outwards.
Definition: CDCSensitiveDetector.h:306
Belle2::CDC::Helix
Helix parameter class.
Definition: Helix.h:51
Belle2::CDC::CDCSensitiveDetector::m_nonUniformField
G4int m_nonUniformField
Magnetic field is uniform or non-uniform.
Definition: CDCSensitiveDetector.h:263
Belle2::CDCSimHit::setMomentum
void setMomentum(TVector3 momentum)
The method to set momentum.
Definition: CDCSimHit.h:123
Belle2::CDC::CDCSensitiveDetector::m_thresholdEnergyDeposit
G4double m_thresholdEnergyDeposit
Threshold energy deposit to be stored.
Definition: CDCSensitiveDetector.h:273
Belle2::RelationElement::getSize
size_t getSize() const
Get number of indices we points to.
Definition: RelationElement.h:81
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:118
Belle2::CDCSimHit::getPosWire
TVector3 getPosWire() const
The method to get position on wire.
Definition: CDCSimHit.h:216
Belle2::CDCSimHit::setPosIn
void setPosIn(TVector3 posIn)
The method to set position of pre-step.
Definition: CDCSimHit.h:139
Belle2::CDC::CDCSensitiveDetector::m_minTrackLength
G4double m_minTrackLength
Min.
Definition: CDCSensitiveDetector.h:284
Belle2::CDC::CDCSensitiveDetector::CCW_NEIGHBOR
const signed short CCW_NEIGHBOR
Constant for counterclockwise.
Definition: CDCSensitiveDetector.h:305
Belle2::RelationElement::getToIndex
index_type getToIndex(size_t n=0) const
Get nth index we point to.
Definition: RelationElement.h:87
Belle2::CDC::CDCSensitiveDetector::setModifiedLeftRightFlag
void setModifiedLeftRightFlag()
set left/right flag modified for tracking
Definition: CDCSensitiveDetector.cc:1446
Belle2::RelationArray::getEntries
int getEntries() const
Get the number of elements.
Definition: RelationArray.h:252
Belle2::CDCSimHit::setGlobalTime
void setGlobalTime(double globalTime)
The method to set global time.
Definition: CDCSimHit.h:114
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:1069
Belle2::CDC::CDCSensitiveDetector::CCW_IN_NEIGHBOR
const signed short CCW_IN_NEIGHBOR
Constant for counterclockwise inwards.
Definition: CDCSensitiveDetector.h:304
Belle2::CDC::CDCGeometryPar::nWiresInLayer
unsigned nWiresInLayer(int layerId) const
Returns wire numbers in a layer.
Definition: CDCGeometryPar.h:1211
Belle2::CDCSimHit::setFlightTime
void setFlightTime(double flightTime)
The method to set flight time.
Definition: CDCSimHit.h:111
Belle2::CDC::CDCSensitiveDetector::Initialize
void Initialize(G4HCofThisEvent *) override
Register CDC hits collection into G4HCofThisEvent.
Definition: CDCSensitiveDetector.cc:96
Belle2::RelationElement::getFromIndex
index_type getFromIndex() const
Get index we point from.
Definition: RelationElement.h:78
Belle2::RelationElement::getWeight
weight_type getWeight(size_t n=0) const
Get nth weight we point to.
Definition: RelationElement.h:90
Belle2::CDC::CDCSimControlPar::getModLeftRightFlag
bool getModLeftRightFlag() const
Get modified left/right flag.
Definition: CDCSimControlPar.h:111
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:1858
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:1618
Belle2::CDCSimHit::getMomentum
TVector3 getMomentum() const
The method to get momentum.
Definition: CDCSimHit.h:210
HepGeom::Point3D< double >
Belle2::StoreArray< MCParticle >
Belle2::Const
This class provides a set of constants for the framework.
Definition: Const.h:36
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:583
Belle2::RelationElement::weight_type
float weight_type
type used for weights.
Definition: RelationElement.h:40
Belle2::WireID::getIWire
unsigned short getIWire() const
Getter for wire within the layer.
Definition: WireID.h:155
Belle2::CDC::CDCSimControlPar::getThresholdEnergyDeposit
double getThresholdEnergyDeposit() const
Get threshold for Energy Deposit;.
Definition: CDCSimControlPar.h:127
Belle2::WireID::getICLayer
unsigned short getICLayer() const
Getter for continuous layer numbering.
Definition: WireID.cc:26
Belle2::CDC::CDCSensitiveDetector::m_hitWithPosWeight
std::multimap< unsigned short, CDCSimHit * > m_hitWithPosWeight
Map containing hits with positive weight.
Definition: CDCSensitiveDetector.h:289
Belle2::CDCSimHit::getWireID
WireID getWireID() const
Getter for WireID object.
Definition: CDCSimHit.h:189
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:102
Belle2::CDCSimHit::setPosWire
void setPosWire(TVector3 posWire)
The method to set position on wire.
Definition: CDCSimHit.h:131
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
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:1419
Belle2::CDC::CDCSimControlPar::getWireSag
bool getWireSag() const
Get wiresag flag.
Definition: CDCSimControlPar.h:103
Belle2::CDC::CDCSensitiveDetector::m_wireSag
G4bool m_wireSag
Switch to activate wire sag effect.
Definition: CDCSensitiveDetector.h:280
Belle2::CDC::CDCSensitiveDetector::CW_NEIGHBOR
const signed short CW_NEIGHBOR
Constant for clockwise.
Definition: CDCSensitiveDetector.h:302
Belle2::CDC::CDCSensitiveDetector::m_modifiedLeftRightFlag
G4bool m_modifiedLeftRightFlag
Switch for left/right flag modified for tracking.
Definition: CDCSensitiveDetector.h:282
Belle2::CDC::CDCSensitiveDetector::m_brot
G4double m_brot[3][3]
a rotation matrix.
Definition: CDCSensitiveDetector.h:265