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