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