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