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