ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DNAParser.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DNAParser.cc
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 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Updated: H. Tran (IRSN), France: 20/12/2018
28 //
29 
30 #include "DNAParser.hh"
31 #include <fstream>
32 #include "G4SystemOfUnits.hh"
33 #include "G4DNAChemistryManager.hh"
34 #include "G4VPhysicalVolume.hh"
35 #include "G4LogicalVolume.hh"
36 #include "G4VPhysicalVolume.hh"
37 #include "G4VSolid.hh"
38 #include "G4PVPlacement.hh"
39 #include "G4NistManager.hh"
40 #include "G4Material.hh"
41 #include "G4Box.hh"
42 #include "G4Orb.hh"
43 #include "G4SubtractionSolid.hh"
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
46 // Molecule struct def
47 struct DNAParser::Molecule
48 {
49  Molecule(std::string name,
50  int copyNumber,
51  const G4ThreeVector& position,
52  double radius,
53  double waterRadius,
54  const std::string& material,
55  int strand)
56  : fName(name),
57  fMaterial(material),
58  fCopyNumber(copyNumber),
59  fPosition(position),
60  fRadius(radius),
61  fRadiusWater(waterRadius),
62  fStrand(strand)
63  {}
64 
65  G4String fName;
66  G4String fMaterial;
67  G4int fCopyNumber;
68  G4ThreeVector fPosition;
69  G4double fRadius;
70  G4double fRadiusWater;
71  G4int fStrand;
72 
73  G4bool operator<(const Molecule& rhs) const
74  {
75  return (fPosition.z() < rhs.fPosition.z());
76  }
77 };
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 DNAParser::DNAParser()
81  : fSize(0)
82  , fGeoName("")
83  , fpWater(nullptr)
84  , fpGun(new G4MoleculeGun())
85 {
86  EnumParser();
87 }
88 
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 
91 DNAParser::~DNAParser() = default;
92 
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
94 
95 G4LogicalVolume* DNAParser::CreateLogicVolume()
96 {
97  G4NistManager * pMan = G4NistManager::Instance();
98  fpWater = pMan->FindOrBuildMaterial("G4_WATER");
99 
100  G4String boxNameSolid = fGeoName + "_solid";
101 
102  G4Box* pBoxSolid = new G4Box(boxNameSolid,
103  fSize/2,
104  fSize/2,
105  fSize/2);
106 
107  G4String boxNameLogic = fGeoName + "_logic";
108 
109  G4LogicalVolume* pBoxLogic = new G4LogicalVolume(pBoxSolid,
110  fpWater,
111  boxNameLogic);
112 
113  std::sort( fMolecules.begin(), fMolecules.end() );
114 
115  for(int i = 0, ie = fMolecules.size(); i < ie; ++i)
116  {
117  G4String name = fMolecules[i].fName;
118  G4double radius = fMolecules[i].fRadius;
119  G4double waterRadius = fMolecules[i].fRadiusWater;
120  G4ThreeVector moleculePosition = fMolecules[i].fPosition;
121 
122  int copyNum = fMolecules[i].fCopyNumber;
123 
124  G4Orb* pMoleculeWaterSolid = nullptr;
125  G4VSolid* pMoleculeWaterCutSolid = nullptr;
126  G4LogicalVolume* pMoleculeWaterLogic = nullptr;
127  G4VPhysicalVolume* pPhysicalVolumeWater = nullptr;
128 
129  G4double tolerence = 1e-4 * nm;
130 
131  if(waterRadius > tolerence)
132  {
133  G4String nameWaterSolid = name + "_water_solid";
134  pMoleculeWaterSolid = new G4Orb(nameWaterSolid, waterRadius);
135  pMoleculeWaterCutSolid = CreateCutSolid(pMoleculeWaterSolid,
136  fMolecules[i],
137  fMolecules,
138  false);
139 
140  G4String nameWaterLogic = name + "_water_logic";
141 
142  pMoleculeWaterLogic = new G4LogicalVolume(pMoleculeWaterCutSolid,
143  fpWater,
144  nameWaterLogic);
145  G4String nameWaterPhys = name + "_water";
146  pPhysicalVolumeWater = new G4PVPlacement(0,
147  moleculePosition,
148  pMoleculeWaterLogic,
149  nameWaterPhys,
150  pBoxLogic,
151  false,
152  copyNum);
153 
154  auto it = fEnumMap.find(nameWaterPhys);
155  if(it != fEnumMap.end())
156  {
157  fGeometryMap.insert(std::make_pair(pPhysicalVolumeWater,
158  it->second));
159  }
160  else
161  {
162  G4ExceptionDescription exceptionDescription;
163  exceptionDescription <<nameWaterPhys
164  <<" could not be exsited";
165  G4Exception("DNAParser::CreateLogicVolume()",
166  "DNAParser001", FatalException,
167  exceptionDescription);
168  }
169  }
170 // Dna volume part
171 
172  G4Orb* pMoleculeSolid = nullptr;
173  G4VSolid* pMoleculeCutSolid = nullptr;
174  G4LogicalVolume* pMoleculeLogic = nullptr;
175  G4VPhysicalVolume* pPhysicalVolumeMolecule = nullptr;
176 
177  G4String nameSolid = fMolecules[i].fName + "_solid";
178  pMoleculeSolid = new G4Orb(nameSolid, radius);
179 
180  pMoleculeCutSolid = CreateCutSolid(pMoleculeSolid,
181  fMolecules[i],
182  fMolecules,
183  true);
184 
185  G4String nameLogic = name + "_logic";
186 
187  pMoleculeLogic = new G4LogicalVolume(pMoleculeCutSolid,
188  fpWater,
189  nameLogic);
190  if(waterRadius > tolerence)
191  {
192  G4ThreeVector position(0,0,0);
193  std::string namePhys = name;
194  pPhysicalVolumeMolecule = new G4PVPlacement(0,
195  position,
196  pMoleculeLogic,
197  namePhys,
198  pMoleculeWaterLogic,
199  false,
200  copyNum);
201 
202 
203  auto it = fEnumMap.find(namePhys);
204  if(it != fEnumMap.end())
205  {
206  fGeometryMap.insert(std::make_pair(pPhysicalVolumeMolecule,
207  it->second));
208  }
209  else
210  {
211  G4ExceptionDescription exceptionDescription;
212  exceptionDescription <<namePhys
213  <<" could not be exsited";
214  G4Exception("DNAParser::CreateLogicVolume()",
215  "DNAParser002", FatalException,
216  exceptionDescription);
217  }
218  }
219  else
220  {
221  G4ThreeVector position = moleculePosition;
222  G4String namePhys = name;
223  pPhysicalVolumeMolecule = new G4PVPlacement(0,
224  position,
225  pMoleculeLogic,
226  namePhys,
227  pBoxLogic,
228  false,
229  copyNum);
230 
231  auto it = fEnumMap.find(namePhys);
232 
233  if(it != fEnumMap.end())
234  {
235  fGeometryMap.insert(std::make_pair(pPhysicalVolumeMolecule,
236  it->second));
237  }
238  else
239  {
240  G4ExceptionDescription exceptionDescription;
241  exceptionDescription <<namePhys
242  <<" could not be exsited";
243  G4Exception("DNAParser::CreateLogicVolume()",
244  "DNAParser003", FatalException,
245  exceptionDescription);
246  }
247  }
248  }
249  fMolecules.clear();
250  fRadiusMap.clear();
251  fWaterRadiusMap.clear();
252  return pBoxLogic;
253 }
254 
255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
256 
257 void DNAParser::ParseFile(const std::string& fileName)
258 {
259  fMolecules.clear();
260  fRadiusMap.clear();
261  fWaterRadiusMap.clear();
262 
263  std::ifstream file(fileName.c_str());
264 
265  if(!file.is_open())
266  {
267  G4ExceptionDescription exceptionDescription;
268  exceptionDescription <<fileName
269  <<"could not be opened";
270  G4Exception("DNAParser::ParseFile()",
271  "DNAParser002", FatalException,
272  exceptionDescription);
273  }
274  std::string line;
275  while(std::getline(file, line) )
276  {
277  if( line.empty() ) continue;
278  std::istringstream issLine(line);
279  std::string firstItem;
280  issLine >> firstItem;
281 
282  if(firstItem == "_Name")
283  {
284  G4String name = "";
285  issLine >> name;
286  fGeoName = name;
287  }
288  if(firstItem == "_Size")
289  {
290  G4double size;
291  issLine >> size;
292  size *= nm;
293  fSize = size;
294  }
295  if(firstItem == "_Radius")
296  {
297  G4String name;
298  issLine >> name;
299  G4double radius;
300  issLine >> radius;
301  radius *= nm;
302  G4double waterRadius;
303  issLine >> waterRadius;
304  waterRadius *= nm;
305  fRadiusMap[name] = radius;
306  fWaterRadiusMap[name] = waterRadius;
307  }
308  if(firstItem == "_pl")
309  {
310  std::string name;
311  issLine >> name;
312  G4String material;
313  issLine >> material;
314 
315  G4int strand;
316  issLine >> strand;
317 
318  G4int copyNumber;
319  issLine >> copyNumber;
320 
321  G4double x;
322  issLine >> x;
323  x *= nm;
324 
325  G4double y;
326  issLine >> y;
327  y *= nm;
328 
329  G4double z;
330  issLine >> z;
331  z *= nm;
332 
333  Molecule molecule(name,
334  copyNumber,
335  G4ThreeVector(x, y, z),
336  fRadiusMap[name],
337  fWaterRadiusMap[name],
338  material,
339  strand);
340  fMolecules.push_back(molecule);
341 
342  auto itt = fEnumMap.find(name);
343 
344  if(itt != fEnumMap.end())
345  {
346 
347  if(itt->second != DNAVolumeType::phosphate1 &&
348  itt->second != DNAVolumeType::phosphate2)
349  {
350  if(itt->second == DNAVolumeType::deoxyribose1 ||
351  itt->second == DNAVolumeType::deoxyribose2)
352  {
353  name = "Deoxyribose";
354  }
355  if(itt->second == DNAVolumeType::base_adenine)
356  {
357  name = "Adenine";
358  }
359  if(itt->second == DNAVolumeType::base_thymine)
360  {
361  name = "Thymine";
362  }
363  if(itt->second == DNAVolumeType::base_cytosine)
364  {
365  name = "Cytosine";
366  }
367  if(itt->second == DNAVolumeType::base_guanine)
368  {
369  name = "Guanine";
370  }
371  if(itt->second == DNAVolumeType::histone)
372  {
373  name = "Histone";
374  }
375 
376  fpGun->AddMolecule(name,
377  G4ThreeVector(x, y, z),
378  1*picosecond);
379  }
380  }
381  else
382  {
383  G4ExceptionDescription exceptionDescription;
384  exceptionDescription <<name
385  <<" could not be exsited";
386  G4Exception("DNAParser::ParseFile()",
387  "DNAParser005", FatalException,
388  exceptionDescription);
389  }
390  }
391  }
392  file.close();
393 }
394 
395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
396 
397 G4VSolid* DNAParser::CreateCutSolid(G4Orb *pSolidOrbRef,
398  Molecule &molRef,
399  std::vector<Molecule> &molList,
400  G4bool in)
401 {
402 // The idea behing this method is to cut overlap volumes by selecting
403 //one of them (the reference) and checking all the other volumes (the targets).
404 // If a reference and a target volumes are close enough to overlap they will be cut.
405 // The reference is already selected when we enter this method.
406 // Use the tiny space to differentiate the frontiers (may not be necessary)
407  G4double tinySpace = 0.001*nm;
408 
409  G4SubtractionSolid* pSolidCut(NULL);
410  G4bool isCutted = false;
411  G4bool isOurVol = false;
412  G4double radiusRef;
413 
414  if(molRef.fRadiusWater == 0)
415  {
416  radiusRef = molRef.fRadius;
417  }
418  else
419  {
420  radiusRef = molRef.fRadiusWater;
421  }
422 
423  G4ThreeVector posRef = molRef.fPosition;
424 
425  if(std::abs(posRef.x() ) + radiusRef > fSize/2 // along x
426  || std::abs(posRef.y() ) + radiusRef > fSize/2 // along y
427  || std::abs(posRef.z() ) + radiusRef > fSize/2) // along z
428  {
429  G4Box* solidBox = new G4Box("solid_box_for_cut",
430  fSize/2,
431  fSize/2,
432  fSize/2);
433  G4ThreeVector posBox;
434  G4RotationMatrix rotMat;
435  rotMat.rotate(0, G4ThreeVector(0,0,1) );
436 
437  if(std::abs( posRef.y() + radiusRef ) > fSize/2 )
438  {
439  posBox = -posRef + G4ThreeVector(0,fSize,0);
440  if(!isCutted)
441  {
442  pSolidCut =
443  new G4SubtractionSolid("solidCut",
444  pSolidOrbRef,
445  solidBox,
446  &rotMat,
447  posBox);
448  isCutted = true;
449  }
450  else
451  {
452  pSolidCut =
453  new G4SubtractionSolid("solidCut",
454  pSolidCut,
455  solidBox,
456  &rotMat,
457  posBox);
458  }
459  }
460 // Down
461  if(std::abs( posRef.y() - radiusRef ) > fSize/2 )
462  {
463  posBox = -posRef + G4ThreeVector(0,-fSize,0);
464 
465  if(!isCutted)
466  {
467  pSolidCut = new G4SubtractionSolid("solidCut",
468  pSolidOrbRef,
469  solidBox,
470  &rotMat,
471  posBox);
472  isCutted = true;
473  }
474  else
475  {
476  pSolidCut = new G4SubtractionSolid("solidCut",
477  pSolidCut,
478  solidBox,
479  &rotMat,
480  posBox);
481  }
482  }
483 // Left
484  if(std::abs( posRef.x() + radiusRef ) > fSize/2 )
485  {
486  posBox = -posRef + G4ThreeVector(fSize,0,0);
487  if(!isCutted)
488  {
489  pSolidCut = new G4SubtractionSolid("solidCut",
490  pSolidOrbRef,
491  solidBox,
492  &rotMat,
493  posBox);
494  isCutted = true;
495  }
496  else
497  {
498  pSolidCut = new G4SubtractionSolid("solidCut",
499  pSolidCut,
500  solidBox,
501  &rotMat,
502  posBox);
503  }
504  }
505 
506 // Right
507  if(std::abs( posRef.x() - radiusRef ) > fSize/2 )
508  {
509  posBox = -posRef + G4ThreeVector(-fSize,0,0);
510  if(!isCutted)
511  {
512  pSolidCut = new G4SubtractionSolid("solidCut",
513  pSolidOrbRef,
514  solidBox,
515  &rotMat,
516  posBox);
517  isCutted = true;
518  }
519  else
520  {
521  pSolidCut = new G4SubtractionSolid("solidCut",
522  pSolidCut,
523  solidBox,
524  &rotMat,
525  posBox);
526  }
527  }
528 // Forward
529  if(std::abs( posRef.z() + radiusRef ) > fSize/2 )
530  {
531  posBox = -posRef + G4ThreeVector(0,0,fSize);
532  if(!isCutted)
533  {
534  pSolidCut = new G4SubtractionSolid("solidCut",
535  pSolidOrbRef,
536  solidBox,
537  &rotMat,
538  posBox);
539  isCutted = true;
540  }
541  else
542  {
543  pSolidCut = new G4SubtractionSolid("solidCut",
544  pSolidCut,
545  solidBox,
546  &rotMat,
547  posBox);
548  }
549  }
550 
551 // Backward
552  if(std::abs( posRef.z() - radiusRef ) > fSize/2 )
553  {
554  posBox = -posRef + G4ThreeVector(0,0,-fSize);
555  if(!isCutted)
556  {
557  pSolidCut = new G4SubtractionSolid("solidCut",
558  pSolidOrbRef,
559  solidBox,
560  &rotMat,
561  posBox);
562  isCutted = true;
563  }
564  else
565  {
566  pSolidCut = new G4SubtractionSolid("solidCut",
567  pSolidCut,
568  solidBox,
569  &rotMat,
570  posBox);
571  }
572  }
573  }
574 
575  for(int i=0, ie=molList.size(); i<ie; ++i)
576  {
577  G4ThreeVector posTar = molList[i].fPosition;
578  G4double rTar = posRef.z();
579  G4double zTar = posTar.z();
580 
581  if(zTar>rTar+20*nm)
582  {
583  break;
584  }
585  else if(zTar<rTar-20*nm)
586  {
587  continue;
588  }
589 
590  G4double radiusTar;
591  if(molList[i].fRadiusWater == 0)
592  {
593  radiusTar = molList[i].fRadius;
594  }
595  else
596  {
597  radiusTar = molList[i].fRadiusWater;
598  }
599 
600  G4double distance = std::abs( (posRef - posTar).getR() );
601 
602  if(distance==0 && !isOurVol)
603  {
604  isOurVol = true;
605  continue;
606  }
607  else if(distance == 0 && isOurVol)
608  {
609  G4ExceptionDescription exceptionDescription;
610  exceptionDescription <<"CreateCutSolid: Two volumes "
611  <<"are placed at the same position.";
612  G4Exception("DNAParser::CreateCutSolid()",
613  "DNAParser004", FatalException,
614  exceptionDescription);
615  }
616  else if(distance <= radiusRef+radiusTar)
617  {
618  G4Box* solidBox = new G4Box("solid_box_for_cut",
619  2*radiusTar,
620  2*radiusTar,
621  2*radiusTar);
622  G4ThreeVector diff = posTar - posRef;
623  G4double d = (std::pow(radiusRef,2) -
624  std::pow(radiusTar,2) +
625  std::pow(distance,2) )/(2*distance) +
626  solidBox->GetZHalfLength() -
627  tinySpace;
628  if(in)
629  {
630  d -= 2*tinySpace;
631  }
632  G4ThreeVector pos = d *( diff/diff.getR() );
633  G4double phi = std::acos(pos.getZ()/pos.getR());
634  G4double theta = std::acos( pos.getX() /
635  ( pos.getR() *
636  std::cos(CLHEP::pi/2. - phi) ) );
637 
638  if(pos.getY()<0)
639  {
640  theta = -theta;
641  }
642 
643  G4ThreeVector rotAxisForPhi(1*nm,0.,0.);
644  rotAxisForPhi.rotateZ(theta + CLHEP::pi/2.);
645 
646  G4RotationMatrix *rotMat = new G4RotationMatrix;
647  rotMat->rotate(-phi, rotAxisForPhi);
648 
649  G4ThreeVector rotZAxis(0.,0.,1*nm);
650  rotMat->rotate(theta, rotZAxis);
651 
652  if(!isCutted)
653  {
654  pSolidCut = new G4SubtractionSolid("solidCut",
655  pSolidOrbRef,
656  solidBox,
657  rotMat,
658  pos);
659  }
660  else
661  {
662  pSolidCut = new G4SubtractionSolid("solidCut",
663  pSolidCut,
664  solidBox,
665  rotMat,
666  pos);
667  }
668  isCutted = true;
669  }
670  }
671 
672  if(isCutted)
673  {
674  return pSolidCut;
675  }
676  else
677  {
678  return pSolidOrbRef;
679  }
680 }
681 
682 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
683 
684 void DNAParser::EnumParser()
685 {
686  fEnumMap["deoxyribose1"] = deoxyribose1;
687  fEnumMap["phosphate1"] = phosphate1;
688  fEnumMap["deoxyribose2"] = deoxyribose2;
689  fEnumMap["phosphate2"] = phosphate2;
690 
691  fEnumMap["base_cytosine"] = base_cytosine;
692  fEnumMap["base_guanine"] = base_guanine;
693  fEnumMap["base_thymine"] = base_thymine;
694  fEnumMap["base_adenine"] = base_adenine;
695 
696  fEnumMap["deoxyribose1_water"] = deoxyribose1_water;
697  fEnumMap["phosphate1_water"] = phosphate1_water;
698  fEnumMap["deoxyribose2_water"] = deoxyribose2_water;
699  fEnumMap["phosphate2_water"] = phosphate2_water;
700 
701 
702  fEnumMap["base_adenine_water"] = base_adenine_water;
703  fEnumMap["base_guanine_water"] = base_guanine_water;
704  fEnumMap["base_cytosine_water"] = base_cytosine_water;
705  fEnumMap["base_thymine_water"] = base_thymine_water;
706 
707  fEnumMap["histone"] = histone;
708  fEnumMap["physWorld"] = physWorld;
709  fEnumMap["VoxelStraight"] = VoxelStraight;
710 }
711