Belle II Software development
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
29typedef HepGeom::Point3D<double> HepPoint3D;
30#endif
31#ifndef ENABLE_BACKWARDS_COMPATIBILITY
32typedef HepGeom::Vector3D<double> HepVector3D;
33#endif
34
35
36
37namespace 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 cout(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
963L10:
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 tentative use)
977
978 if (m_nonUniformField == 0) return;
979
980 G4double bx, by, bz;
981 bx = bfld[0];
982 by = bfld[1];
983 bz = bfld[2];
984
985 //cal. rotation matrix
986 G4double bxz, bfield;
987 bxz = bx * bx + bz * bz;
988 bfield = bxz + by * by;
989 bxz = sqrt(bxz);
990 bfield = sqrt(bfield);
991
992 m_brot[0][0] = bz / bxz;
993 m_brot[1][0] = 0.;
994 m_brot[2][0] = -bx / bxz;
995 m_brot[0][1] = -by * bx / bxz / bfield;
996 m_brot[1][1] = bxz / bfield;
997 m_brot[2][1] = -by * bz / bxz / bfield;
998 m_brot[0][2] = bx / bfield;
999 m_brot[1][2] = by / bfield;
1000 m_brot[2][2] = bz / bfield;
1001
1002 return;
1003
1004 }
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 tentative use) of gsim_cdc_rotat/irotat.F in gsim-cdc
1013 //for Belle
1014
1015 if (m_nonUniformField == 0) return;
1016
1017 G4double x0(x), y0(y), z0(z);
1018
1019 if (mode == 1) {
1020 x = m_brot[0][0] * x0 + m_brot[1][0] * y0 + m_brot[2][0] * z0;
1021 y = m_brot[0][1] * x0 + m_brot[1][1] * y0 + m_brot[2][1] * z0;
1022 z = m_brot[0][2] * x0 + m_brot[1][2] * y0 + m_brot[2][2] * z0;
1023 } else if (mode == -1) {
1024 x = m_brot[0][0] * x0 + m_brot[0][1] * y0 + m_brot[0][2] * z0;
1025 y = m_brot[1][0] * x0 + m_brot[1][1] * y0 + m_brot[1][2] * z0;
1026 z = m_brot[2][0] * x0 + m_brot[2][1] * y0 + m_brot[2][2] * z0;
1027 } else {
1028 //B2ERROR("SensitiveDetector " <<"invalid mode " << mode << "specified");
1029 }
1030 return;
1031
1032 }
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 tentative use) of gsim_cdc_rotat/irotat.F in gsim-cdc
1040 //for Belle
1041
1042 if (m_nonUniformField == 0) return;
1043
1044 G4double x0(x[0]), y0(x[1]), z0(x[2]);
1045
1046 if (mode == 1) {
1047 x[0] = m_brot[0][0] * x0 + m_brot[1][0] * y0 + m_brot[2][0] * z0;
1048 x[1] = m_brot[0][1] * x0 + m_brot[1][1] * y0 + m_brot[2][1] * z0;
1049 x[2] = m_brot[0][2] * x0 + m_brot[1][2] * y0 + m_brot[2][2] * z0;
1050 } else if (mode == -1) {
1051 x[0] = m_brot[0][0] * x0 + m_brot[0][1] * y0 + m_brot[0][2] * z0;
1052 x[1] = m_brot[1][0] * x0 + m_brot[1][1] * y0 + m_brot[1][2] * z0;
1053 x[2] = m_brot[2][0] * x0 + m_brot[2][1] * y0 + m_brot[2][2] * z0;
1054 } else {
1055 //B2ERROR("SensitiveDetector " <<"invalid mode " << mode << "specified");
1056 }
1057 return;
1058
1059 }
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 tentative use)
1078 // ---------------------------------------------------------------------
1079 // Purpose : Calculate closest points between helix and wire.
1080 //
1081 // Input
1082 // xwb4 : x of wire at backward endplate in lab.
1083 // ywb4 : y of wire at backward endplate "
1084 // zwb4 : z of wire at backward endplate "
1085 // xwf4 : x of wire at forward endplate "
1086 // ywf4 : y of wire at forward endplate "
1087 // zwf4 : z of wire at forward endplate "
1088 //
1089 // Output
1090 // q2(1) : x of wire at closest point in lab.
1091 // q2(2) : y of wire at closest point "
1092 // q2(3) : z of wire at closest point "
1093 // q1(1) : x of helix at closest point "
1094 // q1(2) : y of helix at closest point "
1095 // q1(3) : z of helix at closest point "
1096 // ntry :
1097 // ---------------------------------------------------------------------
1098
1099 const G4int ndim = 3;
1100 const G4double delta = 1.e-5;
1101
1102
1103 G4double xwb, ywb, zwb, xwf, ywf, zwf;
1104 G4double xw, yw, zw, xh, yh, zh, pxh, pyh, pzh;
1105 G4double fi, fi_corr;
1106
1107 G4double dr, fi0, cpa, dz, tanl;
1108 G4double x0, y0, z0;
1109 // "chrg" removed by M. U. June, 2nd, 2013
1110 // G4double xc, yc, r, chrg;
1111 G4double xc, yc, r;
1112 G4double xwm, ywm;
1113 G4double sinfi0, cosfi0, sinfi0fi, cosfi0fi;
1114
1115 G4double vx, vy, vz, vv, cx, cy, cz, tt[3][3];
1116 G4double tmp[3];
1117
1118 G4double xx[3], dxx[3], ddxx[3], pp[3];
1119 G4double xxtdxx, dxxtdxx, xxtddxx;
1120
1121
1122 G4double fst = 0.0;
1123 G4double f, fderiv, deltafi, fact, eval;
1124 G4double dx1, dy1, dx2, dy2, crs, dot;
1125
1126 G4int iflg;
1127
1128 //set parameters
1129 xwb = xwb4; ywb = ywb4; zwb = zwb4;
1130 xwf = xwf4; ywf = ywf4; zwf = zwf4;
1131
1132 G4double xxx(xp), yyy(yp), zzz(zp);
1133 G4double pxx(px), pyy(py), pzz(pz);
1134
1135 //rotate z-axis to be parallel to B-field in case of non-uniform B
1136 Rotat(xwb, ywb, zwb, 1);
1137 Rotat(xwf, ywf, zwf, 1);
1138 Rotat(xxx, yyy, zzz, 1);
1139 Rotat(pxx, pyy, pzz, 1);
1140
1141 G4double a[8] = {0.};
1142 G4double pt = sqrt(pxx * pxx + pyy * pyy);
1143 a[1] = atan2(-pxx, pyy);
1144 a[2] = charge / pt;
1145 a[4] = pzz / pt;
1146 a[5] = xxx; a[6] = yyy; a[7] = zzz;
1147
1148 //calculate unit direction vector of the sense wire
1149 vx = xwf - xwb; vy = ywf - ywb; vz = zwf - zwb;
1150 vv = sqrt(vx * vx + vy * vy + vz * vz);
1151 vx /= vv; vy /= vv; vz /= vv;
1152
1153 //flag for distinguishing between stereo and axial wire
1154 iflg = 0;
1155 if (vx == 0. && vy == 0.) iflg = 1;
1156 // std::cout << "iflg= " << iflg << std::endl;
1157 //write(6,*) ' hlx2wir ', xwb, ywb, zwb, vx, vy, vz
1158
1159 //calculate coefficients of f
1160 cx = xwb - vx * (vx * xwb + vy * ywb + vz * zwb);
1161 cy = ywb - vy * (vx * xwb + vy * ywb + vz * zwb);
1162 cz = zwb - vz * (vx * xwb + vy * ywb + vz * zwb);
1163
1164 //calculate tensor for f
1165 tt[0][0] = vx * vx - 1.; tt[1][0] = vx * vy; tt[2][0] = vx * vz;
1166 tt[0][1] = vy * vx; tt[1][1] = vy * vy - 1.; tt[2][1] = vy * vz;
1167 tt[0][2] = vz * vx; tt[1][2] = vz * vy; tt[2][2] = vz * vz - 1.;
1168
1169 //set helix parameters
1170 dr = a[0]; fi0 = a[1]; cpa = a[2];
1171 dz = a[3]; tanl = a[4];
1172 x0 = a[5]; y0 = a[6]; z0 = a[7];
1173
1174 //
1175 // set initial value for phi
1176 //
1177
1178 xwm = xxx;
1179 ywm = yyy;
1180 //r(cm) = alpha/cpa = alpha * pt(GeV); bfield(kG)
1181 G4double bfield = sqrt(B_kG[0] * B_kG[0] +
1182 B_kG[1] * B_kG[1] +
1183 B_kG[2] * B_kG[2]);
1184 G4double alpha = 1.e4 / 2.99792458 / bfield;
1185 r = alpha / cpa;
1186 cosfi0 = cos(fi0);
1187 sinfi0 = sin(fi0);
1188 xc = x0 + (dr + r) * cosfi0;
1189 yc = y0 + (dr + r) * sinfi0;
1190 dx1 = x0 - xc;
1191 dy1 = y0 - yc;
1192 dx2 = xwm - xc;
1193 dy2 = ywm - yc;
1194 crs = dx1 * dy2 - dy1 * dx2;
1195 dot = dx1 * dx2 + dy1 * dy2;
1196 fi = atan2(crs, dot);
1197
1198 //begin iterative procedure for newton 's method '
1199 fact = 1.;
1200 ntry = 0;
1201line1:
1202 ntry += 1;
1203 cosfi0fi = cos(fi0 + fi);
1204 sinfi0fi = sin(fi0 + fi);
1205
1206 //calculate spatial point Q(x,y,z) along the helix
1207 xx[0] = x0 + dr * cosfi0 + r * (cosfi0 - cosfi0fi);
1208 xx[1] = y0 + dr * sinfi0 + r * (sinfi0 - sinfi0fi);
1209 xx[2] = z0 + dz - r * tanl * fi;
1210 pp[0] = -pt * sinfi0fi;
1211 pp[1] = pt * cosfi0fi;
1212 pp[2] = pt * tanl;
1213
1214 if (iflg == 1) {
1215 q2[0] = xwb; q2[1] = ywb; q2[2] = xx[2];
1216 q1[0] = xx[0]; q1[1] = xx[1]; q1[2] = xx[2];
1217 q3[0] = pp[0]; q3[1] = pp[1]; q3[2] = pp[2];
1218 //inverse rotation to lab. frame in case of non-uniform B
1219 Rotat(q1, -1);
1220 Rotat(q2, -1);
1221 Rotat(q3, -1);
1222 distance = sqrt((q2[0] - q1[0]) * (q2[0] - q1[0]) +
1223 (q2[1] - q1[1]) * (q2[1] - q1[1]) +
1224 (q2[2] - q1[2]) * (q2[2] - q1[2]));
1225 return;
1226 }
1227
1228 //calculate direction vector (dx/dphi,dy/dphi,dz/dphi)
1229 //on a point along the helix.
1230 dxx[0] = r * sinfi0fi; dxx[1] = - r * cosfi0fi; dxx[2] = - r * tanl;
1231
1232 // In order to derive the closest point between straight line and helix,
1233 // we can put following two conditions:
1234 // (i) A point H(xh,yh,zh) on the helix given should be on
1235 // the plane which is perpendicular to the straight line.
1236 // (ii) A line HW from W(xw,yw,zw) which is a point on the straight
1237 // line to H(xh,yh,zh) should normal to the direction vector
1238 // on the point H.
1239 //
1240 // Thus, we can make a equation from above conditions.
1241 // f(phi) = cx*(dx/dphi) + cy*(dy/dphi) + cz*(dz/dphi)
1242 // + (x,y,z)*tt(i,j)*(dx/dphi,dy/dphi,dz/dphi)
1243 // = 0,
1244 // where
1245 // cx = xwb - vx*( vx*xwb + vy*ywb + vz*zwb )
1246 // cy = ywb - vy*( vx*xwb + vy*ywb + vz*zwb )
1247 // cz = zwb - vz*( vx*xwb + vy*ywb + vz*zwb )
1248 //
1249 // tt(1,1) = vx*vx - 1 tt(1,2) = vx*vy tt(1,3) = vx*vz
1250 // tt(2,1) = vy*vx tt(2,2) = vy*vy - 1 tt(2,3) = vy*vz
1251 // tt(3,1) = vz*vx tt(3,2) = vz*vy tt(3,3) = vz*vz - 1
1252 //
1253 // and the equation of straight line(stereo wire) is written by
1254 // (x,y,z) = (xwb,ywb,zwb) + beta*(vx,vy,vz), beta is free parameter.
1255
1256 //Now calculate f
1257 Mvopr(ndim, xx, tt, dxx, tmp, 1);
1258 xxtdxx = tmp[0];
1259 f = cx * dxx[0] + cy * dxx[1] + cz * dxx[2] + xxtdxx;
1260 if (std::abs(f) < delta) goto line100;
1261
1262 //evaluate fitting result and prepare some factor to multiply to 1/derivative
1263 if (ntry > 1) {
1264 eval = (1.0 - 0.25 * fact) * std::abs(fst) - std::abs(f);
1265 if (eval <= 0.) fact *= 0.5;
1266 }
1267
1268 //calculate derivative of f
1269 ddxx[0] = r * cosfi0fi; ddxx[1] = r * sinfi0fi; ddxx[2] = 0.;
1270
1271 //Now we have derivative of f
1272 Mvopr(ndim, dxx, tt, dxx, tmp, 1);
1273 dxxtdxx = tmp[0];
1274 Mvopr(ndim, xx, tt, ddxx, tmp, 1);
1275 xxtddxx = tmp[0];
1276 fderiv = cx * ddxx[0] + cy * ddxx[1] + cz * ddxx[2] + dxxtdxx + xxtddxx;
1277 // Commented by M. U. June, 2nd, 2013
1278 // fist = fi;
1279 deltafi = f / fderiv;
1280 fi -= fact * deltafi;
1281 fst = f;
1282
1283 if (ntry > ntryMax) {
1284 //B2DEBUG(" Exceed max. trials HelWir ");
1285 goto line100;
1286 }
1287 //write(6,*) ntry, fist, deltafi
1288 goto line1;
1289
1290 //check if zh is btw zwb and zwf; if not, set zh=zwb or zh=zwf.
1291 //dead regions due to feed-throughs should be considered later.
1292line100:
1293 zh = z0 + dz - r * tanl * fi;
1294 fi_corr = 0.;
1295 if (zh < zwb) fi_corr = (zwb - zh) / (-r * tanl);
1296 if (zh > zwf) fi_corr = (zwf - zh) / (-r * tanl);
1297 fi += fi_corr;
1298
1299 cosfi0fi = cos(fi0 + fi);
1300 sinfi0fi = sin(fi0 + fi);
1301
1302 xh = x0 + dr * cosfi0 + r * (cosfi0 - cosfi0fi);
1303 yh = y0 + dr * sinfi0 + r * (sinfi0 - sinfi0fi);
1304 zh = z0 + dz - r * tanl * fi;
1305 pxh = -pt * sinfi0fi;
1306 pyh = pt * cosfi0fi;
1307 pzh = pt * tanl;
1308
1309 //write(6,*) 'fi_corr, zh, zwb, zwf=', fi_corr, zh, zwb, zwf
1310 //write(6,*) 'zh = ', z0, dz, r, tanl, fi
1311
1312 zw = vx * vz * xh + vy * vz * yh + vz * vz * zh + zwb - vz * (vx * xwb + vy * ywb + vz * zwb);
1313 xw = xwb + vx * (zw - zwb) / vz;
1314 yw = ywb + vy * (zw - zwb) / vz;
1315
1316 q2[0] = xw; q2[1] = yw; q2[2] = zw;
1317 q1[0] = xh; q1[1] = yh; q1[2] = zh;
1318 q3[0] = pxh; q3[1] = pyh; q3[2] = pzh;
1319
1320 //inverse rotation to lab. frame in case of non-uniform B
1321 Rotat(q1, -1);
1322 Rotat(q2, -1);
1323 Rotat(q3, -1);
1324 distance = sqrt((q2[0] - q1[0]) * (q2[0] - q1[0]) +
1325 (q2[1] - q1[1]) * (q2[1] - q1[1]) +
1326 (q2[2] - q1[2]) * (q2[2] - q1[2]));
1327 return;
1328
1329 }
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:155
Belle2::CDCSimHit::setPosIn
void setPosIn(const B2Vector3D &posIn)
The method to set position of pre-step.
Definition CDCSimHit.h:121
Belle2::CDCSimHit::setPosFlag
void setPosFlag(int zeroOrOne)
The method to set position flag.
Definition CDCSimHit.h:148
Belle2::CDCSimHit::getPosWire
B2Vector3D getPosWire() const
The method to get position on wire.
Definition CDCSimHit.h:198
Belle2::CDCSimHit::getWireID
WireID getWireID() const
Getter for WireID object.
Definition CDCSimHit.h:171
Belle2::CDCSimHit::setWireID
void setWireID(int iCLayerID, int iWireID)
Setter for Wire ID.
Definition CDCSimHit.h:78
Belle2::CDCSimHit::setGlobalTime
void setGlobalTime(double globalTime)
The method to set global time.
Definition CDCSimHit.h:96
Belle2::CDCSimHit::setTrackId
void setTrackId(int trackId)
The method to set track id.
Definition CDCSimHit.h:84
Belle2::CDCSimHit::setPDGCode
void setPDGCode(int pdg)
The method to set PDG code.
Definition CDCSimHit.h:87
Belle2::CDCSimHit::setMomentum
void setMomentum(const B2Vector3D &momentum)
The method to set momentum.
Definition CDCSimHit.h:105
Belle2::CDCSimHit::setPosWire
void setPosWire(const B2Vector3D &posWire)
The method to set position on wire.
Definition CDCSimHit.h:113
Belle2::CDCSimHit::setEnergyDep
void setEnergyDep(double edep)
The method to set deposited energy.
Definition CDCSimHit.h:99
Belle2::CDCSimHit::setPosOut
void setPosOut(const B2Vector3D &posOut)
The method to set position of post-step.
Definition CDCSimHit.h:129
Belle2::CDCSimHit::setDriftLength
void setDriftLength(double driftLength)
The method to set drift length.
Definition CDCSimHit.h:90
Belle2::CDCSimHit::setStepLength
void setStepLength(double stepLength)
The method to set step length.
Definition CDCSimHit.h:102
Belle2::CDCSimHit::getPosTrack
B2Vector3D getPosTrack() const
The method to get position on the track.
Definition CDCSimHit.h:216
Belle2::CDCSimHit::getMomentum
B2Vector3D getMomentum() const
The method to get momentum.
Definition CDCSimHit.h:192
Belle2::CDCSimHit::setFlightTime
void setFlightTime(double flightTime)
The method to set flight time.
Definition CDCSimHit.h:93
Belle2::CDCSimHit::setPosTrack
void setPosTrack(const B2Vector3D &posTrack)
The method to set position on the track.
Definition CDCSimHit.h:137
Belle2::CDCSimHit::setLeftRightPassage
void setLeftRightPassage(int minusOneOrZeroOrOne)
The method to set new left/right info. for tracking.
Definition CDCSimHit.h:163
Belle2::CDC::CDCGeometryPar::Instance
static CDCGeometryPar & Instance(const CDCGeometry *=nullptr)
Static method to get a reference to the CDCGeometryPar instance.
Definition CDCGeometryPar.cc:31
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::CDC::Helix::ignoreErrorMatrix
void ignoreErrorMatrix(void)
Unsets error matrix.
Definition Helix.cc:872
Belle2::CDC::Helix::pivot
const HepPoint3D & pivot(void) const
returns pivot position.
Definition Helix.h:356
Belle2::CDC::Helix::bFieldZ
double bFieldZ(double bz)
Sets/returns z component of the magnetic field.
Definition Helix.h:473
Belle2::CDC::Helix::a
const HepVector & a(void) const
Returns helix parameters.
Definition Helix.h:428
Belle2::CDC::Helix::momentum
Hep3Vector momentum(double dPhi=0.) const
returns momentum vector after rotating angle dPhi in phi direction.
Definition Helix.cc:259
Belle2::CDC::Helix::x
HepPoint3D x(double dPhi=0.) const
returns position after rotating angle dPhi in phi direction.
Definition Helix.cc:197
Belle2::Const
This class provides a set of constants for the framework.
Definition Const.h:34
Belle2::MCParticleGraph::GraphParticle::setIgnore
void setIgnore(bool ignore=true)
Set or remove the ignore flag.
Definition MCParticleGraph.h:132
Belle2::MCParticleGraph::GraphParticle::getIgnore
bool getIgnore() const
Get the ignore flag.
Definition MCParticleGraph.h:139
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::SensitiveDetectorBase::SensitiveDetectorBase
SensitiveDetectorBase(const std::string &name, Const::EDetector subdetector)
Create a new Sensitive detecor with a given name and belonging to a given subdetector.
Definition SensitiveDetectorBase.h:35
Belle2::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::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition B2Vector3.h:516
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