Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/medical/dna/pdb4dna/src/PDBlib.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/medical/dna/pdb4dna/src/PDBlib.cc (Version 11.3.0) and /examples/extended/medical/dna/pdb4dna/src/PDBlib.cc (Version 11.0.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // This example is provided by the Geant4-DNA      26 // This example is provided by the Geant4-DNA collaboration
 27 // Any report or published results obtained us     27 // Any report or published results obtained using the Geant4-DNA software
 28 // shall cite the following Geant4-DNA collabo     28 // shall cite the following Geant4-DNA collaboration publication:
 29 // Med. Phys. 37 (2010) 4692-4708                  29 // Med. Phys. 37 (2010) 4692-4708
 30 // Delage et al. PDB4DNA: implementation of DN     30 // Delage et al. PDB4DNA: implementation of DNA geometry from the Protein Data
 31 //                  Bank (PDB) description for     31 //                  Bank (PDB) description for Geant4-DNA Monte-Carlo
 32 //                  simulations (submitted to      32 //                  simulations (submitted to Comput. Phys. Commun.)
 33 // The Geant4-DNA web site is available at htt     33 // The Geant4-DNA web site is available at http://geant4-dna.org
 34 //                                                 34 //
 35 //                                             <<  35 // 
 36 /// \file PDBlib.cc                                36 /// \file PDBlib.cc
 37 /// \brief Implementation file for PDBlib clas     37 /// \brief Implementation file for PDBlib class
 38                                                    38 
 39 #include "PDBlib.hh"                               39 #include "PDBlib.hh"
 40                                                    40 
 41 // define if the program is running with Geant <<  41 //define if the program is running with Geant4
 42 #define GEANT4                                     42 #define GEANT4
 43                                                    43 
 44 #ifdef GEANT4                                      44 #ifdef GEANT4
 45 // Specific to Geant4, globals.hh is used for  <<  45 //Specific to Geant4, globals.hh is used for G4cout
 46 #  include "globals.hh"                        <<  46 #include "globals.hh"
 47 #else                                              47 #else
 48 #  define G4cout std::cout                     <<  48 #define G4cout std::cout
 49 #  define G4cerr std::cerr                     <<  49 #define G4cerr std::cerr 
 50 #  define G4endl std::endl                     <<  50 #define G4endl std::endl 
 51 #  define G4String std::string  // string incl <<  51 #define G4String std::string // string included in PDBatom.hh
 52 #  include <cfloat>                            <<  52 #include <cfloat>
 53 #endif                                             53 #endif
 54 #include <cmath>                               << 
 55 #include <fstream>                                 54 #include <fstream>
 56 #include <iostream>                                55 #include <iostream>
 57 #include <limits>                                  56 #include <limits>
                                                   >>  57 #include <cmath>
 58 #include <sstream>                                 58 #include <sstream>
 59 #include <stdlib.h>                                59 #include <stdlib.h>
 60                                                    60 
 61 using namespace std;                               61 using namespace std;
 62                                                    62 
 63 //....oooOO0OOooo........oooOO0OOooo........oo     63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 64                                                    64 
 65 PDBlib::PDBlib() : fNbNucleotidsPerStrand(0) { <<  65 PDBlib::PDBlib() :
                                                   >>  66     fNbNucleotidsPerStrand(0)
                                                   >>  67 {
                                                   >>  68 }
 66                                                    69 
 67 //....oooOO0OOooo........oooOO0OOooo........oo     70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 68                                                    71 
 69 /**                                                72 /**
 70  * \brief    Load a PDB file into memory           73  * \brief    Load a PDB file into memory
 71  * \details  molecule (polymer,?),             <<  74  * \details  molecule (polymer,?), 
 72  *                  read line, key words           75  *                  read line, key words
 73  * \param    filename    G4String for filename     76  * \param    filename    G4String for filename
 74  * \param    isDNA                                 77  * \param    isDNA
 75  * \param    verbose                               78  * \param    verbose
 76  * \return   List of molecules                     79  * \return   List of molecules
 77  */                                                80  */
 78                                                    81 
 79 Molecule* PDBlib::Load(const std::string& file <<  82 Molecule * PDBlib::Load(const std::string& filename,
 80                        unsigned short int verb <<  83                         unsigned short int& isDNA,
                                                   >>  84                         unsigned short int verbose = 0)
 81 {                                                  85 {
 82   G4String sLine = "";                             86   G4String sLine = "";
 83   ifstream infile;                                 87   ifstream infile;
 84   infile.open(filename.c_str());                   88   infile.open(filename.c_str());
 85   if (!infile) {                               <<  89   if(!infile)
 86     G4cout << "PDBlib::load >> file " << filen <<  90   {
                                                   >>  91     G4cout << "PDBlib::load >> file " << filename << " not found !!!!"
                                                   >>  92     << G4endl;
 87   }                                                93   }
 88   else {                                       <<  94   else
 89     int nbAtom = 0;  // total number of atoms  <<  95   {
 90     int numAtomInRes = 0;  // number of an ato <<  96     int nbAtomTot = 0; // total number of atoms
 91     int nbResidue = 0;  // total number of res <<  97     int nbAtom = 0; // total number of atoms in a residue
 92     int nbMolecule = 0;  // total number of mo <<  98     int numAtomInRes = 0; // number of an atom in a residue
 93                                                <<  99     int nbResidue = 0; // total number of residues
 94     Atom* AtomicOld = nullptr;                 << 100     int nbMolecule = 0; // total number of molecule
 95     Atom* AtomicNext = nullptr;                << 101 
 96                                                << 102     Atom * AtomicOld = nullptr;
 97     int serial;  // Atom serial number         << 103     Atom * AtomicNext = nullptr;
 98     G4String atomName;  // Atom name           << 104 
 99     G4String element;  // Element Symbol       << 105     int serial; //Atom serial number
100     G4String resName;  // Residue name for thi << 106     G4String atomName; //Atom name
101     double x, y, z;  // Orthogonal coordinates << 107     G4String element; //Element Symbol
102     double occupancy;  // Occupancy            << 108     G4String resName; //Residue name for this atom
103                                                << 109     double x, y, z; //Orthogonal coordinates in Angstroms
104     Residue* residueOld = nullptr;             << 110     double occupancy; //Occupancy
105     Residue* residueFirst = nullptr;           << 111 
106     Residue* residueNext = nullptr;            << 112     Residue * residueOld = nullptr;
107                                                << 113     Residue * residueFirst = nullptr;
108     Molecule* moleculeFirst = nullptr;         << 114     Residue * residueNext = nullptr;
109     Molecule* moleculeOld = nullptr;           << 115 
110     Molecule* moleculeNext = nullptr;          << 116     Molecule * moleculeFirst = nullptr;
                                                   >> 117     Molecule * moleculeOld = nullptr;
                                                   >> 118     Molecule * moleculeNext = nullptr;
111                                                   119 
112     /////////////////////////////////             120     /////////////////////////////////
113     // NEW variable to draw a fitting cylinder << 121     //NEW variable to draw a fitting cylinder if z oriented
114     //=> fitting box                              122     //=> fitting box
115     double minGlobZ, maxGlobZ;                    123     double minGlobZ, maxGlobZ;
116     double minGlobX, maxGlobX;                    124     double minGlobX, maxGlobX;
117     double minGlobY, maxGlobY;                    125     double minGlobY, maxGlobY;
118     double minX, maxX, minY, maxY, minZ, maxZ; << 126     double minX, maxX, minY, maxY, minZ, maxZ; //Sort of 'mother volume' box
119                                                   127 
120     minGlobZ = -DBL_MAX;                          128     minGlobZ = -DBL_MAX;
121     minGlobX = -DBL_MAX;                          129     minGlobX = -DBL_MAX;
122     minGlobY = -DBL_MAX;                          130     minGlobY = -DBL_MAX;
123     maxGlobZ = DBL_MAX;                           131     maxGlobZ = DBL_MAX;
124     maxGlobX = DBL_MAX;                           132     maxGlobX = DBL_MAX;
125     maxGlobY = DBL_MAX;                           133     maxGlobY = DBL_MAX;
126     minX = -DBL_MAX;                              134     minX = -DBL_MAX;
127     minY = -DBL_MAX;                              135     minY = -DBL_MAX;
128     minZ = -DBL_MAX;                              136     minZ = -DBL_MAX;
129     maxX = DBL_MAX;                               137     maxX = DBL_MAX;
130     maxY = DBL_MAX;                               138     maxY = DBL_MAX;
131     maxZ = DBL_MAX;                               139     maxZ = DBL_MAX;
132                                                   140 
133     int lastResSeq = -1;                          141     int lastResSeq = -1;
134     int resSeq = 0;                               142     int resSeq = 0;
135                                                   143 
136     G4String nameMol;                             144     G4String nameMol;
137                                                   145 
138     unsigned short int numModel = 0;              146     unsigned short int numModel = 0;
139     unsigned short int model = 0;                 147     unsigned short int model = 0;
140                                                   148 
141     // Number of TER                           << 149     //Number of TER
142     int ter = 0;                                  150     int ter = 0;
143     // Number of TER (chain) max for this file << 151     //Number of TER (chain) max for this file
144                                                   152 
145     int terMax = INT_MAX;                         153     int terMax = INT_MAX;
146                                                   154 
147     if (!infile.eof()) {                       << 155     if(!infile.eof())
                                                   >> 156     {
148       getline(infile, sLine);                     157       getline(infile, sLine);
149       std::size_t found = sLine.find("DNA");      158       std::size_t found = sLine.find("DNA");
150       if (found != G4String::npos) {           << 159       if(found != G4String::npos)
                                                   >> 160       {
151         terMax = 2;                               161         terMax = 2;
152         isDNA = 1;                                162         isDNA = 1;
153       }                                           163       }
154       else                                     << 164       else isDNA = 0;
155         isDNA = 0;                             << 165       //If PDB file have not a header line
156       // If PDB file have not a header line    << 
157       found = sLine.find("HEADER");               166       found = sLine.find("HEADER");
158       if (found == G4String::npos) {           << 167       if(found == G4String::npos)
                                                   >> 168       {
159         infile.close();                           169         infile.close();
160         infile.open(filename.c_str());            170         infile.open(filename.c_str());
161                                                   171 
162         G4cout << "PDBlib::load >> No header f    172         G4cout << "PDBlib::load >> No header found !!!!" << G4endl;
163       }                                           173       }
164     }                                             174     }
165                                                   175 
166     while (!infile.eof()) {                    << 176     while(!infile.eof())
                                                   >> 177     {
167       getline(infile, sLine);                     178       getline(infile, sLine);
168       // In the case of several models         << 179       //In the case of several models
169       if ((sLine.substr(0, 6)).compare("NUMMDL << 180       if((sLine.substr(0, 6)).compare("NUMMDL") == 0)
                                                   >> 181       {
170         istringstream((sLine.substr(10, 4))) >    182         istringstream((sLine.substr(10, 4))) >> numModel;
171       }                                           183       }
172       if ((numModel > 0) && ((sLine.substr(0,  << 184       if((numModel > 0) && ((sLine.substr(0, 6)).compare("MODEL ") == 0))
                                                   >> 185       {
173         istringstream((sLine.substr(10, 4))) >    186         istringstream((sLine.substr(10, 4))) >> model;
174         //////////////////////////////////////    187         //////////////////////////////////////////
175         if (model > 1) break;                  << 188         if(model > 1) break;
176         //////////////////////////////////////    189         //////////////////////////////////////////
177       }                                           190       }
178       // For verification of residue sequence  << 191       //For verification of residue sequence
179       if ((sLine.substr(0, 6)).compare("SEQRES << 192       if((sLine.substr(0, 6)).compare("SEQRES") == 0)
180         // Create list of molecule here        << 193       {
181         if (verbose > 1) G4cout << sLine << G4 << 194         //Create list of molecule here
                                                   >> 195         if(verbose > 1) G4cout << sLine << G4endl;
182       }                                           196       }
183                                                   197 
184       // Coordinate section                    << 198       //Coordinate section
185       if ((numModel > 0) && ((sLine.substr(0,  << 199       if((numModel > 0) && ((sLine.substr(0, 6)).compare("ENDMDL") == 0))
                                                   >> 200       {
186         ;                                         201         ;
187       }                                           202       }
188       else if ((sLine.substr(0, 6)).compare("T << 203       else if((sLine.substr(0, 6)).compare("TER   ") == 0) //3 spaces
189       {                                           204       {
190         //////////////////////////////////////    205         //////////////////////////////////////////
191         // Currently retrieve only the two fir << 206         //Currently retrieve only the two first chains(TER) => two DNA strands
192         ter++;                                    207         ter++;
193         if (ter > terMax) break;               << 208         if(ter > terMax) break;
194         //////////////////////////////////////    209         //////////////////////////////////////////
195                                                   210 
196         // if (verbose>1) G4cout << sLine << G << 211         //if (verbose>1) G4cout << sLine << G4endl;
197         /************ Begin TER **************    212         /************ Begin TER ******************/
198         lastResSeq = -1;                          213         lastResSeq = -1;
199         resSeq = 0;                               214         resSeq = 0;
200                                                   215 
201         AtomicOld->SetNext(NULL);                 216         AtomicOld->SetNext(NULL);
202         residueOld->SetNext(NULL);                217         residueOld->SetNext(NULL);
203         residueOld->fNbAtom = nbAtom;             218         residueOld->fNbAtom = nbAtom;
204                                                   219 
205         // Molecule creation:                  << 220         //Molecule creation:
206         if (moleculeOld == NULL) {             << 221         if(moleculeOld == NULL)
207           nameMol = filename;  //+numModel     << 222         {
                                                   >> 223           nameMol = filename; //+numModel
208           moleculeOld = new Molecule(nameMol,     224           moleculeOld = new Molecule(nameMol, nbMolecule);
209           moleculeOld->SetFirst(residueFirst);    225           moleculeOld->SetFirst(residueFirst);
210           moleculeOld->fNbResidue = nbResidue;    226           moleculeOld->fNbResidue = nbResidue;
211           moleculeFirst = moleculeOld;            227           moleculeFirst = moleculeOld;
212         }                                         228         }
213         else {                                 << 229         else
                                                   >> 230         {
214           moleculeNext = new Molecule(nameMol,    231           moleculeNext = new Molecule(nameMol, nbMolecule);
215           moleculeOld->SetNext(moleculeNext);     232           moleculeOld->SetNext(moleculeNext);
216           moleculeNext->SetFirst(residueFirst)    233           moleculeNext->SetFirst(residueFirst);
217           moleculeNext->fNbResidue = nbResidue    234           moleculeNext->fNbResidue = nbResidue;
218           moleculeOld = moleculeNext;             235           moleculeOld = moleculeNext;
219         }                                         236         }
220                                                   237 
221         nbMolecule++;                             238         nbMolecule++;
222         moleculeOld->SetNext(NULL);               239         moleculeOld->SetNext(NULL);
223         moleculeOld->fCenterX = (int)((minX +  << 240         moleculeOld->fCenterX = (int) ((minX + maxX) / 2.);
224         moleculeOld->fCenterY = (int)((minY +  << 241         moleculeOld->fCenterY = (int) ((minY + maxY) / 2.);
225         moleculeOld->fCenterZ = (int)((minZ +  << 242         moleculeOld->fCenterZ = (int) ((minZ + maxZ) / 2.);
226         moleculeOld->fMaxGlobZ = maxGlobZ;        243         moleculeOld->fMaxGlobZ = maxGlobZ;
227         moleculeOld->fMinGlobZ = minGlobZ;        244         moleculeOld->fMinGlobZ = minGlobZ;
228         moleculeOld->fMaxGlobX = maxGlobX;        245         moleculeOld->fMaxGlobX = maxGlobX;
229         moleculeOld->fMinGlobX = minGlobX;        246         moleculeOld->fMinGlobX = minGlobX;
230         moleculeOld->fMaxGlobY = maxGlobY;        247         moleculeOld->fMaxGlobY = maxGlobY;
231         moleculeOld->fMinGlobY = minGlobY;        248         moleculeOld->fMinGlobY = minGlobY;
232                                                   249 
233         minGlobZ = -DBL_MAX;                      250         minGlobZ = -DBL_MAX;
234         minGlobX = -DBL_MAX;                      251         minGlobX = -DBL_MAX;
235         minGlobY = -DBL_MAX;                      252         minGlobY = -DBL_MAX;
236         maxGlobZ = DBL_MAX;                       253         maxGlobZ = DBL_MAX;
237         maxGlobX = DBL_MAX;                       254         maxGlobX = DBL_MAX;
238         maxGlobY = DBL_MAX;                       255         maxGlobY = DBL_MAX;
239         minX = -DBL_MAX;                          256         minX = -DBL_MAX;
240         minY = -DBL_MAX;                          257         minY = -DBL_MAX;
241         minZ = -DBL_MAX;                          258         minZ = -DBL_MAX;
242         maxX = DBL_MAX;                           259         maxX = DBL_MAX;
243         maxY = DBL_MAX;                           260         maxY = DBL_MAX;
244         maxZ = DBL_MAX;                           261         maxZ = DBL_MAX;
245                                                   262 
246         nbAtom = 0;                               263         nbAtom = 0;
247         numAtomInRes = 0;                         264         numAtomInRes = 0;
248         nbResidue = 0;                            265         nbResidue = 0;
249         AtomicOld = NULL;                         266         AtomicOld = NULL;
250         AtomicNext = NULL;                        267         AtomicNext = NULL;
251         residueOld = NULL;                        268         residueOld = NULL;
252         residueFirst = NULL;                      269         residueFirst = NULL;
253         residueNext = NULL;                       270         residueNext = NULL;
254                                                   271 
255         ///////////// End TER ////////////////    272         ///////////// End TER ///////////////////
256       }                                           273       }
257       else if ((sLine.substr(0, 6)).compare("A << 274       else if((sLine.substr(0, 6)).compare("ATOM  ") == 0)
                                                   >> 275       {
                                                   >> 276 
258         /************ Begin ATOM *************    277         /************ Begin ATOM ******************/
259         // serial                              << 278         //serial
260         istringstream((sLine.substr(6, 5))) >>    279         istringstream((sLine.substr(6, 5))) >> serial;
261                                                   280 
262         // To be improved                      << 281         //To be improved
263         // atomName :                          << 282         //atomName :
264         atomName = sLine.substr(12, 4);           283         atomName = sLine.substr(12, 4);
265         if (atomName.substr(0, 1).compare(" ") << 284         if(atomName.substr(0, 1).compare(" ") == 0) element = sLine.substr(13,
266           element = sLine.substr(13, 1);       << 285                                                                            1);
267         else                                   << 286         else element = sLine.substr(12, 1);
268           element = sLine.substr(12, 1);       << 
269                                                   287 
270         // set Van der Waals radius expressed     288         // set Van der Waals radius expressed in Angstrom
271         double vdwRadius = -1.;                   289         double vdwRadius = -1.;
272         if (element == "H") {                  << 290         if(element == "H")
                                                   >> 291         {
273           vdwRadius = 1.2;                        292           vdwRadius = 1.2;
274         }                                         293         }
275         else if (element == "C") {             << 294         else if(element == "C")
                                                   >> 295         {
276           vdwRadius = 1.7;                        296           vdwRadius = 1.7;
277         }                                         297         }
278         else if (element == "N") {             << 298         else if(element == "N")
                                                   >> 299         {
279           vdwRadius = 1.55;                       300           vdwRadius = 1.55;
280         }                                         301         }
281         else if (element == "O") {             << 302         else if(element == "O")
                                                   >> 303         {
282           vdwRadius = 1.52;                       304           vdwRadius = 1.52;
283         }                                         305         }
284         else if (element == "P") {             << 306         else if(element == "P")
                                                   >> 307         {
285           vdwRadius = 1.8;                        308           vdwRadius = 1.8;
286         }                                         309         }
287         else if (element == "S") {             << 310         else if(element == "S")
                                                   >> 311         {
288           vdwRadius = 1.8;                        312           vdwRadius = 1.8;
289         }                                         313         }
290         else {                                 << 314         else
291 #ifndef GEANT4                                 << 315         {
                                                   >> 316 #ifndef GEANT4   
292           G4cerr << "Element not recognized :     317           G4cerr << "Element not recognized : " << element << G4endl;
293           G4cerr << "Stop now" << G4endl;         318           G4cerr << "Stop now" << G4endl;
294           exit(1);                                319           exit(1);
295 #else                                             320 #else
296           G4ExceptionDescription errMsg;          321           G4ExceptionDescription errMsg;
297           errMsg << "Element not recognized :     322           errMsg << "Element not recognized : " << element << G4endl;
298                                                   323 
299           G4Exception("PDBlib::Load", "ELEM_NO << 324           G4Exception("PDBlib::Load",
300 #endif                                         << 325                       "ELEM_NOT_RECOGNIZED",
                                                   >> 326                       FatalException,
                                                   >> 327                       errMsg);
                                                   >> 328 #endif          
301         }                                         329         }
302                                                   330 
303         {                                         331         {
304           // resName :                         << 332           nbAtomTot++;
                                                   >> 333 
                                                   >> 334           //resName :
305           resName = sLine.substr(17, 3);          335           resName = sLine.substr(17, 3);
306           // resSeq :                          << 336           //resSeq :
307           istringstream((sLine.substr(22, 4)))    337           istringstream((sLine.substr(22, 4))) >> resSeq;
308           // x,y,z :                           << 338           //x,y,z :
309           stringstream((sLine.substr(30, 8)))     339           stringstream((sLine.substr(30, 8))) >> x;
310           stringstream((sLine.substr(38, 8)))     340           stringstream((sLine.substr(38, 8))) >> y;
311           stringstream((sLine.substr(46, 8)))     341           stringstream((sLine.substr(46, 8))) >> z;
312           // occupancy :                       << 342           //occupancy :
313           occupancy = 1.;                         343           occupancy = 1.;
314                                                   344 
315           if (minGlobZ < z) minGlobZ = z;      << 345           if(minGlobZ < z) minGlobZ = z;
316           if (maxGlobZ > z) maxGlobZ = z;      << 346           if(maxGlobZ > z) maxGlobZ = z;
317           if (minGlobX < x) minGlobX = x;      << 347           if(minGlobX < x) minGlobX = x;
318           if (maxGlobX > x) maxGlobX = x;      << 348           if(maxGlobX > x) maxGlobX = x;
319           if (minGlobY < y) minGlobY = y;      << 349           if(minGlobY < y) minGlobY = y;
320           if (maxGlobY > y) maxGlobY = y;      << 350           if(maxGlobY > y) maxGlobY = y;
321           if (minX > x) minX = x;              << 351           if(minX > x) minX = x;
322           if (maxX < x) maxX = x;              << 352           if(maxX < x) maxX = x;
323           if (minY > y) minY = y;              << 353           if(minY > y) minY = y;
324           if (maxY < y) maxY = y;              << 354           if(maxY < y) maxY = y;
325           if (minZ > z) minZ = z;              << 355           if(minZ > z) minZ = z;
326           if (maxZ < z) maxZ = z;              << 356           if(maxZ < z) maxZ = z;
327                                                << 357 
328           // treatment for Atoms:              << 358           //treatment for Atoms:
329           if (AtomicOld == NULL) {             << 359           if(AtomicOld == NULL)
330             AtomicOld =                        << 360           {
331               new Atom(serial, atomName, "", 0 << 361             AtomicOld = new Atom(serial,
332             AtomicOld->SetNext(NULL);  // If o << 362                                  atomName,
333           }                                    << 363                                  "",
334           else {                               << 364                                  0,
335             AtomicNext =                       << 365                                  0,
336               new Atom(serial, atomName, "", 0 << 366                                  x,
                                                   >> 367                                  y,
                                                   >> 368                                  z,
                                                   >> 369                                  vdwRadius,
                                                   >> 370                                  occupancy,
                                                   >> 371                                  0,
                                                   >> 372                                  element);
                                                   >> 373             AtomicOld->SetNext(NULL); //If only one Atom inside residue
                                                   >> 374           }
                                                   >> 375           else
                                                   >> 376           {
                                                   >> 377             AtomicNext = new Atom(serial,
                                                   >> 378                                   atomName,
                                                   >> 379                                   "",
                                                   >> 380                                   0,
                                                   >> 381                                   0,
                                                   >> 382                                   x,
                                                   >> 383                                   y,
                                                   >> 384                                   z,
                                                   >> 385                                   vdwRadius,
                                                   >> 386                                   occupancy,
                                                   >> 387                                   0,
                                                   >> 388                                   element);
337             AtomicOld->SetNext(AtomicNext);       389             AtomicOld->SetNext(AtomicNext);
338             AtomicOld = AtomicNext;               390             AtomicOld = AtomicNext;
339           }                                       391           }
340           nbAtom++;                               392           nbAtom++;
341         }  // END if (element!="H")            << 393         } //END if (element!="H")
342         /****************************Begin Res    394         /****************************Begin Residue************************/
343         // treatment for residues:             << 395         //treatment for residues:
344         if (residueOld == NULL) {              << 396         if(residueOld == NULL)
345           if (verbose > 2) G4cout << "residueO << 397         {
                                                   >> 398           if(verbose > 2) G4cout << "residueOld == NULL" << G4endl;
346                                                   399 
347           AtomicOld->fNumInRes = 0;               400           AtomicOld->fNumInRes = 0;
348           residueOld = new Residue(resName, re    401           residueOld = new Residue(resName, resSeq);
349           residueOld->SetFirst(AtomicOld);        402           residueOld->SetFirst(AtomicOld);
350           residueOld->SetNext(NULL);              403           residueOld->SetNext(NULL);
351           residueFirst = residueOld;              404           residueFirst = residueOld;
352           lastResSeq = resSeq;                    405           lastResSeq = resSeq;
353           nbResidue++;                            406           nbResidue++;
354         }                                         407         }
355         else {                                 << 408         else
356           if (lastResSeq == resSeq) {          << 409         {
                                                   >> 410           if(lastResSeq == resSeq)
                                                   >> 411           {
357             numAtomInRes++;                       412             numAtomInRes++;
358             AtomicOld->fNumInRes = numAtomInRe    413             AtomicOld->fNumInRes = numAtomInRes;
359           }                                       414           }
360           else {                               << 415           else
                                                   >> 416           {
361             numAtomInRes = 0;                     417             numAtomInRes = 0;
362             AtomicOld->fNumInRes = numAtomInRe    418             AtomicOld->fNumInRes = numAtomInRes;
363                                                   419 
364             residueNext = new Residue(resName,    420             residueNext = new Residue(resName, resSeq);
365             residueNext->SetFirst(AtomicOld);     421             residueNext->SetFirst(AtomicOld);
366             residueOld->SetNext(residueNext);     422             residueOld->SetNext(residueNext);
367             residueOld->fNbAtom = nbAtom - 1;     423             residueOld->fNbAtom = nbAtom - 1;
368                                                   424 
369             nbAtom = 1;                           425             nbAtom = 1;
370             residueOld = residueNext;             426             residueOld = residueNext;
371             lastResSeq = resSeq;                  427             lastResSeq = resSeq;
372             nbResidue++;                          428             nbResidue++;
373           }                                       429           }
374         }                                         430         }
375         /////////////////////////End Residue//    431         /////////////////////////End Residue////////////
376         ///////////// End ATOM ///////////////    432         ///////////// End ATOM ///////////////////
377       }  // END if Atom                        << 433       } //END if Atom
378     }  // END while (!infile.eof())            << 434     } //END while (!infile.eof())
379                                                   435 
380     infile.close();                               436     infile.close();
381     return moleculeFirst;                         437     return moleculeFirst;
382   }  // END else if (!infile)                  << 438   } //END else if (!infile)
383   return NULL;                                    439   return NULL;
384 }                                                 440 }
385                                                   441 
386 //....oooOO0OOooo........oooOO0OOooo........oo    442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
387                                                   443 
388 /**                                               444 /**
389  * \brief    Compute barycenters                  445  * \brief    Compute barycenters
390  * \details  Compute barycenters and its coord << 446  * \details  Compute barycenters and its coordinate 
391  *                  for nucleotides               447  *                  for nucleotides
392  * \param    moleculeListTemp *    moleculeLis    448  * \param    moleculeListTemp *    moleculeList
393  * \return   Barycenter *                         449  * \return   Barycenter *
394  */                                               450  */
395 Barycenter* PDBlib::ComputeNucleotideBarycente << 451 Barycenter * PDBlib::ComputeNucleotideBarycenters(Molecule * moleculeListTemp)
396 {                                                 452 {
397   ////////////////////////////////////////////    453   ///////////////////////////////////////////////////////
398   // Placement and physical volume constructio << 454   //Placement and physical volume construction from memory
399   Barycenter* BarycenterFirst = NULL;          << 455   Barycenter * BarycenterFirst = NULL;
400   Barycenter* BarycenterOld = NULL;            << 456   Barycenter * BarycenterOld = NULL;
401   Barycenter* BarycenterNext = NULL;           << 457   Barycenter * BarycenterNext = NULL;
402                                                << 458 
403   // Residue (Base, Phosphate,sugar) list      << 459   //Residue (Base, Phosphate,sugar) list
404   Residue* residueListTemp;                    << 460   Residue *residueListTemp;
405   // Atom list inside a residu                 << 461   //Atom list inside a residu
406   Atom* AtomTemp;                              << 462   Atom *AtomTemp;
407                                                   463 
408   int k = 0;                                      464   int k = 0;
409   int j_old = 0;                                  465   int j_old = 0;
410                                                   466 
411   while (moleculeListTemp) {                   << 467   while(moleculeListTemp)
                                                   >> 468   {
412     residueListTemp = moleculeListTemp->GetFir    469     residueListTemp = moleculeListTemp->GetFirst();
413                                                   470 
414     k++;                                          471     k++;
415     int j = 0;                                    472     int j = 0;
416                                                   473 
417     // Check numerotation style (1->n per stra << 474     //Check numerotation style (1->n per strand or 1->2n for two strand)
418     int correctNumerotationNumber = 0;            475     int correctNumerotationNumber = 0;
419     if (k == 2 && residueListTemp->fResSeq > 1 << 476     if(k == 2 && residueListTemp->fResSeq > 1)
                                                   >> 477     {
420       correctNumerotationNumber = residueListT    478       correctNumerotationNumber = residueListTemp->fResSeq;
421     }                                             479     }
422                                                   480 
423     while (residueListTemp) {                  << 481     while(residueListTemp)
                                                   >> 482     {
424       AtomTemp = residueListTemp->GetFirst();     483       AtomTemp = residueListTemp->GetFirst();
425       j++;                                        484       j++;
426                                                   485 
427       // Correction consequently to numerotati << 486       //Correction consequently to numerotation check
428       if (correctNumerotationNumber != 0) {    << 487       if(correctNumerotationNumber != 0)
429         residueListTemp->fResSeq = residueList << 488       {
                                                   >> 489         residueListTemp->fResSeq = residueListTemp->fResSeq
                                                   >> 490             - correctNumerotationNumber
                                                   >> 491                                    + 1;
430       }                                           492       }
431                                                   493 
432       // Barycenter computation                << 494       //Barycenter computation
433       double baryX = 0., baryY = 0., baryZ = 0    495       double baryX = 0., baryY = 0., baryZ = 0.;
434       double baryBaseX = 0., baryBaseY = 0., b    496       double baryBaseX = 0., baryBaseY = 0., baryBaseZ = 0.;
435       double barySugX = 0., barySugY = 0., bar    497       double barySugX = 0., barySugY = 0., barySugZ = 0.;
436       double baryPhosX = 0., baryPhosY = 0., b    498       double baryPhosX = 0., baryPhosY = 0., baryPhosZ = 0.;
437       unsigned short int nbAtomInBase = 0;        499       unsigned short int nbAtomInBase = 0;
438                                                   500 
439       for (int i = 0; i < residueListTemp->fNb << 501       for(int i = 0; i < residueListTemp->fNbAtom; i++)
440         // Compute barycenter of the nucletoti << 502       {
                                                   >> 503         //Compute barycenter of the nucletotide
441         baryX += AtomTemp->fX;                    504         baryX += AtomTemp->fX;
442         baryY += AtomTemp->fY;                    505         baryY += AtomTemp->fY;
443         baryZ += AtomTemp->fZ;                    506         baryZ += AtomTemp->fZ;
444         // Compute barycenters for Base Sugar  << 507         //Compute barycenters for Base Sugar Phosphat
445         if (residueListTemp->fResSeq == 1) {   << 508         if(residueListTemp->fResSeq == 1)
446           if (i == 0) {                        << 509         {
                                                   >> 510           if(i == 0)
                                                   >> 511           {
447             baryPhosX += AtomTemp->fX;            512             baryPhosX += AtomTemp->fX;
448             baryPhosY += AtomTemp->fY;            513             baryPhosY += AtomTemp->fY;
449             baryPhosZ += AtomTemp->fZ;            514             baryPhosZ += AtomTemp->fZ;
450           }                                       515           }
451           else if (i < 8) {                    << 516           else if(i < 8)
                                                   >> 517           {
452             barySugX += AtomTemp->fX;             518             barySugX += AtomTemp->fX;
453             barySugY += AtomTemp->fY;             519             barySugY += AtomTemp->fY;
454             barySugZ += AtomTemp->fZ;             520             barySugZ += AtomTemp->fZ;
455           }                                       521           }
456           else {                               << 522           else
457             // hydrogen are placed at the end  << 523           {
458             // We don't want them for this cal << 524             //hydrogen are placed at the end of the residue in a PDB file
459             if (AtomTemp->fElement != "H") {   << 525             //We don't want them for this calculation
                                                   >> 526             if(AtomTemp->fElement != "H")
                                                   >> 527             {
460               baryBaseX += AtomTemp->fX;          528               baryBaseX += AtomTemp->fX;
461               baryBaseY += AtomTemp->fY;          529               baryBaseY += AtomTemp->fY;
462               baryBaseZ += AtomTemp->fZ;          530               baryBaseZ += AtomTemp->fZ;
463               nbAtomInBase++;                     531               nbAtomInBase++;
464             }                                     532             }
465           }                                       533           }
466         }                                         534         }
467         else {                                 << 535         else
468           if (i < 4) {                         << 536         {
                                                   >> 537           if(i < 4)
                                                   >> 538           {
469             baryPhosX += AtomTemp->fX;            539             baryPhosX += AtomTemp->fX;
470             baryPhosY += AtomTemp->fY;            540             baryPhosY += AtomTemp->fY;
471             baryPhosZ += AtomTemp->fZ;            541             baryPhosZ += AtomTemp->fZ;
472           }                                       542           }
473           else if (i < 11) {                   << 543           else if(i < 11)
                                                   >> 544           {
474             barySugX += AtomTemp->fX;             545             barySugX += AtomTemp->fX;
475             barySugY += AtomTemp->fY;             546             barySugY += AtomTemp->fY;
476             barySugZ += AtomTemp->fZ;             547             barySugZ += AtomTemp->fZ;
477           }                                       548           }
478           else {                               << 549           else
479             // hydrogen are placed at the end  << 550           {
480             // We don't want them for this cal << 551             //hydrogen are placed at the end of the residue in a PDB file
481             if (AtomTemp->fElement != "H") {   << 552             //We don't want them for this calculation
                                                   >> 553             if(AtomTemp->fElement != "H")
                                                   >> 554             { // break;
482               baryBaseX += AtomTemp->fX;          555               baryBaseX += AtomTemp->fX;
483               baryBaseY += AtomTemp->fY;          556               baryBaseY += AtomTemp->fY;
484               baryBaseZ += AtomTemp->fZ;          557               baryBaseZ += AtomTemp->fZ;
485               nbAtomInBase++;                     558               nbAtomInBase++;
486             }                                     559             }
487           }                                       560           }
488         }                                         561         }
489         AtomTemp = AtomTemp->GetNext();           562         AtomTemp = AtomTemp->GetNext();
490       }  // end of for (  i=0 ; i < residueLis << 563       } //end of for (  i=0 ; i < residueListTemp->nbAtom ; i++)
491                                                   564 
492       baryX = baryX / (double)residueListTemp- << 565       baryX = baryX / (double) residueListTemp->fNbAtom;
493       baryY = baryY / (double)residueListTemp- << 566       baryY = baryY / (double) residueListTemp->fNbAtom;
494       baryZ = baryZ / (double)residueListTemp- << 567       baryZ = baryZ / (double) residueListTemp->fNbAtom;
495                                                   568 
496       if (residueListTemp->fResSeq != 1)  // S << 569       if(residueListTemp->fResSeq != 1) //Special case first Phosphate
497       {                                           570       {
498         baryPhosX = baryPhosX / 4.;               571         baryPhosX = baryPhosX / 4.;
499         baryPhosY = baryPhosY / 4.;               572         baryPhosY = baryPhosY / 4.;
500         baryPhosZ = baryPhosZ / 4.;               573         baryPhosZ = baryPhosZ / 4.;
501       }                                           574       }
502       barySugX = barySugX / 7.;                   575       barySugX = barySugX / 7.;
503       barySugY = barySugY / 7.;                   576       barySugY = barySugY / 7.;
504       barySugZ = barySugZ / 7.;                   577       barySugZ = barySugZ / 7.;
505       baryBaseX = baryBaseX / (double)nbAtomIn << 578       baryBaseX = baryBaseX / (double) nbAtomInBase;
506       baryBaseY = baryBaseY / (double)nbAtomIn << 579       baryBaseY = baryBaseY / (double) nbAtomInBase;
507       baryBaseZ = baryBaseZ / (double)nbAtomIn << 580       baryBaseZ = baryBaseZ / (double) nbAtomInBase;
508                                                << 581 
509       // Barycenter creation:                  << 582       //Barycenter creation:
510       if (BarycenterOld == NULL) {             << 583       if(BarycenterOld == NULL)
511         BarycenterOld = new Barycenter(j + j_o << 584       {
512                                        baryBas << 585         BarycenterOld = new Barycenter(j + j_old, baryX, baryY, baryZ, //j [1..n]
513                                        barySug << 586                                        baryBaseX,
                                                   >> 587                                        baryBaseY,
                                                   >> 588                                        baryBaseZ,
                                                   >> 589                                        barySugX,
                                                   >> 590                                        barySugY,
                                                   >> 591                                        barySugZ,
                                                   >> 592                                        baryPhosX,
                                                   >> 593                                        baryPhosY,
                                                   >> 594                                        baryPhosZ);
514         BarycenterFirst = BarycenterOld;          595         BarycenterFirst = BarycenterOld;
515       }                                           596       }
516       else {                                   << 597       else
517         BarycenterNext =                       << 598       {
518           new Barycenter(j + j_old, baryX, bar << 599         BarycenterNext = new Barycenter(j + j_old,
519                          barySugY, barySugZ, b << 600                                         baryX,
                                                   >> 601                                         baryY,
                                                   >> 602                                         baryZ,
                                                   >> 603                                         baryBaseX,
                                                   >> 604                                         baryBaseY,
                                                   >> 605                                         baryBaseZ,
                                                   >> 606                                         barySugX,
                                                   >> 607                                         barySugY,
                                                   >> 608                                         barySugZ,
                                                   >> 609                                         baryPhosX,
                                                   >> 610                                         baryPhosY,
                                                   >> 611                                         baryPhosZ);
520         BarycenterOld->SetNext(BarycenterNext)    612         BarycenterOld->SetNext(BarycenterNext);
521         BarycenterOld = BarycenterNext;           613         BarycenterOld = BarycenterNext;
522       }                                           614       }
523                                                   615 
524       ////////////////////////////////////////    616       /////////////////////////////////////////////////
525       // distance computation between all atom << 617       //distance computation between all atoms inside
526       // a residue and the barycenter          << 618       //a residue and the barycenter
527       AtomTemp = residueListTemp->GetFirst();     619       AtomTemp = residueListTemp->GetFirst();
528       double dT3Dp;                               620       double dT3Dp;
529       double max = 0.;                            621       double max = 0.;
530       for (int ii = 0; ii < residueListTemp->f << 622       for(int ii = 0; ii < residueListTemp->fNbAtom; ii++)
531         dT3Dp = DistanceTwo3Dpoints(AtomTemp-> << 623       {
532                                     Barycenter << 624         dT3Dp = DistanceTwo3Dpoints(AtomTemp->fX,
                                                   >> 625                                     BarycenterOld->fCenterX,
                                                   >> 626                                     AtomTemp->fY,
                                                   >> 627                                     BarycenterOld->fCenterY,
                                                   >> 628                                     AtomTemp->fZ,
                                                   >> 629                                     BarycenterOld->fCenterZ);
533         BarycenterOld->SetDistance(ii, dT3Dp);    630         BarycenterOld->SetDistance(ii, dT3Dp);
534         if (dT3Dp > max) max = dT3Dp;          << 631         if(dT3Dp > max) max = dT3Dp;
535         AtomTemp = AtomTemp->GetNext();           632         AtomTemp = AtomTemp->GetNext();
536       }  // end of for (  i=0 ; i < residueLis << 633       } //end of for (  i=0 ; i < residueListTemp->nbAtom ; i++)
537                                                   634 
538       BarycenterOld->SetRadius(max + 1.8);        635       BarycenterOld->SetRadius(max + 1.8);
539       residueListTemp = residueListTemp->GetNe    636       residueListTemp = residueListTemp->GetNext();
540                                                   637 
541     }  // end of while sur residueListTemp     << 638     } //end of while sur residueListTemp
542                                                   639 
543     j_old += j;                                   640     j_old += j;
544                                                   641 
545     /// molecs->push_back(*moleculeListTemp);  << 642     ///molecs->push_back(*moleculeListTemp);
546     moleculeListTemp = moleculeListTemp->GetNe    643     moleculeListTemp = moleculeListTemp->GetNext();
547   }  // end of while sur moleculeListTemp      << 644   } //end of while sur moleculeListTemp
548                                                   645 
549   if (BarycenterNext != NULL) {                << 646   if(BarycenterNext != NULL)
                                                   >> 647   {
550     BarycenterNext->SetNext(NULL);                648     BarycenterNext->SetNext(NULL);
551   }                                               649   }
552                                                   650 
553   return BarycenterFirst;                         651   return BarycenterFirst;
554 }                                                 652 }
555                                                   653 
556 //....oooOO0OOooo........oooOO0OOooo........oo    654 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
557                                                   655 
558 /**                                               656 /**
559  * \brief    the corresponding bounding volume    657  * \brief    the corresponding bounding volume parameters
560  * \details  the corresponding bounding volume    658  * \details  the corresponding bounding volume parameters
561  *            to build a box from atoms coordi    659  *            to build a box from atoms coordinates
562  */                                               660  */
563 void PDBlib::ComputeBoundingVolumeParams(Molec << 661 void PDBlib::ComputeBoundingVolumeParams(Molecule *moleculeListTemp,
564                                          doubl << 662                                          double &dX,
565                                          doubl << 663                                          double &dY,
566                                          doubl << 664                                          double &dZ, //Dimensions for bounding volume
                                                   >> 665                                          double &tX,
                                                   >> 666                                          double &tY,
                                                   >> 667                                          double &tZ) //Translation for bounding volume
567 {                                                 668 {
568   double minminX, minminY, minminZ;  // minimu << 669   double minminX, minminY, minminZ; //minimum minimorum
569   double maxmaxX, maxmaxY, maxmaxZ;  // maximu << 670   double maxmaxX, maxmaxY, maxmaxZ; //maximum maximorum
570                                                   671 
571   minminX = DBL_MAX;                              672   minminX = DBL_MAX;
572   minminY = DBL_MAX;                              673   minminY = DBL_MAX;
573   minminZ = DBL_MAX;                              674   minminZ = DBL_MAX;
574   maxmaxX = -DBL_MAX;                             675   maxmaxX = -DBL_MAX;
575   maxmaxY = -DBL_MAX;                             676   maxmaxY = -DBL_MAX;
576   maxmaxZ = -DBL_MAX;                             677   maxmaxZ = -DBL_MAX;
577                                                   678 
578   while (moleculeListTemp) {                   << 679   while(moleculeListTemp)
579     if (minminX > moleculeListTemp->fMaxGlobX) << 680   {
                                                   >> 681     if(minminX > moleculeListTemp->fMaxGlobX)
                                                   >> 682     {
580       minminX = moleculeListTemp->fMaxGlobX;      683       minminX = moleculeListTemp->fMaxGlobX;
581     }                                             684     }
582     if (minminY > moleculeListTemp->fMaxGlobY) << 685     if(minminY > moleculeListTemp->fMaxGlobY)
                                                   >> 686     {
583       minminY = moleculeListTemp->fMaxGlobY;      687       minminY = moleculeListTemp->fMaxGlobY;
584     }                                             688     }
585     if (minminZ > moleculeListTemp->fMaxGlobZ) << 689     if(minminZ > moleculeListTemp->fMaxGlobZ)
                                                   >> 690     {
586       minminZ = moleculeListTemp->fMaxGlobZ;      691       minminZ = moleculeListTemp->fMaxGlobZ;
587     }                                             692     }
588     if (maxmaxX < moleculeListTemp->fMinGlobX) << 693     if(maxmaxX < moleculeListTemp->fMinGlobX)
                                                   >> 694     {
589       maxmaxX = moleculeListTemp->fMinGlobX;      695       maxmaxX = moleculeListTemp->fMinGlobX;
590     }                                             696     }
591     if (maxmaxY < moleculeListTemp->fMinGlobY) << 697     if(maxmaxY < moleculeListTemp->fMinGlobY)
                                                   >> 698     {
592       maxmaxY = moleculeListTemp->fMinGlobY;      699       maxmaxY = moleculeListTemp->fMinGlobY;
593     }                                             700     }
594     if (maxmaxZ < moleculeListTemp->fMinGlobZ) << 701     if(maxmaxZ < moleculeListTemp->fMinGlobZ)
                                                   >> 702     {
595       maxmaxZ = moleculeListTemp->fMinGlobZ;      703       maxmaxZ = moleculeListTemp->fMinGlobZ;
596     }                                             704     }
597                                                   705 
598     moleculeListTemp = moleculeListTemp->GetNe    706     moleculeListTemp = moleculeListTemp->GetNext();
599   }  // end of while sur moleculeListTemp      << 707   } //end of while sur moleculeListTemp
600                                                   708 
601   dX = (maxmaxX - minminX) / 2. + 1.8;  // 1.8 << 709   dX = (maxmaxX - minminX) / 2. + 1.8; //1.8 => size of biggest radius for atoms
602   dY = (maxmaxY - minminY) / 2. + 1.8;            710   dY = (maxmaxY - minminY) / 2. + 1.8;
603   dZ = (maxmaxZ - minminZ) / 2. + 1.8;            711   dZ = (maxmaxZ - minminZ) / 2. + 1.8;
604                                                   712 
605   tX = minminX + (maxmaxX - minminX) / 2.;        713   tX = minminX + (maxmaxX - minminX) / 2.;
606   tY = minminY + (maxmaxY - minminY) / 2.;        714   tY = minminY + (maxmaxY - minminY) / 2.;
607   tZ = minminZ + (maxmaxZ - minminZ) / 2.;        715   tZ = minminZ + (maxmaxZ - minminZ) / 2.;
608 }                                                 716 }
609                                                   717 
610 //....oooOO0OOooo........oooOO0OOooo........oo    718 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
611                                                   719 
612 /**                                               720 /**
613  * \brief    Compute number of nucleotide per     721  * \brief    Compute number of nucleotide per strand
614  * \details  Compute number of nucleotide per  << 722  * \details  Compute number of nucleotide per strand 
615  *                  for DNA                       723  *                  for DNA
616  */                                               724  */
617 void PDBlib::ComputeNbNucleotidsPerStrand(Mole << 725 void PDBlib::ComputeNbNucleotidsPerStrand(Molecule * moleculeListTemp)
618 {                                                 726 {
619   Residue* residueListTemp;                    << 727   Residue *residueListTemp;
620                                                   728 
                                                   >> 729   int k = 0;
621   int j_old = 0;                                  730   int j_old = 0;
622                                                   731 
623   while (moleculeListTemp) {                   << 732   while(moleculeListTemp)
                                                   >> 733   {
624     residueListTemp = moleculeListTemp->GetFir    734     residueListTemp = moleculeListTemp->GetFirst();
625                                                   735 
                                                   >> 736     k++;
626     int j = 0;                                    737     int j = 0;
627                                                   738 
628     while (residueListTemp) {                  << 739     while(residueListTemp)
                                                   >> 740     {
629       j++;                                        741       j++;
630       residueListTemp = residueListTemp->GetNe    742       residueListTemp = residueListTemp->GetNext();
631     }  // end of while sur residueListTemp     << 743     } //end of while sur residueListTemp
632                                                   744 
633     j_old += j;                                   745     j_old += j;
634     moleculeListTemp = moleculeListTemp->GetNe    746     moleculeListTemp = moleculeListTemp->GetNext();
635   }  // end of while sur moleculeListTemp      << 747   } //end of while sur moleculeListTemp
636                                                   748 
637   fNbNucleotidsPerStrand = j_old / 2;             749   fNbNucleotidsPerStrand = j_old / 2;
638 }                                                 750 }
639                                                   751 
640 //....oooOO0OOooo........oooOO0OOooo........oo    752 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
641                                                   753 
642 /**                                               754 /**
643  * \brief    Compute barycenters                  755  * \brief    Compute barycenters
644  * \details  Compute barycenters and its coord << 756  * \details  Compute barycenters and its coordinate 
645  *                  for nucleotides               757  *                  for nucleotides
646  */                                               758  */
647 unsigned short int PDBlib::ComputeMatchEdepDNA << 759 unsigned short int PDBlib::ComputeMatchEdepDNA(Barycenter *BarycenterList,
648                                                << 760                                                Molecule *moleculeListTemp,
649                                                << 761                                                double x,
650                                                << 762                                                double y,
                                                   >> 763                                                double z,
                                                   >> 764                                                int &numStrand,
                                                   >> 765                                                int &numNucleotid,
                                                   >> 766                                                int &codeResidue)
651 {                                                 767 {
652   unsigned short int matchFound = 0;              768   unsigned short int matchFound = 0;
653                                                   769 
654   Molecule* mLTsavedPointer = moleculeListTemp << 770   Molecule *mLTsavedPointer = moleculeListTemp;
655   Barycenter* BLsavedPointer = BarycenterList; << 771   Barycenter *BLsavedPointer = BarycenterList;
656                                                   772 
657   short int strandNum = 0;  // Strand number   << 773   short int strandNum = 0; //Strand number
658   int residueNum = 1;  // Residue (nucleotide) << 774   int residueNum = 1; //Residue (nucleotide) number
659   G4String baseName;  // Base name [A,C,T,G]   << 775   G4String baseName; //Base name [A,C,T,G]
660   unsigned short int BSP = 2;  // Base (defaul << 776   unsigned short int BSP = 2; //Base (default value), Sugar, Phosphat
661                                                   777 
662   double smallestDist;  // smallest dist Atom  << 778   double smallestDist; //smallest dist Atom <-> edep coordinates
663   double distEdepDNA;                             779   double distEdepDNA;
664   double distEdepAtom;                            780   double distEdepAtom;
665                                                   781 
666   // Residue (Base, Phosphate,suggar) list     << 782   //Residue (Base, Phosphate,suggar) list
667   Residue* residueListTemp;                    << 783   Residue *residueListTemp;
668   // Atom list inside a residue                << 784   //Atom list inside a residue
669   Atom* AtomTemp;                              << 785   Atom *AtomTemp;
670                                                   786 
671   int k = 0;  // Molecule number               << 787   int k = 0; //Molecule number
672   moleculeListTemp = mLTsavedPointer;             788   moleculeListTemp = mLTsavedPointer;
673   BarycenterList = BLsavedPointer;                789   BarycenterList = BLsavedPointer;
674                                                   790 
675   smallestDist = 33.0;  // Sufficiently large  << 791   smallestDist = 33.0; //Sufficiently large value
676   while (moleculeListTemp) {                   << 792   while(moleculeListTemp)
                                                   >> 793   {
677     k++;                                          794     k++;
678     residueListTemp = moleculeListTemp->GetFir    795     residueListTemp = moleculeListTemp->GetFirst();
679                                                   796 
680     int j = 0;  // Residue number              << 797     int j = 0; //Residue number
681                                                   798 
682     int j_save = INT_MAX;  // Saved res. numbe << 799     int j_save = INT_MAX; //Saved res. number if match found
683                                                   800 
684     while (residueListTemp) {                  << 801     while(residueListTemp)
                                                   >> 802     {
685       j++;                                        803       j++;
686                                                   804 
687       if (j - j_save > 2) break;               << 805       if(j - j_save > 2) break;
688                                                   806 
689       distEdepDNA = DistanceTwo3Dpoints(x, Bar << 807       distEdepDNA = DistanceTwo3Dpoints(x,
                                                   >> 808                                         BarycenterList->fCenterX,
                                                   >> 809                                         y,
                                                   >> 810                                         BarycenterList->fCenterY,
                                                   >> 811                                         z,
690                                         Baryce    812                                         BarycenterList->fCenterZ);
691       if (distEdepDNA < BarycenterList->GetRad << 813       if(distEdepDNA < BarycenterList->GetRadius())
692         // Find the closest atom               << 814       {
693         // Compute distance between energy dep << 815         //Find the closest atom
694         // if distance <1.8 then match OK but  << 816         //Compute distance between energy deposited and atoms for a residue
                                                   >> 817         //if distance <1.8 then match OK but search inside 2 next residues
695         AtomTemp = residueListTemp->GetFirst()    818         AtomTemp = residueListTemp->GetFirst();
696         for (int iii = 0; iii < residueListTem << 819         for(int iii = 0; iii < residueListTemp->fNbAtom; iii++)
697           distEdepAtom =                       << 820         {
698             DistanceTwo3Dpoints(x, AtomTemp->G << 
699                                                   821 
700           if ((distEdepAtom < AtomTemp->GetVan << 822           distEdepAtom = DistanceTwo3Dpoints(x,
                                                   >> 823                                              AtomTemp->GetX(),
                                                   >> 824                                              y,
                                                   >> 825                                              AtomTemp->GetY(),
                                                   >> 826                                              z,
                                                   >> 827                                              AtomTemp->GetZ());
                                                   >> 828 
                                                   >> 829           if((distEdepAtom < AtomTemp->GetVanDerWaalsRadius()) && (smallestDist
                                                   >> 830               > distEdepAtom))
                                                   >> 831           {
701             strandNum = k;                        832             strandNum = k;
702                                                   833 
703             if (k == 2) {                      << 834             if(k == 2)
704               residueNum = fNbNucleotidsPerStr << 835             {
                                                   >> 836               residueNum = fNbNucleotidsPerStrand + 1
                                                   >> 837                   - residueListTemp->fResSeq;
705             }                                     838             }
706             else {                             << 839             else
                                                   >> 840             {
707               residueNum = residueListTemp->fR    841               residueNum = residueListTemp->fResSeq;
708             }                                     842             }
709                                                   843 
710             baseName = (residueListTemp->fResN    844             baseName = (residueListTemp->fResName)[2];
711             if (residueListTemp->fResSeq == 1) << 845             if(residueListTemp->fResSeq == 1)
712               if (iii == 0)                    << 846             {
713                 BSP = 0;  //"Phosphate"        << 847               if(iii == 0) BSP = 0; //"Phosphate"
714               else if (iii < 8)                << 848               else if(iii < 8) BSP = 1; //"Sugar"
715                 BSP = 1;  //"Sugar"            << 849               else BSP = 2; //"Base"
716               else                             << 
717                 BSP = 2;  //"Base"             << 
718             }                                     850             }
719             else {                             << 851             else
720               if (iii < 4)                     << 852             {
721                 BSP = 0;  //"Phosphate"        << 853               if(iii < 4) BSP = 0; //"Phosphate"
722               else if (iii < 11)               << 854               else if(iii < 11) BSP = 1; //"Sugar"
723                 BSP = 1;  //"Sugar"            << 855               else BSP = 2; //"Base"
724               else                             << 
725                 BSP = 2;  //"Base"             << 
726             }                                     856             }
727                                                   857 
728             smallestDist = distEdepAtom;          858             smallestDist = distEdepAtom;
729                                                   859 
730             int j_max_value = INT_MAX;            860             int j_max_value = INT_MAX;
731                                                   861 
732             if (j_save == j_max_value) j_save  << 862             if(j_save == j_max_value) j_save = j;
733             matchFound = 1;                       863             matchFound = 1;
734           }                                       864           }
735           AtomTemp = AtomTemp->GetNext();         865           AtomTemp = AtomTemp->GetNext();
736         }  // end of for (  iii=0 ; iii < resi << 866         } //end of for (  iii=0 ; iii < residueListTemp->nbAtom ; iii++)
737       }  // end for if (distEdepDNA < Barycent << 867       } //end for if (distEdepDNA < BarycenterList->GetRadius())
738       BarycenterList = BarycenterList->GetNext    868       BarycenterList = BarycenterList->GetNext();
739       residueListTemp = residueListTemp->GetNe    869       residueListTemp = residueListTemp->GetNext();
740     }  // end of while sur residueListTemp     << 870     } //end of while sur residueListTemp
741     moleculeListTemp = moleculeListTemp->GetNe    871     moleculeListTemp = moleculeListTemp->GetNext();
742   }  // end of while sur moleculeListTemp      << 872   } //end of while sur moleculeListTemp
743                                                   873 
744   numStrand = strandNum;                          874   numStrand = strandNum;
745   numNucleotid = residueNum;                      875   numNucleotid = residueNum;
746   codeResidue = BSP;                              876   codeResidue = BSP;
747                                                   877 
748   return matchFound;                              878   return matchFound;
749 }                                                 879 }
750                                                   880 
751 //....oooOO0OOooo........oooOO0OOooo........oo    881 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
752                                                   882 
753 double PDBlib::DistanceTwo3Dpoints(double xA,  << 883 double PDBlib::DistanceTwo3Dpoints(double xA,
                                                   >> 884                                    double xB,
                                                   >> 885                                    double yA,
                                                   >> 886                                    double yB,
                                                   >> 887                                    double zA,
                                                   >> 888                                    double zB)
754 {                                                 889 {
755   return sqrt((xA - xB) * (xA - xB) + (yA - yB << 890   return sqrt((xA - xB) * (xA - xB) + (yA - yB) * (yA - yB)
                                                   >> 891               + (zA - zB) * (zA - zB));
756 }                                                 892 }
757                                                   893