1 #include <framework/logging/Logger.h>
2 #include <tracking/modules/genfitUtilities/Geant4MaterialInterface.h>
3 #include <geometry/GeometryManager.h>
5 #include "genfit/Exception.h"
6 #include "genfit/RKTrackRep.h"
11 #include "G4ThreeVector.hh"
12 #include "G4Navigator.hh"
13 #include "G4VPhysicalVolume.hh"
14 #include "G4LogicalVolume.hh"
15 #include "G4Material.hh"
16 #include "G4TouchableHistory.hh"
19 #include <G4VExceptionHandler.hh>
20 #include <G4StateManager.hh>
23 static const bool debug =
false;
48 virtual bool Notify(
const char* origin,
const char* code, G4ExceptionSeverity serv,
const char* description)
50 B2WARNING(
"Geant4 exception during material lookup -- origin: "
51 << origin <<
"\n code:"
52 << code <<
"\n description:"
55 if (serv == EventMustBeAborted || serv == RunMustBeAborted) {
59 if (serv == FatalException || serv == FatalErrorInArgument) {
115 nav_.SetWorldVolume(pWorld);
116 worldsolid_ = pWorld->GetLogicalVolume()->GetSolid();
124 const G4ThreeVector* direction = 0,
125 const G4bool pRelativeSearch =
true,
126 const G4bool ignoreDirection =
true);
132 const G4ThreeVector& pDirection,
133 const G4double pCurrentProposedStepLength,
134 G4double& pNewSafety);
141 const G4ThreeVector& direction,
142 const G4TouchableHistory& h);
154 return nav_.CreateTouchableHistory();
165 otherHandler = G4StateManager::GetStateManager()->GetExceptionHandler();
176 G4StateManager::GetStateManager()->SetExceptionHandler(
otherHandler);
180 genfit::Exception exc(
"Geant4MaterialInterface::findNextBoundary ==> Geant4 exception during geometry navigation", __LINE__,
205 nav_.SetGeometricallyLimitedStep();
209 const G4ThreeVector* direction,
210 const G4bool pRelativeSearch,
211 const G4bool ignoreDirection)
219 G4VPhysicalVolume* volume =
nav_.LocateGlobalPointAndSetup(point, direction, pRelativeSearch, ignoreDirection);
223 volume =
nav_.GetWorldVolume();
232 const G4ThreeVector& direction,
233 const G4TouchableHistory& h)
240 G4VPhysicalVolume* volume =
nav_.ResetHierarchyAndLocate(point, direction, h);
244 volume =
nav_.GetWorldVolume();
253 const G4ThreeVector& direction,
254 const G4double pCurrentProposedStepLength,
255 G4double& pNewSafety)
260 return worldsolid_->DistanceToIn(point, direction);
264 const auto distance =
nav_.CheckNextStep(point, direction, pCurrentProposedStepLength, pNewSafety);
271 Geant4MaterialInterface::Geant4MaterialInterface()
275 nav_->SetWorldVolume(world);
278 Geant4MaterialInterface::~Geant4MaterialInterface()
285 double dirX,
double dirY,
double dirZ)
287 G4ThreeVector pos(posX * CLHEP::cm, posY * CLHEP::cm, posZ * CLHEP::cm);
288 G4ThreeVector dir(dirX, dirY, dirZ);
292 nav_->SetGeometricallyLimitedStep();
295 const G4VPhysicalVolume* newVolume =
nav_->LocateGlobalPointAndSetup(pos, &dir);
308 const G4Material* mat =
currentVolume_->GetLogicalVolume()->GetMaterial();
310 double density, Z, A, radiationLength, mEE;
311 if (mat->GetNumberOfElements() == 1) {
317 for (
unsigned i = 0; i < mat->GetNumberOfElements(); ++i) {
318 const G4Element* element = (*mat->GetElementVector())[i];
319 Z += element->GetZ() * mat->GetFractionVector()[i];
320 A += element->GetA() * mat->GetFractionVector()[i];
324 density = mat->GetDensity() / CLHEP::g * CLHEP::cm3;
326 A *= CLHEP::mole / CLHEP::g;
327 radiationLength = mat->GetRadlen() / CLHEP::cm;
328 mEE = mat->GetIonisation()->GetMeanExcitationEnergy() / CLHEP::eV;
342 const double delta(1.E-2);
343 const double epsilon(1.E-1);
347 oldState7 = stateOrig;
349 int stepSign(sMax < 0 ? -1 : 1);
351 G4ThreeVector pointOld(stateOrig[0] * CLHEP::cm,
352 stateOrig[1] * CLHEP::cm,
353 stateOrig[2] * CLHEP::cm);
354 G4ThreeVector dirOld(stepSign * stateOrig[3],
355 stepSign * stateOrig[4],
356 stepSign * stateOrig[5]);
360 const unsigned maxIt = 300;
367 double slDist =
nav_->CheckNextStep(pointOld, dirOld, fabs(sMax) * CLHEP::cm, safety);
368 if (slDist == kInfinity) {
377 if (safety > fabs(sMax)) {
379 std::cout <<
" next boundary is farther away than sMax \n";
380 return stepSign * safety;
384 if (slDist < delta) {
386 std::cout <<
" very close to the boundary -> return @ it " << it
387 <<
" stepSign*slDist = "
388 << stepSign <<
"*" << slDist <<
"\n";
389 return stepSign * slDist;
391 double step = slDist;
395 genfit::Exception exc(
"Geant4MaterialInterface::findNextBoundary ==> maximum number of iterations exceeded", __LINE__, __FILE__);
401 return stepSign * (s + step);
409 rep->
RKPropagate(state7, NULL, SA, stepSign * (s + step), varField);
411 G4ThreeVector pos(state7[0] * CLHEP::cm, state7[1] * CLHEP::cm, state7[2] * CLHEP::cm);
412 G4ThreeVector dir(stepSign * state7[3], stepSign * state7[4], stepSign * state7[5]);
416 double dist2 = (pow(state7[0] - oldState7[0], 2)
417 + pow(state7[1] - oldState7[1], 2)
418 + pow(state7[2] - oldState7[2], 2));
422 if (dist2 > safety * safety) {
431 double maxDeviation2 = 0.25 * (step * step - dist2);
432 if (maxDeviation2 > epsilon * epsilon) {
437 step = safety + 0.5 * (step - safety);
448 nav_->SetGeometricallyLimitedStep();
451 std::unique_ptr<G4TouchableHistory> hist(
nav_->CreateTouchableHistory());
452 G4VPhysicalVolume* newVolume =
nav_->LocateGlobalPointAndSetup(pos, &dir);
463 return stepSign * (s + step);
467 G4ThreeVector dirCloser(pos - pointOld);
468 dirCloser.setMag(1.);
469 nav_->ResetHierarchyAndLocate(pointOld, dirCloser, *hist);
474 double secantDist =
nav_->CheckNextStep(pointOld, dirCloser,
475 step * CLHEP::cm, safety) / CLHEP::cm;
477 if (secantDist >= step) {
481 step = std::max(0.9 * step, safety);
483 slDist = step = std::max(secantDist, safety);
487 if (slDist < delta) {
489 std::cout <<
" very close to the boundary -> return @ it " << it
490 <<
" stepSign*(s + slDist) = "
491 << stepSign <<
"*(" << s + slDist <<
")\n";
492 return stepSign * (s + slDist);
508 nav_->SetGeometricallyLimitedStep();
510 nav_->LocateGlobalPointAndSetup(pos, &dir);
511 slDist =
nav_->CheckNextStep(pos, dir,
512 (fabs(sMax) - s) * CLHEP::cm, safety);
513 if (slDist == kInfinity) {
515 slDist = fabs(sMax) - s;
523 if (s + safety > fabs(sMax)) {
525 std::cout <<
" next boundary is farther away than sMax \n";
526 return stepSign * (s + safety);
530 if (slDist < delta) {
532 std::cout <<
" very close to the boundary -> return @ it " << it
533 <<
" stepSign*(s + slDist) = "
534 << stepSign <<
"*(" << s + slDist <<
")\n";
535 return stepSign * (s + slDist);