repairmodels/src/LEMIVModel.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 LEMIVModel.cc
0028 /// \brief Implementation of the LEMIVModel class
0029
0030 #include "LEMIVModel.hh"
0031
0032 #include "ClassifiedDamage.hh"
0033 #include "Damage.hh"
0034 #include "DamageClassifier.hh"
0035
0036 #include <cmath>
0037 #include <iostream>
0038 #include <fstream>
0039 #include <algorithm>
0040
0041 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0042
0043 LEMIVModel::LEMIVModel(double pLoopLength,double pNi, double pNc,
0044 double pNDSB,double pFunrej,double pTfast,double pTslow):
0045 fLoopLength(pLoopLength),
0046 fNi(pNi),
0047 fNc(pNc),
0048 fNDSB(pNDSB),
0049 fFunrej(pFunrej),
0050 fTfast(pTfast),
0051 fTslow(pTslow)
0052 {
0053     fBpForDSB = 10;
0054     fDefaultsChromosomeSizes={250, 250, 242, 242, 198, 198, 190, 190, 182,
0055                     182, 171, 171, 159, 159, 145, 145, 138, 138, 134,
0056                     134, 135, 135, 133, 133, 114, 114, 107, 107, 102,
0057                     102, 90, 90, 83, 83, 80, 80, 59, 59, 64, 64, 47, 47, 51, 51, 156, 57}; // in Mbp
0058 }
0059
0060 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0061
0062 double LEMIVModel::ComputeUnrej(double pTime)
0063 {
0064     double lambdac = (fNDSB-fNi)/fNc;
0065     double Ffast = fNi/fNDSB;
0066     double Fslow = fNc*lambdac/fNDSB;
0067
0068     return Ffast*std::exp(-std::log(2)*pTime/fTfast)+
0069     (Fslow-fFunrej)*std::exp(-std::log(2)*pTime/fTslow)+
0070     fFunrej;
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 void LEMIVModel::ComputeAndSetDamageInput(std::vector<Damage> vecDamage)
0076 {
0077
0078     double nidsb=0;
0079     double ncdsb=0;
0080     double ndsb=0;
0081
0082     DamageClassifier damclass = DamageClassifier();
0083
0084     auto sortedDamage = damclass.SortDamageByChromo(vecDamage);
0085     // Check chromosome sizes:
0086     if (fChromosomeBpMap.size() == 0) {// using default values
0087         std::cout<<"=====> LEMIV Calculation will Using Default Choromosome sizes"<<std::endl;
0088         for (int ii=0;ii<fDefaultsChromosomeSizes.size();ii++) {
0089             unsigned long long int nBp = fDefaultsChromosomeSizes[ii]*1E6; // convert MBp tp Bp
0090             fChromosomeBpMap.insert({ii,nBp});
0091         }
0092     }
0093     // Loop on each chromosome
0094     int i=0;
0095     for(auto it=sortedDamage.begin();it!=sortedDamage.end();it++)
0096     {
0097         i++;
0098         // Damage are now sorted by event, push all the damage in the same vector
0099         std::vector<Damage> chromoDamage;
0100         for(auto itt=it->second.begin();itt!=it->second.end();itt++)
0101         {
0102             std::move(itt->second.begin(), itt->second.end(), std::back_inserter(chromoDamage));
0103         }
0104
0105         // sort the list of damage by ascending bp
0106         std::sort(chromoDamage.begin(), chromoDamage.end(),
0107             [](const Damage& a, const Damage& b) {
0108                 return a.GetCopyNb() < b.GetCopyNb();
0109             });
0110
0111         // for each loop inside the chromosome
0112         auto chromID = it->first;
0113         if (fChromosomeBpMap.find(chromID) == fChromosomeBpMap.end()) {
0114             std::cerr<<"**** Fatal Error *****"<<std::endl;
0115             std::cerr<<"LEMIVModel::ComputeAndSetDamageInput: Cannot find size info for chrom ID "
0116                     <<chromID<<std::endl;
0117             std::cerr<<"*************** *****"<<std::endl;
0118             exit(EXIT_FAILURE);
0119         }
0120         for(auto startLoop=0;startLoop<fChromosomeBpMap[chromID];startLoop+=fLoopLength)
0121         {
0122
0123             int n = GetDSBPerLoop(chromoDamage,startLoop);
0124
0125             ndsb+=n;
0126
0127             if(n==1)
0128                 nidsb+=1.0;
0129             if(n>=2)
0130                 ncdsb+=1.0;
0131         }
0132     }
0133
0134     // Set in the model input parameters
0135     SetNumDSB(ndsb/fDose);
0136     SetNumDomainIsolated(nidsb/fDose);
0137     SetNumDomainClustered(ncdsb/fDose);
0138 }
0139
0140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0141
0142 int LEMIVModel::GetDSBPerLoop(std::vector<Damage> vecDamage,unsigned int startLoop)
0143 {
0144     // Start to fill a vector with damage having bp between startLopp and startLoop+2Mbp
0145     std::vector<Damage> loopDamage;
0146
0147     for(int i=0;i<vecDamage.size();i++)
0148     {
0149         if((vecDamage[i].GetCopyNb()>startLoop)&&(vecDamage[i].GetCopyNb()<startLoop+fLoopLength))
0150         {
0151             loopDamage.push_back(vecDamage[i]);
0152         }
0153     }
0154
0155     // Make cluster
0156     DamageClassifier dam;
0157     auto classifiedDamage = dam.MakeCluster(loopDamage,fBpForDSB,false);
0158
0159     // Return the number of DSB in this loop
0160     return dam.GetNumDSB(classifiedDamage);
0161 }
0162
0163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0164
0165 void LEMIVModel::CalculateRepair(double pTMax, double pDeltaT)
0166 {
0167     if (pTMax <= 0.) {
0168         std::cout<<"LEMIVModel::CalculateRepair() wrong value for timeMax !!!\n"
0169         <<"Plese check the input macro file!!!"<<std::endl;
0170         exit(0);
0171     }
0172     if (pDeltaT <= 0.) {
0173         std::cout<<"LEMIVModel::CalculateRepair() wrong value for deltaTime  !!!\n"
0174         <<"Plese check the input macro file!!!"<<std::endl;
0175         exit(0);
0176     }
0177     fUCurve.clear();
0178
0179     for(double time=0.;time<=pTMax;time+=pDeltaT)
0180     {
0181         fUCurve.push_back(std::make_pair(time,ComputeUnrej(time)));
0182     }
0183
0184 }
0185
0186 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0187
0188 void LEMIVModel::WriteOutput(std::string pFileName)
0189 {
0190     std::fstream file;
0191     file.open(pFileName.c_str(), std::ios_base::out);
0192     //Header part
0193     file <<"#============================================= LEMIV  MODEL =============================================#\n";
0194     file << "                                  LEMIV Model, CalculateRepair with:\n";
0195     file << "#Number of DSBs: " << fNDSB * fDose    << " DSBs.\n";
0196     file << "#Number of domains with clustered DSB, Nc         = " << fNc * fDose    << " domains.\n";
0197     file << "#Number of domains with isolated DSB, Ni    = " << fNi * fDose<< " domains.\n";
0198     file << "#Funrej     = " << fFunrej << "\n";
0199     file << "#Tfast      = " << fTfast  << " h-1              " << "#Tslow = " << fTslow << " h-1\n";
0200     file << "#LoopLength (length of domain) = " << fLoopLength<<" bp \n";
0201     file <<"#========================================================================================================#\n";
0202     file << "Time (h)\tU\n";
0203     //End Header part
0204     for(int i=0;i<fUCurve.size();i++)
0205     {
0206         file << fUCurve[i].first << "\t" << fUCurve[i].second << "\n";
0207     }
0208
0209     file.close();
0210 }