BlineTracer/src/G4BlineTracer.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file field/BlineTracer/src/G4BlineTracer.cc
0027 /// \brief Implementation of the G4BlineTracer class
0028 //
0029 //
0030 //
0031 //
0032 // --------------------------------------------------------------------
0033 //
0034 // G4BlineTracer implementation
0035 //
0036 // --------------------------------------------------------------------
0037 // Author: Laurent Desorgher (desorgher@phim.unibe.ch)
0038 //         Created - 2003-10-06
0039 // --------------------------------------------------------------------
0040
0041 #include "G4BlineTracer.hh"
0042
0043 #include "G4BlineEquation.hh"
0044 #include "G4BlineEventAction.hh"
0045 #include "G4BlinePrimaryGeneratorAction.hh"
0046 #include "G4BlineSteppingAction.hh"
0047 #include "G4BlineTracerMessenger.hh"
0048 #include "G4CashKarpRKF45.hh"
0049 #include "G4ChordFinder.hh"
0050 #include "G4FieldManager.hh"
0051 #include "G4LogicalVolumeStore.hh"
0052 #include "G4MagIntegratorDriver.hh"
0053 #include "G4PropagatorInField.hh"
0054 #include "G4RunManager.hh"
0055 #include "G4SystemOfUnits.hh"
0056 #include "G4TransportationManager.hh"
0057
0058 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0059
0060 G4BlineTracer::G4BlineTracer()
0061 {
0062   fMessenger = new G4BlineTracerMessenger(this);
0063   fSteppingAction = new G4BlineSteppingAction(this);
0064   fEventAction = new G4BlineEventAction(this);
0065   fPrimaryGeneratorAction = new G4BlinePrimaryGeneratorAction();
0066 }
0067
0068 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0069
0070 G4BlineTracer::~G4BlineTracer()
0071 {
0072   delete fMessenger;
0073   delete fSteppingAction;
0074   delete fEventAction;
0075   delete fPrimaryGeneratorAction;
0076   for (size_t i = 0; i < fVecEquationOfMotion.size(); i++) {
0077     if (fVecEquationOfMotion[i]) delete fVecEquationOfMotion[i];
0078     if (fVecChordFinders[i]) delete fVecChordFinders[i];
0079   }
0080 }
0081
0082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0083
0084 void G4BlineTracer::BeginOfRunAction(const G4Run*) {}
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 void G4BlineTracer::EndOfRunAction(const G4Run*) {}
0089
0090 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0091
0092 void G4BlineTracer::ComputeBlines(G4int n_of_lines)
0093 {
0094   // the first time ResetChordFinders should be called
0095   //
0096   if (!fWas_ResetChordFinders_already_called) {
0097     ResetChordFinders();
0098     fWas_ResetChordFinders_already_called = true;
0099   }
0100
0101   // Replace the user action by the ad-hoc actions for Blines
0102
0103   G4RunManager* theRunManager = G4RunManager::GetRunManager();
0104   auto user_run_action = (G4UserRunAction*)theRunManager->GetUserRunAction();
0105   theRunManager->SetUserAction(this);
0106
0107   auto user_stepping_action = (G4UserSteppingAction*)theRunManager->GetUserSteppingAction();
0108   theRunManager->SetUserAction(fSteppingAction);
0109
0110   auto userPrimaryAction =
0111     (G4VUserPrimaryGeneratorAction*)theRunManager->GetUserPrimaryGeneratorAction();
0112   if (userPrimaryAction) fPrimaryGeneratorAction->SetUserPrimaryAction(userPrimaryAction);
0113   theRunManager->SetUserAction(fPrimaryGeneratorAction);
0114
0115   auto user_event_action = (G4UserEventAction*)theRunManager->GetUserEventAction();
0116   theRunManager->SetUserAction(fEventAction);
0117
0118   auto user_tracking_action = (G4UserTrackingAction*)theRunManager->GetUserTrackingAction();
0119   G4UserTrackingAction* aNullTrackingAction = nullptr;
0120   theRunManager->SetUserAction(aNullTrackingAction);
0121
0122   auto user_stacking_action = (G4UserStackingAction*)theRunManager->GetUserStackingAction();
0123   G4UserStackingAction* aNullStackingAction = nullptr;
0124   theRunManager->SetUserAction(aNullStackingAction);
0125
0126   // replace the user defined chordfinder by the element of fVecChordFinders
0127
0128   std::vector<G4ChordFinder*> user_chord_finders;
0129   std::vector<G4double> user_largest_acceptable_step;
0130   for (size_t i = 0; i < fVecChordFinders.size(); i++) {
0131     user_largest_acceptable_step.push_back(-1.);
0132     if (fVecChordFinders[i]) {
0133       user_chord_finders.push_back(fVecFieldManagers[i]->GetChordFinder());
0134       fVecChordFinders[i]->SetDeltaChord(user_chord_finders[i]->GetDeltaChord());
0135       fVecFieldManagers[i]->SetChordFinder(fVecChordFinders[i]);
0136     }
0137     else
0138       user_chord_finders.push_back(nullptr);
0139   }
0140
0141   // I have tried to use the smooth line filter ability but I could not obtain
0142   // a smooth trajectory in the G4TrajectoryContainer after an event
0143   // Another solution for obtaining a smooth trajectory is to limit
0144   // the LargestAcceptableStep in the G4PropagatorInField object.
0145   // This is the solution I used.
0146
0147   // Old solution:
0148   // G4TransportationManager::GetTransportationManager()
0149   //     ->GetPropagatorInField()->SetTrajectoryFilter(fTrajectoryFilter);
0150
0151   // New solution:
0152   // set the largest_acceptable_step to max_step:length
0153
0154   G4TransportationManager* tmanager = G4TransportationManager::GetTransportationManager();
0155   G4double previous_largest_acceptable_step =
0156     tmanager->GetPropagatorInField()->GetLargestAcceptableStep();
0157
0158   tmanager->GetPropagatorInField()->SetLargestAcceptableStep(fMaxTrackingStep);
0159
0160   // Start the integration of n_of_lines different magnetic field lines
0161
0162   for (G4int il = 0; il < n_of_lines; il++) {
0163     // for each magnetic field line we integrate once backward and once
0164     // forward from the same starting point
0165
0166     // backward integration
0167
0168     for (size_t i = 0; i < fVecEquationOfMotion.size(); i++) {
0169       if (fVecEquationOfMotion[i]) fVecEquationOfMotion[i]->SetBackwardDirectionOfIntegration(true);
0170     }
0171     theRunManager->BeamOn(1);
0172
0173     // forward integration
0174
0175     for (size_t i = 0; i < fVecEquationOfMotion.size(); i++) {
0176       if (fVecEquationOfMotion[i])
0177         fVecEquationOfMotion[i]->SetBackwardDirectionOfIntegration(false);
0178     }
0179     theRunManager->BeamOn(1);
0180   }
0181
0182   // Remove trajectory filter to PropagatorInField
0183   // It was for old solution when using smooth trajectory filter
0184
0185   // tmanager->GetPropagatorInField()->SetTrajectoryFilter(0);
0186
0187   // back to User defined actions and other parameters
0188   // -------------------------------------------------
0189
0190   tmanager->GetPropagatorInField()->SetLargestAcceptableStep(previous_largest_acceptable_step);
0191
0192   // return to User actions
0193
0194   theRunManager->SetUserAction(user_run_action);
0195   theRunManager->SetUserAction(user_event_action);
0196   theRunManager->SetUserAction(userPrimaryAction);
0197   theRunManager->SetUserAction(user_stepping_action);
0198   theRunManager->SetUserAction(user_tracking_action);
0199   theRunManager->SetUserAction(user_stacking_action);
0200
0201   // set user defined chord finders and largest acceptable step
0202
0203   for (size_t i = 0; i < fVecFieldManagers.size(); i++) {
0204     if (user_chord_finders[i]) fVecFieldManagers[i]->SetChordFinder(user_chord_finders[i]);
0205   }
0206 }
0207
0208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0209 ////////////////////////////////////////////////////////////////
0210
0211 /*
0212 G4bool G4BlineTracer::CheckMagneticFields()
0213 {
0214   // Check FieldManagers
0215
0216   G4TransportationManager* tmanager =
0217     G4TransportationManager::GetTransportationManager();
0218
0219   if (fVecFieldManagers[0] != tmanager->GetFieldManager())
0220     return false;
0221   if (fVecMagneticFields[0] != tmanager->GetFieldManager()->GetDetectorField())
0222     return false;
0223   G4LogicalVolumeStore* theVolumeStore = G4LogicalVolumeStore::GetInstance();
0224
0225   std::vector<G4FieldManagers*> LogicalVolumeFields;
0226   size_t j=0;
0227   for (size_t i=0; i<theVolumeStore.size();i++)
0228   {
0229     if (theVolumeStore[i]->GetFieldManager())
0230     {
0231       j++;
0232       if (j >= fVecFieldManagers.size()) return false;
0233       if (fVecFieldManagers[j] != theVolumeStore[i]->GetFieldManager())
0234         return false;
0235       if (fVecMagneticFields[j] !=
0236           theVolumeStore[i]->GetFieldManager()->GetDetectorField())
0237         return false;
0238     }
0239   }
0240   if (j<fVecFieldManagers.size()) return false;
0241
0242  return true;
0243 }
0244 */
0245
0246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0247
0248 void G4BlineTracer::ResetChordFinders()
0249 {
0250   for (size_t i = 0; i < fVecEquationOfMotion.size(); i++) {
0251     delete fVecEquationOfMotion[i];
0252     delete fVecChordFinders[i];
0253   }
0254
0255   fVecChordFinders.clear();
0256   fVecFieldManagers.clear();
0257   fVecMagneticFields.clear();
0258   fVecEquationOfMotion.clear();
0259
0260   // global field
0261
0262   fVecChordFinders.push_back(nullptr);
0263   fVecMagneticFields.push_back(nullptr);
0264   fVecEquationOfMotion.push_back(nullptr);
0265   fVecFieldManagers.push_back(
0266     G4TransportationManager::GetTransportationManager()->GetFieldManager());
0267   if (fVecFieldManagers[0]) {
0268     fVecMagneticFields[0] = (G4MagneticField*)fVecFieldManagers[0]->GetDetectorField();
0269     if (fVecMagneticFields[0]) {
0270       fVecEquationOfMotion[0] = new G4BlineEquation(fVecMagneticFields[0]);
0271       auto pStepper = new G4CashKarpRKF45(fVecEquationOfMotion[0]);
0272       auto pIntgrDriver =
0273         new G4MagInt_Driver(0.01 * mm, pStepper, pStepper->GetNumberOfVariables());
0274       fVecChordFinders[0] = new G4ChordFinder(pIntgrDriver);
0275     }
0276   }
0277
0278   // local fields
0279
0280   G4LogicalVolumeStore* theVolumeStore = G4LogicalVolumeStore::GetInstance();
0281
0282   size_t j = 0;
0283   for (size_t i = 0; i < theVolumeStore->size(); i++) {
0284     if ((*theVolumeStore)[i]->GetFieldManager()) {
0285       j++;
0286       fVecFieldManagers.push_back(((*theVolumeStore)[i])->GetFieldManager());
0287       fVecMagneticFields.push_back((G4MagneticField*)fVecFieldManagers[j]->GetDetectorField());
0288       fVecEquationOfMotion.push_back(nullptr);
0289       fVecChordFinders.push_back(nullptr);
0290       if (fVecMagneticFields[j]) {
0291         fVecEquationOfMotion[j] = new G4BlineEquation(fVecMagneticFields[j]);
0292         auto pStepper = new G4CashKarpRKF45(fVecEquationOfMotion[j]);
0293         auto pIntgrDriver =
0294           new G4MagInt_Driver(.01 * mm, pStepper, pStepper->GetNumberOfVariables());
0295         fVecChordFinders[j] = new G4ChordFinder(pIntgrDriver);
0296       }
0297     }
0298   }
0299 }