| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// \file TimeStepAction.cc
28 /// \brief Implementation of the TimeStepAction class
29
30 #include "TimeStepAction.hh"
31 #include "G4UnitsTable.hh"
32 #include "G4SystemOfUnits.hh"
33
34 #include "G4DNAChemistryManager.hh"
35 #include "G4MoleculeCounter.hh"
36 #include "G4MoleculeTable.hh"
37 #include "G4Scheduler.hh"
38 #include "UserMolecule.hh"
39 #include "G4ChemTimeStepModel.hh"
40 #include "G4EmParameters.hh"
41 #include "Analysis.hh"
42 #include "G4DNAMolecule.hh"
43 #include "G4H3O.hh"
44 #include "G4OH.hh"
45 #include "G4H2O2.hh"
46 #include "G4H2.hh"
47 #include "G4Hydrogen.hh"
48 #include "G4Electron_aq.hh"
49 #include "G4ParticleTable.hh"
50 #include "G4ParticleDefinition.hh"
51 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
52
53 TimeStepAction::TimeStepAction()
54 : G4UserTimeStepAction()
55 {
56 auto param = G4EmParameters::Instance();
57 if (param->GetTimeStepModel() == G4ChemTimeStepModel::SBS) {
58 AddTimeStep(1*picosecond,0.1*picosecond);
59 AddTimeStep(10*picosecond,1*picosecond);
60 AddTimeStep(100*picosecond,3*picosecond);
61 AddTimeStep(1000*picosecond,10*picosecond);
62 AddTimeStep(10000*picosecond,100*picosecond);
63 }
64 }
65
66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
67
68
69 TimeStepAction::TimeStepAction(const TimeStepAction& other)
70 : G4UserTimeStepAction(other)
71 {}
72
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
74
75 TimeStepAction& TimeStepAction::operator=(const TimeStepAction& rhs)
76 {
77 if (this == &rhs) return *this; // handle self assignment
78 //assignment operator
79 return *this;
80 }
81
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83
84 void TimeStepAction::StartProcessing()
85 {
86
87 G4cout<<"Start chemistry "<<G4Scheduler::Instance()->GetGlobalTime()/s<<" s"<<G4endl;
88 G4cout<<"Chemical endTime "<<G4Scheduler::Instance()->GetEndTime()/nanosecond<<" ns"<<G4endl;
89 }
90
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
92
93 void TimeStepAction::UserReactionAction(const G4Track& a, const G4Track& b,
94 const std::vector<G4Track*>* products)
95 {
96
97 // ***************************************
98 // Flag Part
99 // ***************************************
100 //
101 // set the flags
102
103 fReactif1 = 0;
104 fReactif2 = 0;
105 fProduct1 = 0;
106 fProduct2 = 0;
107
108 fReactif1 = SetParticleFlag(a.GetDynamicParticle()->GetDefinition());
109 fReactif2 = SetParticleFlag(b.GetDynamicParticle()->GetDefinition());
110 if(products)
111 {
112 for(unsigned int i=0;i<products->size();++i)
113 {
114 switch(i)
115 {
116 case 0:
117 fProduct1 = SetParticleFlag(((*products)[i])->GetDynamicParticle()->GetDefinition());
118 break;
119 case 1:
120 fProduct2 = SetParticleFlag(((*products)[i])->GetDynamicParticle()->GetDefinition());;
121 break;
122 default:
123 break;
124 }
125 }
126 }
127
128 // ***************************************
129 // Save Part
130 // ***************************************
131 // For DBScan
132 G4int base;
133
134 if( (fReactif1==3 || fReactif2 ==3) // OH
135 && (
136 (fReactif1==8 || fReactif2 ==8) // Desoxyribose
137 || (fReactif1==9 || fReactif2 ==9) // Phosphate
138 || (fReactif1==10|| fReactif2 ==10) // Adenine
139 || (fReactif1==11|| fReactif2 ==11) // Thymine
140 || (fReactif1==12|| fReactif2 ==12) // Guanine
141 || (fReactif1==13|| fReactif2 ==13) // Cytosine
142 ) )
143 {
144 // Retrieve the dna molecule
145 const G4Track* dnaMolecule = nullptr;
146 const G4Track* radical = nullptr;
147 //
148 if (a.GetDynamicParticle()->GetDefinition() == G4OH::Definition()) {
149 dnaMolecule = &b;
150 radical = &a;
151 }
152 else {
153 dnaMolecule = &a;
154 radical = &b;
155 }
156
157 auto _dnaMolecule = dynamic_cast<UserMolecule*> (dnaMolecule->GetUserInformation());
158 if (_dnaMolecule) {
159 G4int copyNo = _dnaMolecule->GetCopyNumber();
160 G4int strand = _dnaMolecule->GetStrand();
161
162 if((fReactif1==8 || fReactif2 ==8) // Desoxyribose
163 || (fReactif1==9 || fReactif2 ==9)) // Phosphate
164 {
165 base=0;
166 }
167 else
168 {
169 base=1;
170 }
171
172 // Retrieve the molecule coordinates
173 auto analysisManager = Analysis::GetAnalysis()->GetAnalysisManager();
174 analysisManager->FillNtupleIColumn(1, 0, strand);
175 analysisManager->FillNtupleIColumn(1, 1, copyNo);
176 analysisManager->FillNtupleDColumn(1, 2, radical->GetPosition().getX()/nm);
177 analysisManager->FillNtupleDColumn(1, 3, radical->GetPosition().getY()/nm);
178 analysisManager->FillNtupleDColumn(1, 4, radical->GetPosition().getZ()/nm);
179 analysisManager->FillNtupleDColumn(1, 5, G4Scheduler::Instance()->GetGlobalTime()/ns );
180 analysisManager->FillNtupleIColumn(1, 6, base);
181 analysisManager->AddNtupleRow(1);
182 } else {
183 G4String msg = "Error in Downcasting";
184 G4Exception("TimeStepAction::UserReactionAction", "", FatalException, msg);
185 }
186
187 }
188
189 }
190
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
192
193 G4int TimeStepAction::SetParticleFlag(const G4ParticleDefinition* partDef)
194 {
195 G4int flag (0);
196
197 if(partDef == G4H3O::Definition()) flag = 1;
198 else if(partDef == G4OH::Definition()) flag = 3;
199 else if(partDef == G4Electron_aq::Definition()) flag = 4;
200 else if(partDef == G4Hydrogen::Definition()) flag = 5;
201 else if(partDef == G4H2::Definition()) flag = 6;
202 else if(partDef == G4H2O2::Definition()) flag = 7;
203 else if(partDef == G4Deoxyribose::Definition()) flag = 8;
204 else if(partDef == G4Phosphate::Definition()) flag = 9;
205 else if(partDef == G4Adenine::Definition()) flag = 10;
206 else if(partDef == G4Thymine::Definition()) flag = 11;
207 else if(partDef == G4Guanine::Definition()) flag = 12;
208 else if(partDef == G4Cytosine::Definition()) flag = 13;
209 else if(partDef == G4Histone::Definition()) flag = 14;
210 else if(partDef == G4DamagedDeoxyribose::Definition()) flag = 15;
211 else if(partDef == G4DamagedAdenine::Definition()) flag = 16;
212 else if(partDef == G4DamagedThymine::Definition()) flag = 17;
213 else if(partDef == G4DamagedCytosine::Definition()) flag = 18;
214 else if(partDef == G4DamagedGuanine::Definition()) flag = 19;
215 else {
216 G4ParticleDefinition* OHm =G4ParticleTable::GetParticleTable()->FindParticle("OHm");
217 if (OHm && partDef == OHm) {
218 flag = 2;
219 }
220 }
221 return flag;
222 }
223
224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
225