analysis/src/AnalysisHandler.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 AnalysisHandler.cc
0028 /// \brief Implementation of the AnalysisHandler class
0029
0030 #include "AnalysisHandler.hh"
0031 #include "ScanDamage.hh"
0032 #include "DamageClassifier.hh"
0033 #include "ParametersParser.hh"
0034 #include "SDDData.hh"
0035
0036 #include <cmath>
0037 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0038
0039 AnalysisHandler::AnalysisHandler(): pTLKDoseMax(6.0), pTLKDeltaDose(0.25), fNBp (-1)
0040 {
0041     fScanDamage = std::make_unique<ScanDamage>();
0042     fTLKModel = std::make_unique<TLKModel>();
0043     fLEMIVModel = std::make_unique<LEMIVModel>();
0044     fBelovModel = std::make_unique<BelovModel>();
0045     fBpForDSB = 10;
0046
0047     if (ParametersParser::Instance()->GetThresholdE() != "") {
0048         auto e = std::stod(ParametersParser::Instance()->GetThresholdE());
0049         if (e > 0) fScanDamage->SetThresholdEnergy(e);
0050     }
0051     if (ParametersParser::Instance()->GetProbabilityForIndirectSB() != "") {
0052         auto p = std::stod(ParametersParser::Instance()->GetProbabilityForIndirectSB());
0053         if (p > 0 && p <= 100) fScanDamage->SetProbabilityForIndirectSBSelection(p/100.);
0054     }
0055 }
0056
0057 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0058
0059 void AnalysisHandler::SetThresholdEnergy(double e)
0060 {
0061     fScanDamage->SetThresholdEnergy(e);
0062 }
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0065
0066 void AnalysisHandler::GetAllDamageAndScanSB()
0067 {
0068     std::map<unsigned int,std::map<unsigned int,std::vector<Damage> > > dmMap;
0069     if (ParametersParser::Instance()->WannaLoadDamagesFromSDD()) {
0070         SDDData sdddata(ParametersParser::Instance()->GetSDDFileName());
0071         dmMap = sdddata.GetAllDamage();
0072         fDose = sdddata.GetDose();
0073         fChromosomeBpMap = sdddata.GetChromosomeBpSizesMap(fNBp);
0074     } else {
0075         if (ParametersParser::Instance()->WannaSkipScanningIndirectDamage()) {
0076             fScanDamage->SkipScanningIndirectDamage();
0077         }
0078         dmMap = fScanDamage->ExtractDamage();
0079         fEdepInNucleus = fScanDamage->GetEdepSumInNucleus();//eV
0080         double nuclesumass = fScanDamage->GetNucleusMass(); // kg
0081         double eVtoJ = 1.60E-19;
0082         fDose = fEdepInNucleus*eVtoJ/nuclesumass;
0083         fNBp = fScanDamage->GetTotalNbBpPlacedInGeo();
0084         fChromosomeBpMap = fScanDamage->GetChromosomeBpSizesMap();
0085     }
0086     DamageClassifier damClass;
0087     std::map<int,int> ndsbMap, ncdsbMap, nssbMap, nsbMap;
0088     std::map<int,int> ndirsbMap, ndsbdirMap, ndsbdirIMap, ndsbInMap;
0089     for (const auto& [chromo,evtDm] : dmMap) {
0090         for (const auto& [evt, dmV] : evtDm) {
0091             fAllDamage.insert(fAllDamage.end(),dmV.begin(),dmV.end()); // to write SDD file and for LEM-IV
0092             std::vector<Damage> tmpV{dmV};
0093             auto classifiedDamage = damClass.MakeCluster(tmpV,fBpForDSB,false);
0094
0095             for (auto dm : dmV) {
0096                 if (dm.GetDamageType() == Damage::Damage::fBackbone ) {
0097                     if ( (nsbMap.find(evt) == nsbMap.end()) ) {
0098                         nsbMap.insert({evt,1});
0099                     } else {
0100                         nsbMap[evt] ++;
0101                     }
0102                     if (dm.GetCause() == Damage::Damage::fDirect) {
0103                         if ( (ndirsbMap.find(evt) == ndirsbMap.end()) ) {
0104                             ndirsbMap.insert({evt,1});
0105                         } else {
0106                             ndirsbMap[evt] ++;
0107                         }
0108                     }
0109                 }
0110             }
0111
0112             int NumDSBForThisCluster = damClass.GetNumDSB(classifiedDamage);
0113             if ( (ndsbMap.find(evt) == ndsbMap.end()) ) {
0114                 if (NumDSBForThisCluster > 0) ndsbMap.insert({evt,NumDSBForThisCluster});
0115             } else {
0116                 ndsbMap[evt] += NumDSBForThisCluster;
0117             }
0118
0119             int NumcDSBForThisCluster = damClass.GetNumComplexDSB(classifiedDamage);
0120             if ( (ncdsbMap.find(evt) == ncdsbMap.end()) ) {
0121                 if (NumcDSBForThisCluster > 0) ncdsbMap.insert({evt,NumcDSBForThisCluster});
0122             } else {
0123                 ncdsbMap[evt] += NumcDSBForThisCluster;
0124             }
0125
0126             int NumSSBForThisCluster = damClass.GetNumSSB(classifiedDamage);
0127             if ( (nssbMap.find(evt) == nssbMap.end()) ) {
0128                 if (NumSSBForThisCluster > 0) nssbMap.insert({evt,NumSSBForThisCluster});
0129             } else {
0130                 nssbMap[evt] += NumSSBForThisCluster;
0131             }
0132
0133             int NumDSBdirForThisCluster = damClass.GetNumDSBwithDirectDamage(classifiedDamage);
0134             if ( (ndsbdirMap.find(evt) == ndsbdirMap.end()) ) {
0135                 if (NumDSBdirForThisCluster > 0) ndsbdirMap.insert({evt,NumDSBdirForThisCluster});
0136             } else {
0137                 ndsbdirMap[evt] += NumDSBdirForThisCluster;
0138             }
0139
0140             int NumDSBInForThisCluster = damClass.GetNumDSBwithIndirectDamage(classifiedDamage);
0141             if ( (ndsbInMap.find(evt) == ndsbInMap.end()) ) {
0142                 if (NumDSBInForThisCluster > 0) ndsbInMap.insert({evt,NumDSBInForThisCluster});
0143             } else {
0144                 ndsbInMap[evt] += NumDSBInForThisCluster;
0145             }
0146
0147             int NumDSBdirInForThisCluster = damClass.GetNumDSBwithBothDirectIndirectDamage(classifiedDamage);
0148             if ( (ndsbdirIMap.find(evt) == ndsbdirIMap.end()) ) {
0149                 if (NumDSBdirInForThisCluster > 0) ndsbdirIMap.insert({evt,NumDSBdirInForThisCluster});
0150             } else {
0151                 ndsbdirIMap[evt] += NumDSBdirInForThisCluster;
0152             }
0153         }
0154     }
0155
0156     // DSB and its error:
0157     float xxtotal=0, rms, xtotal = 0;
0158     if (ndsbMap.size() >0) {
0159         for (auto const& [evt,numdsb] : ndsbMap) {
0160         xtotal += (float)numdsb;
0161         xxtotal += float(numdsb*numdsb);
0162         }
0163         if (ndsbMap.size() == 1) {
0164             // try to estimate error using poisson distribution
0165             rms = std::sqrt(xtotal);
0166         } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndsbMap.size()));
0167         fNDSBandError.first = xtotal;
0168         fNDSBandError.second = rms;
0169     }
0170
0171
0172     // cDSB and its error:
0173     if (ncdsbMap.size() > 0) {
0174         xxtotal=0, rms = 0, xtotal = 0;
0175         for (auto const& [evt,numcdsb] : ncdsbMap) {
0176             xtotal += (float)numcdsb;
0177             xxtotal += float(numcdsb*numcdsb);
0178         }
0179         if (ncdsbMap.size() == 1) {
0180             // try to estimate error using poisson distribution
0181             rms = std::sqrt(xtotal);
0182         } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ncdsbMap.size()));
0183         fNcDSBandError.first = xtotal;
0184         fNcDSBandError.second = rms;
0185     }
0186     // sDSB and its error, using error propagation method:
0187     if (fNDSBandError.first > 0) {
0188         fNsDSBandError.first = fNDSBandError.first -  fNcDSBandError.first;
0189         if (fNcDSBandError.first > 0) {
0190             fNsDSBandError.second = (fNsDSBandError.first)*std::sqrt(
0191                 (fNDSBandError.second/fNDSBandError.first)*(fNDSBandError.second/fNDSBandError.first) +
0192                 (fNcDSBandError.second/fNcDSBandError.first)*(fNcDSBandError.second/fNcDSBandError.first));
0193         }
0194         else fNsDSBandError.second = (fNsDSBandError.first)*(fNDSBandError.second/fNDSBandError.first);
0195     }
0196
0197     // DSBdir and its error:
0198     if (ndsbdirMap.size() > 0) {
0199         xxtotal=0, rms = 0, xtotal = 0;
0200         for (auto const& [evt,numdsbdir] : ndsbdirMap) {
0201             xtotal += (float)numdsbdir;
0202             xxtotal += float(numdsbdir*numdsbdir);
0203         }
0204         if (ndsbdirMap.size() == 1) {
0205             // try to estimate error using poisson distribution
0206             rms = std::sqrt(xtotal);
0207         } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndsbdirMap.size()));
0208         fNDSBdirandError.first = xtotal;
0209         fNDSBdirandError.second = rms;
0210     }
0211
0212     // DSBIn and its error:
0213     if (ndsbInMap.size() > 0) {
0214         xxtotal=0, rms = 0, xtotal = 0;
0215         for (auto const& [evt,numdsbIn] : ndsbInMap) {
0216             xtotal += (float)numdsbIn;
0217             xxtotal += float(numdsbIn*numdsbIn);
0218         }
0219         if (ndsbInMap.size() == 1) {
0220             // try to estimate error using poisson distribution
0221             rms = std::sqrt(xtotal);
0222         } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndsbInMap.size()));
0223         fNDSBIndandError.first = xtotal;
0224         fNDSBIndandError.second = rms;
0225     }
0226
0227     // DSBdirIn and its error:
0228     if (ndsbdirIMap.size() > 0) {
0229         xxtotal=0, rms = 0, xtotal = 0;
0230         for (auto const& [evt,numdsbdirIn] : ndsbdirIMap) {
0231             xtotal += (float)numdsbdirIn;
0232             xxtotal += float(numdsbdirIn*numdsbdirIn);
0233         }
0234         if (ndsbdirIMap.size() == 1) {
0235             // try to estimate error using poisson distribution
0236             rms = std::sqrt(xtotal);
0237         } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndsbdirIMap.size()));
0238         fNDSBdirIandError.first = xtotal;
0239         fNDSBdirIandError.second = rms;
0240     }
0241
0242
0243     // SSB and its error:
0244     if (nssbMap.size() > 0) {
0245         xxtotal=0, rms = 0, xtotal = 0;
0246         for (auto const& [evt,numssb] : nssbMap) {
0247             xtotal += (float)numssb;
0248             xxtotal += float(numssb*numssb);
0249         }
0250         if (nssbMap.size() == 1) {
0251             // try to estimate error using poisson distribution
0252             rms = std::sqrt(xtotal);
0253         } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(nssbMap.size()));
0254         fNSSBandError.first = xtotal;
0255         fNSSBandError.second = rms;
0256     }
0257
0258     // SB and its error:
0259     if (nsbMap.size() > 0) {
0260         xxtotal=0, rms = 0, xtotal = 0;
0261         for (auto const& [evt,numsb] : nsbMap) {
0262             xtotal += (float)numsb;
0263             xxtotal += float(numsb*numsb);
0264         }
0265         if (nsbMap.size() == 1) {
0266             // try to estimate error using poisson distribution
0267             rms = std::sqrt(xtotal);
0268         } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(nsbMap.size()));
0269         fNSBandError.first = xtotal;
0270         fNSBandError.second = rms;
0271     }
0272
0273     // direct SB and its error:
0274     if (ndirsbMap.size() > 0) {
0275         xxtotal=0, rms = 0, xtotal = 0;
0276         for (auto const& [evt,numdirsb] : ndirsbMap) {
0277             xtotal += (float)numdirsb;
0278             xxtotal += float(numdirsb*numdirsb);
0279         }
0280         if (ndirsbMap.size() == 1) {
0281             // try to estimate error using poisson distribution
0282             rms = std::sqrt(xtotal);
0283         } else rms = std::sqrt(std::fabs(xxtotal - xtotal*xtotal)/float(ndirsbMap.size()));
0284         fNdirSBandError.first = xtotal;
0285         fNdirSBandError.second = rms;
0286     }
0287
0288     // indirect SB its error, using error propagation method:
0289     if (fNSBandError.first > 0) {
0290         fNindirSBandError.first = fNSBandError.first -  fNdirSBandError.first;
0291         if (fNdirSBandError.first > 0) {
0292             fNindirSBandError.second = (fNindirSBandError.first)*std::sqrt(
0293                 (fNSBandError.second/fNSBandError.first)*(fNSBandError.second/fNSBandError.first) +
0294                 (fNdirSBandError.second/fNdirSBandError.first)*(fNdirSBandError.second/fNdirSBandError.first));
0295         }
0296         else fNindirSBandError.second = (fNindirSBandError.first)*(fNSBandError.second/fNSBandError.first);
0297     }
0298
0299     // clear Maps
0300     ndsbMap.clear();
0301     ncdsbMap.clear();
0302     nssbMap.clear();
0303     nsbMap.clear();
0304     ndirsbMap.clear();
0305     ndsbdirMap.clear();
0306     ndsbdirIMap.clear();
0307     ndsbInMap.clear();
0308     dmMap.clear();
0309
0310     fIsSBScanned = true;
0311 }
0312
0313 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0314
0315 void AnalysisHandler::GiveMeSBs()
0316 {
0317     if (!fIsSBScanned) GetAllDamageAndScanSB();
0318     std::string outName = ParametersParser::Instance()->GetOutputName();
0319     std::fstream file;
0320     file.open(outName.c_str(), std::ios_base::out);
0321     double norm = 1.0;
0322     std::string normunit = "";
0323     if (ParametersParser::Instance()->GetUnitTypeOfNormalization() == 2) {
0324         norm = 1.0/fDose;
0325         normunit = "[SB/Gy]";
0326     } else {
0327         double BbToGb = 1e-9; // convert Bb to Gb
0328         norm = 1.0/(fDose*fNBp*BbToGb);
0329         normunit = "[SB/Gy/Gbp]";
0330     }
0331     file <<"#=========================== Strand Breaks ============================#\n";
0332     file <<"Name of Cell Nucleus: "<<ParametersParser::Instance()->GetCellNucleusName()<<"\n";
0333     if (ParametersParser::Instance()->WannaLoadDamagesFromSDD()) {
0334         std::string outstrtmp = "No info from SDD file!!!\n";
0335         file <<"Volume of Cell Nucleus: "<<outstrtmp;
0336         file <<"Mass Density of Cell Nucleus: "<<outstrtmp;
0337         file <<"Mass of Cell Nucleus: "<<outstrtmp;
0338         file <<"Energy deposited in Cell Nucleus: "<<outstrtmp;
0339         file <<"Dose delivered in Cell Nucleus: "<<fDose <<" (Gy)\n";
0340         file <<"Minimum Distance between two clusters: "<<fBpForDSB <<" (bp)\n";
0341         file <<"Number of basepairs in Cell Nucleus: "<<fNBp <<" (bp)\n";
0342         file <<"Threshold Energy for direct damage selection: "<<outstrtmp;
0343         file <<"Propability for indirect damage selection: "<<outstrtmp;
0344     } else {
0345         file <<"Volume of Cell Nucleus: "<<fScanDamage->GetNucleusVolume()<<" (m3)\n";
0346         file <<"Mass Density of Cell Nucleus: "<<fScanDamage->GetNucleusMassDensity()<<" (kg/m3)\n";
0347         file <<"Mass of Cell Nucleus: "<<fScanDamage->GetNucleusMass()<<" (kg)\n";
0348         file <<"Energy deposited in Cell Nucleus: "<<fEdepInNucleus <<" (eV)\n";
0349         file <<"Dose delivered in Cell Nucleus: "<<fDose <<" (Gy)\n";
0350         file <<"Minimum Distance between two clusters: "<<fBpForDSB <<" (bp)\n";
0351         file <<"Number of basepairs in Cell Nucleus: "<<fNBp <<" (bp)\n";
0352         file <<"Threshold Energy for direct damage selection: "<<fScanDamage->GetThresholdEnergy() <<" (eV)\n";
0353         if (fScanDamage->SkippedScanningIndirectDamage()) file <<"Propability for indirect SB selection: "
0354                                                                 <<" Skipped the indirect analysis\n";
0355         else file <<"Propability for indirect damage selection: "
0356                 <<fScanDamage->GetProbabilityForIndirectSBSelection()*100.<<" (%)\n";
0357     }
0358
0359     file <<"#======================================================================#\n";
0360     file << "\n";
0361     file <<"#Un-normalized results:\n";
0362     file << "TotalSB  [SB]   \t" << fNSBandError.first <<"\t+/-\t"<<fNSBandError.second<< "\n";
0363     file << "DirSB    [SB]   \t" << fNdirSBandError.first <<"\t+/-\t"<<fNdirSBandError.second<< "\n";
0364     file << "IndirSB  [SB]   \t" << fNindirSBandError.first <<"\t+/-\t"<<fNindirSBandError.second<< "\n";
0365     file << "SSB      [SB]   \t" << fNSSBandError.first <<"\t+/-\t"<<fNSSBandError.second<< "\n";
0366     file << "DSB      [SB]   \t" << fNDSBandError.first <<"\t+/-\t"<<fNDSBandError.second<< "\n";
0367     file << "cDSB     [SB]   \t" << fNcDSBandError.first <<"\t+/-\t"<<fNcDSBandError.second<< "\n";
0368     file << "sDSB     [SB]   \t" << fNsDSBandError.first <<"\t+/-\t"<<fNsDSBandError.second<< "\n";
0369     file << "DSBdir   [SB]   \t" << fNDSBdirandError.first <<"\t+/-\t"<<fNDSBdirandError.second<< "\n";
0370     file << "DSBind   [SB]   \t" << fNDSBIndandError.first <<"\t+/-\t"<<fNDSBIndandError.second<< "\n";
0371     file << "DSBdirIn [SB]   \t" << fNDSBdirIandError.first <<"\t+/-\t"<<fNDSBdirIandError.second<< "\n";
0372     file << "\n";
0373     file <<"#Normalized results:\n";
0374     file << "TotalSB  " + normunit +"    \t" << fNSBandError.first * norm <<"\t+/-\t"<<fNSBandError.second * norm<< "\n";
0375     file << "DirSB    " + normunit +"    \t" << fNdirSBandError.first * norm<<"\t+/-\t"<<fNdirSBandError.second * norm<< "\n";
0376     file << "IndirSB  " + normunit +"    \t" << fNindirSBandError.first * norm<<"\t+/-\t"<<fNindirSBandError.second * norm<< "\n";
0377     file << "SSB      " + normunit +"    \t" << fNSSBandError.first * norm<<"\t+/-\t"<<fNSSBandError.second * norm<< "\n";
0378     file << "DSB      " + normunit +"    \t" << fNDSBandError.first * norm<<"\t+/-\t"<<fNDSBandError.second * norm<< "\n";
0379     file << "cDSB     " + normunit +"    \t" << fNcDSBandError.first * norm<<"\t+/-\t"<<fNcDSBandError.second * norm<< "\n";
0380     file << "sDSB     " + normunit +"    \t" << fNsDSBandError.first * norm<<"\t+/-\t"<<fNsDSBandError.second * norm<< "\n";
0381     file << "DSBdir   " + normunit +"    \t" << fNDSBdirandError.first * norm<<"\t+/-\t"<<fNDSBdirandError.second * norm<< "\n";
0382     file << "DSBind   " + normunit +"    \t" << fNDSBIndandError.first * norm<<"\t+/-\t"<<fNDSBIndandError.second * norm<< "\n";
0383     file << "DSBdirIn " + normunit +"    \t" << fNDSBdirIandError.first * norm<<"\t+/-\t"<<fNDSBdirIandError.second * norm<< "\n";
0384     file <<"#======================================================================#\n";
0385     file << "where: \n";
0386     file << "-----> TotalSB: Total strand-breaks\n";
0387     file << "-----> DirSB: Direct strand-breaks\n";
0388     file << "-----> IndirSB: Indirect strand-breaks\n";
0389     file << "-----> SSB: Single strand-breaks\n";
0390     file << "-----> DSB: Double strand-breaks\n";
0391     file << "-----> cDSB: Complex DSB\n";
0392     file << "-----> sDSB: Simple DSB\n";
0393     file << "-----> DSBdir: DSB that contains at least one direct SB\n";
0394     file << "-----> DSBdind: DSB that contains at least one indirect SB\n";
0395     file << "-----> DSBdirIn: DSB that contains at both direct and indirect SB\n";
0396     file <<"#============================== End ===================================#\n";
0397     file.close();
0398 }
0399
0400 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0401
0402 void AnalysisHandler::ApplyDNAModel(std::string dnaModel)
0403 {
0404     if (!fIsSBScanned) GetAllDamageAndScanSB();
0405
0406     if (dnaModel == "TLK") {
0407         std::cout << "Invoking TLK Model" << std::endl;
0408         //fTLKModel->SetDose(fDose);
0409         SetParametersForTLKModel();
0410         //fTLKModel->ComputeAndSetDamageInput(fAllDamage);
0411         fTLKModel->SetSingleDSBYield(fNsDSBandError.first/fDose);
0412         fTLKModel->SetComplexDSBYield(fNcDSBandError.first/fDose);
0413         if (ParametersParser::Instance()->GetTLKdoseMax() != "") {
0414             auto val = std::stod(ParametersParser::Instance()->GetTLKdoseMax());
0415             if (val != pTLKDoseMax) pTLKDoseMax = val;
0416         }
0417         if (ParametersParser::Instance()->GetTLKdeltaDose() != "") {
0418             auto val = std::stod(ParametersParser::Instance()->GetTLKdeltaDose());
0419             if (val != pTLKDeltaDose) pTLKDeltaDose = val;
0420         }
0421         fTLKModel->CalculateRepair(pTLKDoseMax,pTLKDeltaDose);
0422         std::string outname = "TLK_"+ParametersParser::Instance()->GetOutputName();
0423         fTLKModel->WriteOutput(outname);
0424     }
0425
0426     if (dnaModel == "LEMIV") {
0427         std::cout << "Invoking LEMIV Model" << std::endl;
0428         fLEMIVModel->SetChromosomeBpSizesMap(fChromosomeBpMap);
0429         fLEMIVModel->SetDose(fDose);
0430         SetParametersForLEMIVModel();
0431         fLEMIVModel->ComputeAndSetDamageInput(fAllDamage);
0432         if (ParametersParser::Instance()->GetLEMtimeMax() != "") {
0433             auto val = std::stod(ParametersParser::Instance()->GetLEMtimeMax());
0434             if (val != pLEMIVtimeMax) pLEMIVtimeMax = val;
0435         }
0436         if (ParametersParser::Instance()->GetLEMdeltaTime() != "") {
0437             auto val = std::stod(ParametersParser::Instance()->GetLEMdeltaTime());
0438             if (val != pLEMIVdeltaTime) pLEMIVdeltaTime = val;
0439         }
0440         fLEMIVModel->CalculateRepair(pLEMIVtimeMax,pLEMIVdeltaTime);
0441         std::string outname = "LEMIV_"+ParametersParser::Instance()->GetOutputName();
0442         fLEMIVModel->WriteOutput(outname);
0443     }
0444
0445     if (dnaModel == "BELOV") {
0446         std::cout << "Invoking Belov's Model" << std::endl;
0447         fBelovModel->SetDSBandComDSBandDose(fNDSBandError.first,fNcDSBandError.first,fDose);
0448         if (ParametersParser::Instance()->GetBELOVNirrep() != "") {
0449             auto Nirrep = std::stod(ParametersParser::Instance()->GetBELOVNirrep());
0450             fBelovModel->SetNirrep(Nirrep);
0451         }
0452         double Dz = 1.0;
0453         if (ParametersParser::Instance()->GetBELOVDz() != "") {
0454             Dz = std::stod(ParametersParser::Instance()->GetBELOVDz());
0455         }
0456         fBelovModel->CalculateRepair(Dz);
0457         std::string outname = "BELOV_"+ParametersParser::Instance()->GetOutputName();
0458         fBelovModel->WriteOutput(outname);
0459     }
0460 }
0461
0462 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0463
0464 void AnalysisHandler::CreateSDD(std::string filename)
0465 {
0466     std::string str_tmp;
0467     std::fstream ofile;
0468     ofile.open(filename.c_str(), std::ios_base::out);
0469     if (ofile.is_open()) {
0470         // Create header
0471         ofile
0472             <<"SDD version, SDDv1.0;\n"
0473             <<"Software, dsbandrepair;\n"
0474             <<"Author contact, Le Tuan Anh - anh.letuan@irsn.fr, , ;\n"
0475             <<"Simulation Details, DNA damages from direct and indirect effects;\n"
0476             <<"Source, ;\n"
0477             <<"Source type, ;\n"
0478             <<"Incident particles, "<<0<<";\n"
0479             <<"Mean particle energy ("<<ParametersParser::Instance()->GetEnergyUnit()<<"), "
0480             <<ParametersParser::Instance()->GetParticleEnergy()<<";\n"
0481             <<"Energy distribution, , ;\n"
0482             <<"Particle fraction, 0;\n"
0483             <<"Dose or fluence, 1, "<<fDose<<";\n"
0484             <<"Dose rate, 0;\n"
0485             <<"Irradiation target, ;\n"
0486             <<"Volumes, 0;\n";
0487         ofile<<"Chromosome sizes, "<<fChromosomeBpMap.size();
0488         for (auto const& [chroID, nBps] :fChromosomeBpMap) {
0489             float nMBps = nBps*1E-6;// convert from Bp to MBp
0490             ofile<<", "<<nMBps;
0491         }
0492         ofile<<";\n";
0493         ofile<<"DNA Density, 0;\n"
0494             <<"Cell Cycle Phase, 0;\n"
0495             <<"DNA Structure, 0;\n"
0496             <<"In vitro / in vivo, ;\n"
0497             <<"Proliferation status, ;\n"
0498             <<"Microenvironment, 0, 0;\n"
0499             <<"Damage definition, 0;\n"
0500             <<"Time, 0;\n"
0501             <<"Damage and primary count, "+std::to_string(fAllDamage.size())+", 0;\n"
0502             <<"Data entries, 1, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0;\n"
0503             <<"Data field explaination, Field 1: [1]-eventID, Field 3: [0]-Chromatin "
0504             <<"type [1]-ChromosomeID [3]-strand, Field 4:chrom position (copynb), Field 5: "
0505             <<"Cause (direct: [0]=0)  (indirect: [0]=1), Field 6: Damage types (Base:[0]>0) (Backbone: [1]>0);\n"
0506             <<"\n"
0507             <<"***EndOfHeader***;\n"
0508             <<"\n";
0509
0510         // Data Section
0511         int prevEvt = -1;
0512         for (auto &damage : fAllDamage) {
0513             //Field 1 Calassification
0514             int newEvtFlag = 0; // = 2 if new event;
0515             if (prevEvt != damage.GetEvt()) {
0516                 newEvtFlag = 2;
0517                 prevEvt = damage.GetEvt();
0518             }
0519             ofile<<newEvtFlag<<", "<<damage.GetEvt()<<"; ";
0520             //Field 3 Chromosome IDs
0521             ofile<<damage.GetDamageChromatin()<<", "<<damage.GetChromo()<<", "<<0<<", "<<damage.GetStrand()<<"; ";
0522             //Field 4, Chromosome position
0523             ofile<<damage.GetCopyNb()<<"; ";
0524             //Field 5, Cause: Unknown = -1, Direct = 0, Indirect = 1
0525             ofile<<damage.GetCause()<<", "<<0<<", "<<0<<"; ";
0526             //Field 6, Damage types:
0527             int firstval = 0, secval = 0;
0528             if (damage.GetDamageType() == Damage::DamageType::fBase) firstval = 1;
0529             if (damage.GetDamageType() == Damage::DamageType::fBackbone) secval = 1;
0530             ofile<<firstval<<", "<<secval<<", "<<0<<"; ";
0531             ofile<<"\n";
0532         }
0533     }
0534     ofile.close();
0535 }
0536
0537 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0538
0539 void AnalysisHandler::SetBpForDSB(unsigned int pVal)
0540 {
0541     if (pVal == fBpForDSB) return;
0542     fBpForDSB = pVal;
0543     fTLKModel->SetBpForDSB(fBpForDSB);
0544     fLEMIVModel->SetBpForDSB(fBpForDSB);
0545     fBelovModel->SetBpForDSB(fBpForDSB);
0546 }
0547
0548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0549
0550 void AnalysisHandler::SetParametersForTLKModel(double pLambda1,double pLambda2, double pBeta1, double pBeta2,double pEta)
0551 {
0552     if (ParametersParser::Instance()->GetTLKLambda1() != "")
0553         pLambda1 = std::stod(ParametersParser::Instance()->GetTLKLambda1());
0554     if (ParametersParser::Instance()->GetTLKLambda2() != "")
0555         pLambda2 = std::stod(ParametersParser::Instance()->GetTLKLambda2());
0556     if (ParametersParser::Instance()->GetTLKBeta1() != "")
0557         pBeta1 = std::stod(ParametersParser::Instance()->GetTLKBeta1());
0558     if (ParametersParser::Instance()->GetTLKBeta2() != "")
0559         pBeta2 = std::stod(ParametersParser::Instance()->GetTLKBeta2());
0560     if (ParametersParser::Instance()->GetTLKEta() != "")
0561         pEta = std::stod(ParametersParser::Instance()->GetTLKEta());
0562     if (pBeta1 != fTLKModel->GetBeta1()) fTLKModel->SetBeta1(pBeta1);
0563     if (pBeta2 != fTLKModel->GetBeta2()) fTLKModel->SetBeta2(pBeta2);
0564     if (pLambda1 != fTLKModel->GetLambda1()) fTLKModel->SetLambda1(pLambda1);
0565     if (pLambda2 != fTLKModel->GetLambda2()) fTLKModel->SetLambda2(pLambda2);
0566     if (pEta != fTLKModel->GetEta()) fTLKModel->SetEta(pEta);
0567 }
0568
0569 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
0570
0571 void AnalysisHandler::SetParametersForLEMIVModel(double pLoopLength,double pFunrej,double pTfast,double pTslow)
0572 {
0573     if (ParametersParser::Instance()->GetEMIVLoopLength() != "")
0574         pLoopLength = std::stod(ParametersParser::Instance()->GetEMIVLoopLength());
0575     if (ParametersParser::Instance()->GetEMIVFunrej() != "")
0576         pFunrej = std::stod(ParametersParser::Instance()->GetEMIVFunrej());
0577     if (ParametersParser::Instance()->GetEMIVTFast() != "")
0578         pTfast = std::stod(ParametersParser::Instance()->GetEMIVTFast());
0579     if (ParametersParser::Instance()->GetEMIVTSlow() != "")
0580         pTslow = std::stod(ParametersParser::Instance()->GetEMIVTSlow());
0581
0582     if (pLoopLength != fLEMIVModel->GetLoopLength()) fLEMIVModel->SetLoopLength(pLoopLength);
0583     if (pFunrej != fLEMIVModel->GetFunrej()) fLEMIVModel->SetFunrej(pFunrej);
0584     if (pTfast != fLEMIVModel->GetTfast()) fLEMIVModel->SetTfast(pTfast);
0585     if (pTslow != fLEMIVModel->GetTslow()) fLEMIVModel->SetTslow(pTslow);
0586 }
0587
0588 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....