dsbandrepair/src/ChemTimeStepAction.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 //
0027 /// \file ChemTimeStepAction.cc
0028 /// \brief Implementation of the ChemTimeStepAction class
0029
0030 #include "ChemTimeStepAction.hh"
0031 #include "G4UnitsTable.hh"
0032 #include "G4SystemOfUnits.hh"
0033
0034 #include "G4DNAChemistryManager.hh"
0035 #include "G4MoleculeCounter.hh"
0036 #include "G4MoleculeTable.hh"
0037 #include "G4Scheduler.hh"
0038 #include "UserMolecule.hh"
0039 #include "G4ChemTimeStepModel.hh"
0040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0041
0042 ChemTimeStepAction::ChemTimeStepAction(ChemNtupleManager* histo, G4ChemTimeStepModel timeStepModel)
0043 : G4UserTimeStepAction(),fpHisto(histo)
0044 {
0045     if (timeStepModel == G4ChemTimeStepModel::SBS) {
0046         AddTimeStep(1*picosecond,0.1*picosecond);
0047         AddTimeStep(10*picosecond,1*picosecond);
0048         AddTimeStep(100*picosecond,3*picosecond);
0049         AddTimeStep(1000*picosecond,10*picosecond);
0050         AddTimeStep(10000*picosecond,100*picosecond);
0051     }
0052 }
0053
0054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0055
0056
0057 ChemTimeStepAction::ChemTimeStepAction(const ChemTimeStepAction& other)
0058 : G4UserTimeStepAction(other)
0059 {}
0060
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0062
0063 ChemTimeStepAction& ChemTimeStepAction::operator=(const ChemTimeStepAction& rhs)
0064 {
0065     if (this == &rhs) return *this; // handle self assignment
0066     //assignment operator
0067     return *this;
0068 }
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071
0072 void ChemTimeStepAction::StartProcessing()
0073 {
0074
0075     G4cout<<"Start chemistry "<<G4Scheduler::Instance()->GetGlobalTime()/s<<" s"<<G4endl;
0076     G4cout<<"Chemical endTime "<<G4Scheduler::Instance()->GetEndTime()/nanosecond<<" ns"<<G4endl;
0077 }
0078
0079 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0080
0081 void ChemTimeStepAction::UserReactionAction(const G4Track& a, const G4Track& b,
0082     const std::vector<G4Track*>* products)
0083 {
0084
0085     // ***************************************
0086     // Flag Part
0087     // ***************************************
0088     //
0089     // set the flags
0090
0091     fReactif1 = 0.;
0092     fReactif2 = 0.;
0093     fProduct1 = 0.;
0094     fProduct2 = 0.;
0095
0096     G4String aName = a.GetDynamicParticle()->GetDefinition()->GetParticleName();
0097     fReactif1 = SetFlag(aName);
0098
0099     G4String bName = b.GetDynamicParticle()->GetDefinition()->GetParticleName();
0100     fReactif2 = SetFlag(bName);
0101
0102     if(products)
0103     {
0104         for(unsigned int i=0;i<products->size();++i)
0105         {
0106             switch(i)
0107             {
0108                 case 0:
0109                 fProduct1 = SetFlag(((*products)[i])->GetDynamicParticle()->GetDefinition()->GetParticleName());
0110                 break;
0111                 case 1:
0112                 fProduct2 = SetFlag(((*products)[i])->GetDynamicParticle()->GetDefinition()->GetParticleName());;
0113                 break;
0114                 default:
0115                 break;
0116             }
0117         }
0118     }
0119
0120     // ***************************************
0121     // Save Part
0122     // ***************************************
0123     // For DBScan
0124     G4int base;
0125
0126     if( (fReactif1==3 || fReactif2 ==3) // OH
0127         && (
0128             (fReactif1==8 || fReactif2 ==8) // Desoxyribose
0129                 || (fReactif1==9 || fReactif2 ==9) // Phosphate
0130         || (fReactif1==10|| fReactif2 ==10) // Adenine
0131             || (fReactif1==11|| fReactif2 ==11) // Thymine
0132             || (fReactif1==12|| fReactif2 ==12) // Guanine
0133             || (fReactif1==13|| fReactif2 ==13) // Cytosine
0134                 ) )
0135     {
0136         // Retrieve the dna molecule
0137         const G4Track* dnaMolecule;
0138         const G4Track* radical;
0139         //
0140         if (aName == "OH") {
0141             dnaMolecule = &b;
0142             radical = &a;
0143         }
0144         else {
0145             dnaMolecule = &a;
0146             radical = &b;
0147         }
0148
0149         auto _dnaMolecule = dynamic_cast<UserMolecule*> (dnaMolecule->GetUserInformation());
0150         if (_dnaMolecule) {
0151             G4int copyNo = _dnaMolecule->GetCopyNumber();
0152             G4int strand = _dnaMolecule->GetStrand();
0153
0154             if((fReactif1==8 || fReactif2 ==8) // Desoxyribose
0155                 || (fReactif1==9 || fReactif2 ==9)) // Phosphate
0156                 {
0157                     base=0;
0158                 }
0159             else
0160                 {
0161                     base=1;
0162                 }
0163
0164             // Retrieve the molecule coordinates
0165             fpHisto->FillNtupleDColumn(1, 0, strand);
0166             fpHisto->FillNtupleDColumn(1, 1, copyNo);
0167             fpHisto->FillNtupleDColumn(1, 2, radical->GetPosition().getX()/nm);
0168             fpHisto->FillNtupleDColumn(1, 3, radical->GetPosition().getY()/nm);
0169             fpHisto->FillNtupleDColumn(1, 4, radical->GetPosition().getZ()/nm);
0170             fpHisto->FillNtupleDColumn(1, 5, G4Scheduler::Instance()->GetGlobalTime()/ns );
0171             fpHisto->FillNtupleDColumn(1, 6, base);
0172             fpHisto->AddNtupleRow(1);
0173         } else {
0174             G4String msg = "Error in Downcasting";
0175             G4Exception("ChemTimeStepAction::UserReactionAction", "", FatalException, msg);
0176         }
0177
0178     }
0179
0180 }
0181
0182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0183
0184 G4double ChemTimeStepAction::SetFlag(const G4String& val)
0185 {
0186     G4double flag (0);
0187
0188     if(val == "H3O") flag = 1;
0189     else if(val == "OHm") flag = 2;
0190     else if(val == "OH") flag = 3;
0191     else if(val == "e_aq") flag = 4;
0192     else if(val == "H") flag = 5;
0193     else if(val == "H_2") flag = 6;
0194     else if(val == "H2O2") flag = 7;
0195     else if(val == "Deoxyribose") flag = 8;
0196     else if(val == "Phosphate") flag = 9;
0197     else if(val == "Adenine") flag = 10;
0198     else if(val == "Thymine") flag = 11;
0199     else if(val == "Guanine") flag = 12;
0200     else if(val == "Cytosine") flag = 13;
0201     else if(val == "Histone") flag = 14;
0202     else if(val == "Damaged_Deoxyribose") flag = 15;
0203     else if(val == "Damaged_Adenine") flag = 16;
0204     else if(val == "Damaged_Thymine") flag = 17;
0205     else if(val == "Damaged_Cytosine") flag = 18;
0206     else if(val == "Damaged_Guanine") flag = 19;
0207
0208     return flag;
0209 }
0210
0211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......