pdb4dna/src/EventAction.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 EventAction.cc
0027 /// \brief Implementation of the EventAction class
0028
0029 // This example is provided by the Geant4-DNA collaboration
0030 // Any report or published results obtained using the Geant4-DNA software
0031 // shall cite the following Geant4-DNA collaboration publication:
0032 // Med. Phys. 37 (2010) 4692-4708
0033 // Delage et al. PDB4DNA: implementation of DNA geometry from the Protein Data
0034 //                  Bank (PDB) description for Geant4-DNA Monte-Carlo
0035 //                  simulations (submitted to Comput. Phys. Commun.)
0036 // The Geant4-DNA web site is available at http://geant4-dna.org
0037 //
0038 //
0039
0040 #include "EventAction.hh"
0041
0042 #include "EventActionMessenger.hh"
0043
0044 #include "G4AnalysisManager.hh"
0045 #include "G4Event.hh"
0046 #include "G4SystemOfUnits.hh"
0047 #include "G4UnitsTable.hh"
0048 #include "Randomize.hh"
0049
0050 #include <algorithm>
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 EventAction::EventAction() : G4UserEventAction()
0055 {
0056   // default parameter values
0057   //
0058   fThresEdepForSSB = 8.22 * eV;
0059   fThresDistForDSB = 10;
0060   fTotalEnergyDeposit = 0;
0061
0062   // create commands
0063   //
0064   fpEventMessenger = new EventActionMessenger(this);
0065 }
0066
0067 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0068
0069 EventAction::~EventAction()
0070 {
0071   delete fpEventMessenger;
0072 }
0073
0074 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0075
0076 void EventAction::BeginOfEventAction(const G4Event*)
0077 {
0078   // Initialization of parameters
0079   //
0080   fTotalEnergyDeposit = 0.;
0081   fEdepStrand1.clear();
0082   fEdepStrand2.clear();
0083 }
0084
0085 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0086
0087 void EventAction::EndOfEventAction(const G4Event*)
0088 {
0089   // At the end of an event, compute the number of strand breaks
0090   //
0091   G4int sb[2] = {0, 0};
0092   ComputeStrandBreaks(sb);
0093   // Fill histograms
0094   //
0095   G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
0096
0097   if (fTotalEnergyDeposit > 0.) {
0098     analysisManager->FillH1(1, fTotalEnergyDeposit);
0099   }
0100   if (sb[0] > 0) {
0101     analysisManager->FillH1(2, sb[0]);
0102   }
0103   if (sb[1] > 0) {
0104     analysisManager->FillH1(3, sb[1]);
0105   }
0106 }
0107
0108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0109
0110 void EventAction::ComputeStrandBreaks(G4int* sb)
0111 {
0112   // sb quantities
0113   //
0114   G4int ssb1 = 0;
0115   G4int ssb2 = 0;
0116   G4int dsb = 0;
0117
0118   // nucleotide id and energy deposit for each strand
0119   G4int nucl1;
0120   G4int nucl2;
0121   G4double edep1;
0122   G4double edep2;
0123
0124   // Read strand1
0125   //
0126   while (!fEdepStrand1.empty()) {
0127     nucl1 = fEdepStrand1.begin()->first;
0128     edep1 = fEdepStrand1.begin()->second;
0129     fEdepStrand1.erase(fEdepStrand1.begin());
0130
0131     // SSB in strand1
0132     //
0133     if (edep1 >= fThresEdepForSSB / eV) {
0134       ssb1++;
0135     }
0136
0137     // Look at strand2
0138     //
0139     if (!fEdepStrand2.empty()) {
0140       do {
0141         nucl2 = fEdepStrand2.begin()->first;
0142         edep2 = fEdepStrand2.begin()->second;
0143         if (edep2 >= fThresEdepForSSB / eV) {
0144           ssb2++;
0145         }
0146         fEdepStrand2.erase(fEdepStrand2.begin());
0147       } while (((nucl1 - nucl2) > fThresDistForDSB) && (!fEdepStrand2.empty()));
0148
0149       // no dsb
0150       //
0151       if (nucl2 - nucl1 > fThresDistForDSB) {
0152         fEdepStrand2[nucl2] = edep2;
0153         if (edep2 >= fThresEdepForSSB / eV) {
0154           ssb2--;
0155         }
0156       }
0157
0158       // one dsb
0159       //
0160       if (std::abs(nucl2 - nucl1) <= fThresDistForDSB) {
0161         if ((edep2 >= fThresEdepForSSB / eV) && (edep1 >= fThresEdepForSSB / eV)) {
0162           ssb1--;
0163           ssb2--;
0164           dsb++;
0165         }
0166       }
0167     }
0168   }
0169
0170   // End with not processed data
0171   //
0172   while (!fEdepStrand1.empty()) {
0173     nucl1 = fEdepStrand1.begin()->first;
0174     edep1 = fEdepStrand1.begin()->second;
0175     if (edep1 >= fThresEdepForSSB / eV) {
0176       ssb1++;
0177     }
0178     fEdepStrand1.erase(fEdepStrand1.begin());
0179   }
0180
0181   while (!fEdepStrand2.empty()) {
0182     nucl2 = fEdepStrand2.begin()->first;
0183     edep2 = fEdepStrand2.begin()->second;
0184     if (edep2 >= fThresEdepForSSB / eV) {
0185       ssb2++;
0186     }
0187     fEdepStrand2.erase(fEdepStrand2.begin());
0188   }
0189
0190   sb[0] = ssb1 + ssb2;
0191   sb[1] = dsb;
0192 }