Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4DNAMolecularReaction.cc 64057 2012-10-30 15:04:49Z gcosmo $ 26 // 27 // 27 // Author: Mathieu Karamitros (kara@cenbg.in2p 28 // Author: Mathieu Karamitros (kara@cenbg.in2p3.fr) 28 // 29 // 29 // WARNING : This class is released as a proto 30 // WARNING : This class is released as a prototype. 30 // It might strongly evolve or even disapear i 31 // It might strongly evolve or even disapear in the next releases. 31 // 32 // 32 // History: 33 // History: 33 // ----------- 34 // ----------- 34 // 10 Oct 2011 M.Karamitros created 35 // 10 Oct 2011 M.Karamitros created 35 // 36 // 36 // ------------------------------------------- 37 // ------------------------------------------------------------------- 37 38 38 #include "G4DNAMolecularReaction.hh" 39 #include "G4DNAMolecularReaction.hh" 39 #include "G4DNAMolecularReactionTable.hh" 40 #include "G4DNAMolecularReactionTable.hh" 40 #include "G4VDNAReactionModel.hh" 41 #include "G4VDNAReactionModel.hh" 41 #include "G4MolecularConfiguration.hh" << 42 #include "G4ITManager.hh" 42 #include "G4Molecule.hh" << 43 #include "G4UnitsTable.hh" 43 #include "G4MoleculeFinder.hh" << 44 44 #include "G4ITReactionChange.hh" << 45 using namespace std; 45 #include "G4ITReaction.hh" << 46 46 << 47 G4DNAMolecularReaction::G4DNAMolecularReaction():G4VITReactionProcess(), 47 #include "G4ITTrackHolder.hh" << 48 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)) 48 << 49 G4DNAMolecularReaction::G4DNAMolecularReaction << 50 : << 51 fMolReactionTable(reference_cast<const G4 << 52 , fpReactionModel(nullptr) << 53 { 49 { >> 50 //ctor >> 51 fVerbose = 0; >> 52 fReactionModel = 0; >> 53 fReactionRadius = -1; >> 54 fDistance = -1; 54 } 55 } 55 56 56 G4DNAMolecularReaction::G4DNAMolecularReaction << 57 G4DNAMolecularReaction::~G4DNAMolecularReaction() 57 : G4DNAMolecularReaction() << 58 { 58 { 59 fpReactionModel = pReactionModel; << 59 //dtor >> 60 if(fpChanges) delete fpChanges; 60 } 61 } 61 62 62 G4bool G4DNAMolecularReaction::TestReactibilit << 63 G4DNAMolecularReaction::G4DNAMolecularReaction(const G4DNAMolecularReaction& other):G4VITReactionProcess(other), 63 << 64 fMolReactionTable(reference_cast<const G4DNAMolecularReactionTable*>(fpReactionTable)) 64 << 65 << 66 { 65 { 67 const auto pMoleculeA = GetMolecule(trackA << 66 //copy ctor 68 const auto pMoleculeB = GetMolecule(trackB << 67 fVerbose = other.fVerbose ; 69 << 68 fMolReactionTable = other.fMolReactionTable; 70 const G4double reactionRadius = fpReaction << 69 fReactionModel = 0; >> 70 fReactionRadius = -1; >> 71 fDistance = -1; >> 72 } 71 73 72 G4double separationDistance = -1.; << 74 G4DNAMolecularReaction& G4DNAMolecularReaction::operator=(const G4DNAMolecularReaction& rhs) >> 75 { >> 76 if (this == &rhs) return *this; // handle self assignment 73 77 74 if (currentStepTime == 0.) << 78 fVerbose = rhs.fVerbose ; 75 { << 79 fMolReactionTable = rhs.fMolReactionTable; 76 userStepTimeLimit = false; << 80 fReactionRadius = -1; 77 } << 81 fDistance = -1; 78 82 79 G4bool output = fpReactionModel->FindReact << 83 //assignment operator 80 << 84 return *this; 81 return output; << 82 } 85 } 83 86 84 std::unique_ptr<G4ITReactionChange> G4DNAMolec << 87 G4bool G4DNAMolecularReaction::TestReactibility(const G4Track& trackA, 85 << 88 const G4Track& trackB, >> 89 const double currentStepTime, >> 90 const double /*previousStepTime*/, >> 91 bool userStepTimeLimit) /*const*/ 86 { 92 { 87 std::unique_ptr<G4ITReactionChange> pChang << 93 G4Molecule* moleculeA = GetMolecule(trackA); 88 pChanges->Initialize(trackA, trackB); << 94 G4Molecule* moleculeB = GetMolecule(trackB); 89 95 90 const auto pMoleculeA = GetMolecule(trackA << 96 if(!fReactionModel) 91 const auto pMoleculeB = GetMolecule(trackB << 97 { >> 98 G4ExceptionDescription exceptionDescription ("You have to give a reaction model to the molecular reaction process"); >> 99 G4Exception("G4DNAMolecularReaction::TestReactibility","MolecularReaction001", >> 100 FatalErrorInArgument,exceptionDescription); >> 101 return false; // makes coverity happy >> 102 } >> 103 if(fMolReactionTable==0) >> 104 { >> 105 G4ExceptionDescription exceptionDescription ("You have to give a reaction table to the molecular reaction process"); >> 106 G4Exception("G4DNAMolecularReaction::TestReactibility","MolecularReaction002", >> 107 FatalErrorInArgument,exceptionDescription); >> 108 return false; // makes coverity happy >> 109 } 92 110 93 const auto pReactionData = fMolReactionTab << 111 // Retrieve reaction range >> 112 fReactionRadius = -1 ; // reaction Range >> 113 fReactionRadius = fReactionModel -> GetReactionRadius(moleculeA, moleculeB); 94 114 95 const G4int nbProducts = pReactionData->Ge << 115 fDistance = -1 ; // separation distance 96 116 97 if (nbProducts != 0) << 117 if(currentStepTime == 0.) 98 { 118 { 99 const G4double D1 = pMoleculeA->GetDif << 119 userStepTimeLimit = false; 100 const G4double D2 = pMoleculeB->GetDif << 120 } 101 const G4double sqrD1 = D1 == 0. ? 0. : << 102 const G4double sqrD2 = D2 == 0. ? 0. : << 103 const G4double inv_numerator = 1./(sqr << 104 const G4ThreeVector reactionSite = sqr << 105 + sqr << 106 << 107 for (G4int j = 0; j < nbProducts; ++j) << 108 { << 109 auto pProduct = new G4Molecule(pRe << 110 auto pProductTrack = pProduct->Bui << 111 << 112 pProductTrack->SetTrackStatus(fAli << 113 121 114 G4ITTrackHolder::Instance()->Push( << 122 G4bool output = fReactionModel->FindReaction(trackA, trackB, fReactionRadius,fDistance, userStepTimeLimit); 115 123 116 pChanges->AddSecondary(pProductTra << 124 #ifdef G4VERBOSE 117 G4MoleculeFinder::Instance()->Push << 125 // DEBUG 118 } << 126 if(fVerbose > 1) >> 127 { >> 128 G4cout << "\033[1;39;36m" << "G4MolecularReaction "<< G4endl; >> 129 G4cout << "FindReaction returned : " << G4BestUnit(output,"Length") << G4endl; >> 130 G4cout << " reaction radius : " << G4BestUnit(fReactionRadius,"Length") >> 131 << " real distance : " << G4BestUnit((trackA.GetPosition() - trackB.GetPosition()).mag(), "Length") >> 132 << " calculated distance by model (= -1 if real distance > reaction radius and the user limitation step is not reached) : " >> 133 << G4BestUnit(fDistance,"Length") >> 134 << G4endl; >> 135 >> 136 G4cout << "TrackID A : " << trackA.GetTrackID() >> 137 << ", TrackID B : " << trackB.GetTrackID() >> 138 << " | MolA " << moleculeA->GetName() >> 139 << ", MolB " << moleculeB->GetName() >> 140 <<"\033[0m\n" >> 141 << G4endl; >> 142 G4cout << "--------------------------------------------" << G4endl; 119 } 143 } 120 << 144 #endif 121 pChanges->KillParents(true); << 145 return output ; 122 return pChanges; << 123 } 146 } 124 147 125 void G4DNAMolecularReaction::SetReactionModel( << 126 { << 127 fpReactionModel = pReactionModel; << 128 } << 129 148 130 std::vector<std::unique_ptr<G4ITReactionChange << 149 G4ITReactionChange* G4DNAMolecularReaction::MakeReaction(const G4Track& trackA, const G4Track& trackB) 131 G4ITReactionSet* pReactionSet, << 132 const double currentStepTime, << 133 const double /*fGlobalTime*/, << 134 const bool reachedUserStepTimeLimit) << 135 { 150 { 136 std::vector<std::unique_ptr<G4ITReactionCh << 151 fpChanges = new G4ITReactionChange(); 137 fReactionInfo.clear(); << 152 fpChanges->Initialize(trackA, trackB); 138 153 139 if (pReactionSet == nullptr) << 154 G4Molecule* moleculeA = GetMolecule(trackA); >> 155 G4Molecule* moleculeB = GetMolecule(trackB); >> 156 >> 157 #ifdef G4VERBOSE >> 158 // DEBUG >> 159 if(fVerbose) 140 { 160 { 141 return fReactionInfo; << 161 G4cout << "G4DNAMolecularReaction::MakeReaction" << G4endl; >> 162 G4cout<<"TrackA n°" << trackA.GetTrackID() >> 163 <<"\t | Track B n°" << trackB.GetTrackID() << G4endl; >> 164 >> 165 G4cout<<"Track A : Position : " << G4BestUnit(trackA.GetPosition(),"Length") >> 166 <<"\t Global Time : " << G4BestUnit(trackA.GetGlobalTime(), "Time")<< G4endl; >> 167 >> 168 G4cout<<"Track B : Position : " << G4BestUnit(trackB.GetPosition() ,"Length") >> 169 <<"\t Global Time : " << G4BestUnit(trackB.GetGlobalTime(), "Time")<< G4endl; >> 170 >> 171 G4cout<<"Reaction range : " << G4BestUnit(fReactionRadius,"Length") >> 172 << " \t Separation distance : " << G4BestUnit(fDistance,"Length")<< G4endl; >> 173 G4cout << "--------------------------------------------" << G4endl; 142 } 174 } >> 175 #endif 143 176 144 G4ITReactionPerTrackMap& reactionPerTrackM << 177 const G4DNAMolecularReactionData* reactionData = fMolReactionTable->GetReactionData(moleculeA, moleculeB); 145 for (auto tracks_i = reactionPerTrackMap.c << 146 tracks_i != reactionPerTrackMap.cend( << 147 tracks_i = reactionPerTrackMap.cbegin << 148 { << 149 G4Track* pTrackA = tracks_i->first; << 150 if (pTrackA->GetTrackStatus() == fStop << 151 { << 152 continue; << 153 } << 154 178 155 G4ITReactionPerTrackPtr reactionPerTra << 179 G4int nbProducts = reactionData->GetNbProducts(); 156 G4ITReactionList& reactionList = react << 157 180 158 assert(reactionList.cbegin() != reacti << 181 if (nbProducts) >> 182 { >> 183 G4double D1 = moleculeA->GetDiffusionCoefficient(); >> 184 G4double D2 = moleculeB->GetDiffusionCoefficient(); >> 185 G4double sqrD1 = sqrt(D1); >> 186 G4double sqrD2 = sqrt(D2); >> 187 G4double numerator = sqrD1 + sqrD2; >> 188 G4ThreeVector reactionSite = sqrD1/numerator * trackA.GetPosition() >> 189 + sqrD2/numerator * trackB.GetPosition() ; 159 190 160 for (auto it = reactionList.cbegin(); << 191 for (G4int j=0 ; j < nbProducts; j++) 161 { 192 { 162 G4ITReactionPtr reaction(*it); << 193 G4Molecule* product = new G4Molecule(*reactionData->GetProduct(j)); 163 G4Track* pTrackB = reaction->GetRe << 194 G4Track* productTrack = product->BuildTrack(trackA.GetGlobalTime(), 164 if (pTrackB->GetTrackStatus() == f << 195 reactionSite); 165 { << 196 166 continue; << 197 productTrack->SetTrackStatus(fAlive); 167 } << 198 168 << 199 fpChanges->AddSecondary(productTrack); 169 if (pTrackB == pTrackA) << 200 G4ITManager<G4Molecule>::Instance()->Push(productTrack); 170 { << 171 G4ExceptionDescription excepti << 172 exceptionDescription << 173 << "The IT reaction proces << 174 exceptionDescription << "The p << 175 G4Exception("G4ITModelProcesso << 176 "ITModelProcessor0 << 177 FatalErrorInArgume << 178 exceptionDescripti << 179 } << 180 << 181 pReactionSet->SelectThisReaction(s << 182 << 183 if (TestReactibility(*pTrackA, *pT << 184 { << 185 auto pReactionChange = MakeRea << 186 << 187 if (pReactionChange) << 188 { << 189 fReactionInfo.push_back(st << 190 break; << 191 } << 192 } << 193 } 201 } 194 } 202 } 195 203 196 pReactionSet->CleanAllReaction(); << 204 fpChanges->KillParents(true); 197 return fReactionInfo; << 205 >> 206 return fpChanges; 198 } 207 } 199 208