Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/dna/dsbandrepair/src/PhysGeoImport.cc

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

  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 //
 27 /// \file PhysGeoImport.cc
 28 /// \brief Implementation of the PhysGeoImport class
 29 
 30 #include "PhysGeoImport.hh"
 31 #include "Randomize.hh"
 32 #include "G4Sphere.hh"
 33 #include "G4PhysicalConstants.hh"
 34 
 35 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 36 
 37 PhysGeoImport::PhysGeoImport() :
 38     fFactor(1.),
 39     fGeoName("")
 40 {
 41     DefineMaterial();
 42 }
 43 
 44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 45 
 46 PhysGeoImport::PhysGeoImport(G4bool isVisu) :
 47     fIsVisu(isVisu),
 48     fFactor(1.),
 49     fGeoName("")
 50 {
 51     if (fIsVisu) G4cout<<" For Checking Visualization!"<<G4endl;
 52     DefineMaterial();
 53 }
 54 
 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 56 
 57 G4LogicalVolume* PhysGeoImport::CreateLogicVolume(const G4String& fileName,G4String& voxelName)
 58 {
 59     // The idea is to return a pointer to the mother volume of the input file as output.
 60     // Create a temporary material map to assign materials
 61 
 62     std::map<G4String, G4Material*> materials;
 63 
 64     // Loop on each material to build the map
 65     for(G4int i=0, ie=fMaterialVect.size(); i<ie; ++i)
 66     {
 67         G4String matName("");
 68 
 69         if(fMaterialVect[i]->GetName()=="adenine_PU") matName = "adenine";
 70         else if(fMaterialVect[i]->GetName()=="thymine_PY") matName = "thymine";
 71         else if(fMaterialVect[i]->GetName()=="guanine_PU") matName = "guanine";
 72         else if(fMaterialVect[i]->GetName()=="cytosine_PY") matName = "cytosine";
 73         else if(fMaterialVect[i]->GetName()=="backbone_THF") matName = "deoxyribose";
 74         else if(fMaterialVect[i]->GetName()=="backbone_TMP") matName = "phosphate";
 75         else if(fMaterialVect[i]->GetName()=="G4_WATER") matName = "water";
 76         else if(fMaterialVect[i]->GetName()=="homogeneous_dna") matName = "homogeneous_dna";
 77         else if(fMaterialVect[i]->GetName()=="G4_Galactic") matName = "vacuum";
 78         else
 79         {
 80             matName = fMaterialVect[i]->GetName();
 81         }
 82 
 83         materials[matName] = fMaterialVect[i];
 84     }
 85 
 86 
 87     // Parse the input file
 88 
 89     voxelName = ParseFile(fileName);
 90 
 91     // Create the volumes
 92 
 93     G4String boxNameSolid = fGeoName+"_solid";
 94     G4Box* box_solid = new G4Box(boxNameSolid, fSize/2, fSize/2, fSize/2);
 95 
 96     G4String boxNameLogic = fGeoName+"_logic";
 97     G4LogicalVolume* box_logic = new G4LogicalVolume(box_solid, fpWater, boxNameLogic);
 98 
 99     // sort the molecules
100     std::sort(fMolecules.begin(), fMolecules.end() );
101    
102     // Loop on all the parsed molecules
103     for(G4int i=0, ie=fMolecules.size(); i<ie; ++i)
104     {
105         // Retrieve general molecule informations
106         //
107         G4String name = fMolecules[i].fName;
108         G4String materialName = "water";// fMolecules[i].fMaterial;
109         G4double radius = fMolecules[i].fRadius;
110         G4double waterRadius = fMolecules[i].fRadiusWater;
111         G4ThreeVector moleculePosition = fMolecules[i].fPosition;
112         G4int copyNum = fMolecules[i].fCopyNumber;
113 
114         // Water hydration shell volume part
115 
116         G4Orb* moleculeWater_solid = 0;
117         G4VSolid* moleculeWaterCut_solid = 0;
118         G4LogicalVolume* moleculeWater_logic = 0;
119 
120         // If water radius != 0 then we have a water hydration shell
121         G4double tol = 0.0001;
122         if(waterRadius > (0 + tol)*nm)
123         {
124             G4Material* mat = 0;
125             
126             mat=materials[materialName];
127 
128             G4String nameWaterSolid = name+"_water_solid";
129             moleculeWater_solid = new G4Orb(nameWaterSolid, waterRadius);
130             moleculeWaterCut_solid = CreateCutSolid(moleculeWater_solid, fMolecules[i], fMolecules, false);
131 
132             G4String nameWaterLogic = name+"_water_logic";
133             moleculeWater_logic = new G4LogicalVolume(moleculeWaterCut_solid, mat, nameWaterLogic);
134 
135             G4String nameWaterPhys = name+"_water_phys";
136             new G4PVPlacement(0, moleculePosition, moleculeWater_logic, nameWaterPhys, box_logic, false, copyNum);
137         }
138 
139         // Dna volume part
140 
141         G4Orb* molecule_solid = 0;
142         G4VSolid* moleculeCut_solid = 0;
143         G4LogicalVolume* molecule_logic = 0;
144 
145         G4String nameSolid = fMolecules[i].fName+"_solid";
146         molecule_solid = new G4Orb(nameSolid, radius);
147         moleculeCut_solid = CreateCutSolid(molecule_solid, fMolecules[i], fMolecules, true);
148 
149         G4String nameLogic = name+"_logic";
150         molecule_logic = new G4LogicalVolume(moleculeCut_solid, materials[materialName], nameLogic);
151 
152         // If there was a water hydration shell volume then the current dna volume is
153         // placed within it and, thus, its relative coordinates are 0,0,0.
154         if(waterRadius > (0 + tol)*nm)
155         {
156             G4ThreeVector position(0.,0.,0.);
157             G4String namePhys = name+"_phys";
158             new G4PVPlacement(0, position, molecule_logic, namePhys, moleculeWater_logic, false, copyNum);
159         }
160         // If not, coordinates are those of the molecule
161         else
162         {
163             G4ThreeVector position = moleculePosition;
164             G4String namePhys = name+"_phys";
165             new G4PVPlacement(0, position, molecule_logic, namePhys, box_logic, false, copyNum);
166         }
167     }
168 
169     // Clear the containers
170     fMolecules.clear();
171     fRadiusMap.clear();
172     fWaterRadiusMap.clear();
173 
174     return box_logic;
175 }
176 
177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
178 
179 G4String PhysGeoImport::ParseFile(const G4String& fileName)
180 {
181     G4cout<<"Start parsing of "<<fileName<<G4endl;
182 
183     // Clear the containers
184     fMolecules.clear();
185     fRadiusMap.clear();
186     fWaterRadiusMap.clear();
187 
188     // Setup the input stream
189     std::ifstream file(fileName.c_str());
190 
191     // Check if the file was correctly opened
192     if(!file.is_open())
193     {
194         // Geant4 exception
195         G4String msg = fileName+" could not be opened";
196         G4Exception("PhysGeoImport::ParseFile()", "Geo_InputFileNotOpened", FatalException, msg);
197     }
198 
199     // Define the line string variable
200     G4String line;
201     G4int preBpCpNo = -1;
202     G4int preHistoneCpNo =-1;
203     G4int NoBp = 0;
204     G4int NoHistone = 0;
205     // Read the file line per line
206     while(std::getline(file, line) )
207     {
208         // Check the line to determine if it is empty
209         if(line.empty()) continue; // skip the line if it is empty
210 
211         // Data string stream
212         std::istringstream issLine(line);
213 
214         // String to determine the first letter/word
215         G4String firstItem;
216 
217         // Put the first letter/word within the string
218         issLine >> firstItem;
219 
220         // Check first letter to determine if the line is data or comment
221         if(firstItem=="#") continue; // skip the line if it is comment
222 
223         // Use the file
224         else if(firstItem=="_Name")
225         {
226             G4String name;
227             issLine >> name;
228 
229             fGeoName = name;
230         }
231         else if(firstItem=="_Size")
232         {
233             G4double size;
234             issLine >> size;
235             size *= fFactor*nm;
236 
237             fSize = size;
238         }
239         else if(firstItem == "_Version")
240         {
241 
242         }
243         else if(firstItem=="_Number")
244         {
245 
246         }
247         else if(firstItem=="_Radius")
248         {
249             G4String name;
250             issLine >> name;
251 
252             G4double radius;
253             issLine >> radius;
254             radius *= fFactor*nm;
255 
256             G4double waterRadius;
257             issLine >> waterRadius;
258             waterRadius *= fFactor*nm;
259 
260             fRadiusMap[name] = radius;
261             fWaterRadiusMap[name] = waterRadius;
262         }
263         else if(firstItem=="_pl")
264         {
265             G4String name;
266             issLine >> name;
267 
268             G4String material;
269             issLine >> material;
270 
271             G4int strand;
272             issLine >> strand;
273 
274             G4int copyNumber;
275             issLine >> copyNumber;
276 
277             if (name == "histone") {
278                 if (copyNumber != preHistoneCpNo ) {
279                     NoHistone ++;
280                     preHistoneCpNo = copyNumber;
281                 }
282             } else {
283                 if (copyNumber != preBpCpNo) {
284                     NoBp++;
285                     preBpCpNo = copyNumber;
286                 }
287             }
288 
289             G4double x;
290             issLine >> x;
291             x *= fFactor*nm;
292 
293             G4double y;
294             issLine >> y;
295             y *= fFactor*nm;
296 
297             G4double z;
298             issLine >> z;
299             z *= fFactor*nm;
300 
301             Molecule molecule(name, copyNumber, G4ThreeVector(x, y, z), fRadiusMap[name], fWaterRadiusMap[name], material, strand);
302 
303             fMolecules.push_back(molecule);
304         }
305         else
306         {
307             // Geant4 exception
308             G4String msg = firstItem+" is not defined in the parser. Check the input file: "+fileName+".";
309             G4Exception("PhysGeoImport::ParseFile()", "Geo_WrongParse", FatalException, msg);
310         }
311     }
312 
313     // Close the file once the reading is done
314     file.close();
315 
316     G4cout<<"End parsing of "<<fileName<<G4endl;
317     fVoxelNbHistoneMap.insert({fGeoName,NoHistone});
318     fVoxelNbBpMap.insert({fGeoName,NoBp});
319     return fGeoName;
320 }
321 
322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
323 
324 G4VSolid* PhysGeoImport::CreateCutSolid(G4Orb *solidOrbRef,
325                                                Molecule &molRef,
326                                                std::vector<Molecule> &molList,
327                                                G4bool in)
328 {
329     // The idea behing this method is to cut overlap volumes by selecting one of them (the reference) and checking all the other volumes (the targets).
330     // If a reference and a target volumes are close enough to overlap they will be cut.
331     // The reference is already selected when we enter this method.
332 
333     // Use the tiny space to differentiate the frontiers (may not be necessary)
334     G4double tinySpace = 0.001*fFactor*nm;
335 
336     // Cutted solid to be returned
337     G4SubtractionSolid* solidCut(NULL);
338 
339     // Some flags
340     G4bool isCutted = false;
341     G4bool isOurVol = false;
342 
343     // Radius of the molecule to cut
344     G4double radiusRef;
345     if(molRef.fRadiusWater==0) radiusRef = molRef.fRadius;
346     else radiusRef = molRef.fRadiusWater;
347 
348     // Reference volume position
349     G4ThreeVector posRef = molRef.fPosition;
350     // Before looping on all the volumes we check if the current 
351     //reference volume overlaps with its container voxel boundaries.
352     if(std::abs(posRef.x() ) + radiusRef > fSize/2 // along x
353             || std::abs(posRef.y() ) + radiusRef > fSize/2 // along y
354             || std::abs(posRef.z() ) + radiusRef > fSize/2) // along z
355     {
356         // If we enter here, then the reference volume overlaps with the boundaries of the container voxel
357         // Box used to cut
358         G4Box* solidBox = new G4Box("solid_box_for_cut", fSize/2, fSize/2, fSize/2);
359         G4ThreeVector posBox;
360 
361         // Create a dummy rotation matrix
362         G4RotationMatrix *rotMat = new G4RotationMatrix;
363         rotMat->rotate(0, G4ThreeVector(0,0,1) );
364 
365         // To choose the cut direction
366 
367         // Up
368         if(std::abs( posRef.y() + radiusRef ) > fSize/2 )
369         {
370            posBox = -posRef +  G4ThreeVector(0,fSize,0);
371 
372             // If the volume is cutted for the first time
373             if(!isCutted)
374             {
375                 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox);
376                 isCutted = true;
377             }
378             // For the other times
379             else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox);
380         }
381 
382         // Down
383         if(std::abs( posRef.y() - radiusRef ) > fSize/2 )
384         {
385             posBox = -posRef + G4ThreeVector(0,-fSize,0);
386 
387             // If the volume is cutted for the first time
388             if(!isCutted)
389             {
390                 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef,solidBox,rotMat, posBox);
391             isCutted = true;
392         }
393             // For the other times
394             else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox);
395         }
396 
397         // Left
398         if(std::abs( posRef.x() + radiusRef ) > fSize/2 )
399         {
400             posBox = -posRef + G4ThreeVector(fSize,0,0);
401 
402             // If the volume is cutted for the first time
403             if(!isCutted)
404             {
405                 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox);
406                 isCutted = true;
407             }
408             // For the other times
409             else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox);
410         }
411 
412         // Right
413         if(std::abs( posRef.x() - radiusRef ) > fSize/2 )
414         {
415             posBox = -posRef + G4ThreeVector(-fSize,0,0);
416 
417             // If the volume is cutted for the first time
418             if(!isCutted)
419             {
420                 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox);
421                 isCutted = true;
422             }
423             // For the other times
424             else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox);
425         }
426 
427         // Forward
428         if(std::abs( posRef.z() + radiusRef ) > fSize/2 )
429         {
430             posBox = -posRef + G4ThreeVector(0,0,fSize);
431 
432             // If the volume is cutted for the first time
433             if(!isCutted)
434             {
435                 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox);
436                 isCutted = true;
437             }
438             // For the other times
439             else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox);
440         }
441 
442         // Backward
443         if(std::abs( posRef.z() - radiusRef ) > fSize/2 )
444         {
445             posBox = -posRef + G4ThreeVector(0,0,-fSize);
446 
447             // If the volume is cutted for the first time
448             if(!isCutted)
449             {
450                 solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, solidBox, rotMat, posBox);
451                 isCutted = true;
452             }
453             // For the other times
454             else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, posBox);
455         }
456     }
457 
458     // Look the other volumes of the voxel
459     // Loop on all the target volumes (other volumes with potential overlaps)
460     for(G4int i=0, ie=molList.size(); i<ie; ++i)
461     {
462         G4ThreeVector posTar = molList[i].fPosition;
463 
464         G4double rTar = posRef.z();
465         G4double zTar = posTar.z();
466 
467         if(zTar>rTar+20*fFactor*nm)
468         {
469             break;
470         }
471         else if(zTar<rTar-20*fFactor*nm)
472         {
473             continue;
474         }
475 
476         // Retrieve current target sphere informations
477         G4double radiusTar;
478         if(molList[i].fRadiusWater==0) radiusTar = molList[i].fRadius;
479         else radiusTar = molList[i].fRadiusWater;
480 
481         // Compute the distance reference-target
482         G4double distance = std::abs( (posRef - posTar).getR() );
483 
484         // Use the distance to check if the current target is also the reference.
485         // This can only happen once per loop.
486         if(distance==0 && !isOurVol)
487         {
488             // Target volume is also reference volume.
489 
490             // Set the flag
491             isOurVol = true;
492 
493             // Next iteration
494             continue;
495         }
496         // If the condition is correct more than one time then there is a mistake somewhere.
497         else if(distance == 0 && isOurVol)
498         {
499             G4ExceptionDescription desmsg;
500             desmsg << "DetectorConstruction::CreateCutSolid: Two volumes are placed at the same position.";
501             G4Exception("ChemGeoImport::BuildGeometry", "", FatalException, desmsg);
502         }
503 
504         // If the volumes are differents then we want to know if they are
505         // close enough to overlap and, thus, to intiate a cut.
506         else if(distance <= radiusRef+radiusTar)
507         {
508             // Volumes are close enough, there will be a cut
509 
510             // Box used to cut
511             G4Box* solidBox = new G4Box("solid_box_for_cut", 2*radiusTar, 2*radiusTar, 2*radiusTar);
512 
513             // This part is tricky.
514             // The goal is to calculate the position of the intersection center
515 
516             // diff vector to from ref to tar
517             G4ThreeVector diff = posTar - posRef;
518 
519             // Find the intersection point and add to it half the length of the box used to cut
520             G4double d = (std::pow(radiusRef,2)-std::pow(radiusTar,2)+std::pow(distance,2) ) / 
521                         (2*distance) + solidBox->GetZHalfLength() - tinySpace;
522 
523             // If we are in another volume we choose to double the tinySpace to 
524             //differentiate without ambiguities the inner and outer volume frontiers.
525             // (may not be necessary)
526             if(in) d -= 2*tinySpace;
527 
528             // Position of the box in order to achieve the cut.
529             // "* ( diff/diff.getR() )" is necessary to get a vector in the right direction as output
530             G4ThreeVector pos = d *( diff/diff.getR() );
531 
532             // Get the rotation angles because the box used to cut needs to be rotated
533             // to give the right "cut angle".
534             G4double phi = std::acos(pos.getZ()/pos.getR());
535             G4double theta = std::acos( pos.getX() / ( pos.getR()*std::cos(pi/2.-phi) ) );
536 
537             if(pos.getY()<0) theta = -theta;
538 
539             G4ThreeVector rotAxisForPhi(1*fFactor*nm,0.,0.);
540             rotAxisForPhi.rotateZ(theta+pi/2);
541 
542             // Create the rotation matrix
543             G4RotationMatrix *rotMat = new G4RotationMatrix;
544             rotMat->rotate(-phi, rotAxisForPhi);
545 
546             // Rotate it again
547             G4ThreeVector rotZAxis(0.,0.,1*fFactor*nm);
548             rotMat->rotate(theta, rotZAxis);
549 
550             // If the volume is cutted for the first time
551             if(!isCutted) solidCut = new G4SubtractionSolid("solidCut", solidOrbRef, 
552             solidBox, rotMat, pos);
553 
554             // For the other times
555             else solidCut = new G4SubtractionSolid("solidCut", solidCut, solidBox, rotMat, pos);
556 
557             // Set the cut flag
558             isCutted = true;
559         }
560     }
561 
562     // If there was at least one cut then we return the cutted volume
563     if(isCutted) return solidCut;
564 
565     // Otherwise, we return the original volume
566     else return solidOrbRef;
567     
568 }
569 
570 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
571 
572 void PhysGeoImport::DefineMaterial()
573 {
574     G4NistManager * man = G4NistManager::Instance();
575     fpWater = man->FindOrBuildMaterial("G4_WATER");
576     fMaterialVect.push_back(fpWater);
577     fVacuum = man->FindOrBuildMaterial("G4_Galactic");
578     fMaterialVect.push_back(fVacuum);
579 
580     G4double z, a, density;
581     G4String name, symbol;
582     G4int nComponents, nAtoms;
583 
584     a = 12.0107*g/mole;
585     G4Element* elC = new G4Element(name="Carbon", symbol="C", z=6., a);
586     a = 1.00794*g/mole;
587     G4Element* elH  = new G4Element(name="Hydrogen",symbol="H" , z= 1., a);
588     a = 15.9994*g/mole;
589     G4Element* elO  = new G4Element(name="Oxygen"  ,symbol="O" , z= 8., a);
590     a = 14.00674*g/mole;
591     G4Element* elN  = new G4Element(name="Nitrogen"  ,symbol="N" , z= 7., a);
592     a = 30.973762*g/mole;
593     G4Element* elP  = new G4Element(name="Phosphorus"  ,symbol="P" , z= 15., a);
594 
595     density = 1.*g/cm3;
596 
597     fTHF = new G4Material("THF", density, nComponents=3);
598     fTHF->AddElement(elC, nAtoms=4);
599     fTHF->AddElement(elH, nAtoms=8);
600     fTHF->AddElement(elO, nAtoms=1);
601 
602     fPY = new G4Material("PY", density, nComponents=3);
603     fPY->AddElement(elC, nAtoms=4);
604     fPY->AddElement(elH, nAtoms=4);
605     fPY->AddElement(elN, nAtoms=2);
606 
607     fPU = new G4Material("PU", density, nComponents=3);
608     fPU->AddElement(elC, nAtoms=5);
609     fPU->AddElement(elH, nAtoms=4);
610     fPU->AddElement(elN, nAtoms=4);
611 
612     fTMP = new G4Material("TMP", density,4);
613     fTMP->AddElement(elC, nAtoms=3);
614     fTMP->AddElement(elH, nAtoms=6);
615     fTMP->AddElement(elP, nAtoms=1);
616     fTMP->AddElement(elO, nAtoms=4);
617 
618     fMaterialVect.push_back(fTHF);
619     fMaterialVect.push_back(fPY);
620     fMaterialVect.push_back(fPU);
621     fMaterialVect.push_back(fTMP);
622 
623     // Mixture creation
624     fSugarMixt = new G4Material("sugarMixt", density, nComponents=2);
625     fSugarMixt->AddMaterial(fTHF,0.5);
626     fSugarMixt->AddMaterial(fTMP,0.5);
627     fMaterialVect.push_back(fSugarMixt);
628 
629     // Real DNA materials creation
630     //
631     // Density calculation:
632     //
633     // d=m/V and M=m/n <=> m = M*n
634     // <=> d = M*n/V
635     // <=> d = (M*nb)/(Na*V)
636 
637     G4double MBackbone = (5.*12.0104 + 10.*1.00794 + 4.*15.9994)*g/mole;
638     G4double VBackbone = 1.823912/20. * std::pow(1.e-9, 3.) *m3;
639     G4double densityBaTHF = (MBackbone/(g/mole) * 1. )/(6.022e23*VBackbone/cm3) *g/cm3; // g/cm3
640 
641     fDeoxyribose = new G4Material("backbone_THF", densityBaTHF, nComponents=3);
642     fDeoxyribose->AddElement(elC, nAtoms=5);
643     fDeoxyribose->AddElement(elH, nAtoms=10);
644     fDeoxyribose->AddElement(elO, nAtoms=4);
645     fMaterialVect.push_back(fDeoxyribose);
646 
647     MBackbone = (3.*1.00794 + 1.*30.973762 + 4.*15.9994)*g/mole;
648     VBackbone = 1.19114/20. * std::pow(1.e-9, 3.) *m3;
649     G4double densityBaTMP = (MBackbone/(g/mole) * 1. )/(6.022e23*VBackbone/cm3) *g/cm3; // g/cm3
650 
651     fPhosphate = new G4Material("backbone_TMP", densityBaTMP,3);
652     fPhosphate->AddElement(elH, nAtoms=3);
653     fPhosphate->AddElement(elP, nAtoms=1);
654     fPhosphate->AddElement(elO, nAtoms=4);
655     fMaterialVect.push_back(fPhosphate);
656 
657     G4double MCytosine = (4.*12.0104 + 5.*1.00794 + 3.*14.0067 + 1.*15.9994)*g/mole;
658     G4double VBase = 1.8527205/20. * std::pow(1.e-9, 3.) *m3;
659     G4double densityCy = (MCytosine/(g/mole) * 1. )/(6.022e23*VBase/cm3) *g/cm3; // g/cm3
660 
661     fCytosine_PY = new G4Material("cytosine_PY", densityCy, nComponents=4);
662     fCytosine_PY->AddElement(elC, nAtoms=4);
663     fCytosine_PY->AddElement(elH, nAtoms=5);
664     fCytosine_PY->AddElement(elN, nAtoms=3);
665     fCytosine_PY->AddElement(elO, nAtoms=1);
666     fMaterialVect.push_back(fCytosine_PY);
667 
668     G4double MThy = (5.*12.0104 + 6.*1.00794 + 2.*14.0067 + 2.*15.9994)*g/mole;
669     VBase = 1.8527205/20. * std::pow(1.e-9, 3.) *m3;
670     densityCy = (MThy/(g/mole) * 1. )/(6.022e23*VBase/cm3) *g/cm3; // g/cm3
671 
672     fThymine_PY = new G4Material("thymine_PY", densityCy, nComponents=4);
673     fThymine_PY->AddElement(elC, nAtoms=5);
674     fThymine_PY->AddElement(elH, nAtoms=6);
675     fThymine_PY->AddElement(elN, nAtoms=2);
676     fThymine_PY->AddElement(elO, nAtoms=2);
677     fMaterialVect.push_back(fThymine_PY);
678 
679     G4double MAdenine = (5.*12.0104 + 5.*1.00794 + 5.*14.0067)*g/mole;
680     VBase = 1.8527205/20. * std::pow(1.e-9, 3.) *m3;
681     G4double densityAd = (MAdenine/(g/mole) * 1. )/(6.022e23*VBase/cm3) *g/cm3; // g/cm3
682 
683     fAdenine_PU = new G4Material("adenine_PU", densityAd, nComponents=3);
684     fAdenine_PU->AddElement(elC, nAtoms=5);
685     fAdenine_PU->AddElement(elH, nAtoms=5);
686     fAdenine_PU->AddElement(elN, nAtoms=5);
687     fMaterialVect.push_back(fAdenine_PU);
688 
689     MAdenine = (5.*12.0104 + 5.*1.00794 + 5.*14.0067 + 1.*15.999)*g/mole;
690     VBase = 1.8527205/20. * std::pow(1.e-9, 3.) *m3;
691     densityAd = (MAdenine/(g/mole) * 1. )/(6.022e23*VBase/cm3) *g/cm3; // g/cm3
692 
693     fGuanine_PU = new G4Material("guanine_PU", densityAd, nComponents=4);
694     fGuanine_PU->AddElement(elC, nAtoms=5);
695     fGuanine_PU->AddElement(elH, nAtoms=5);
696     fGuanine_PU->AddElement(elN, nAtoms=5);
697     fGuanine_PU->AddElement(elO, nAtoms=1);
698     fMaterialVect.push_back(fGuanine_PU);
699 
700     fHomogeneous_dna = new G4Material("homogeneous_dna", 19.23e-21/(12.30781*1.e-21)  *g/cm3, nComponents=7);//21
701     fHomogeneous_dna->AddMaterial(fpWater, 0.37);
702     fHomogeneous_dna->AddMaterial(fDeoxyribose, 0.23);
703     fHomogeneous_dna->AddMaterial(fPhosphate, 0.17);
704     fHomogeneous_dna->AddMaterial(fAdenine_PU, 0.06);
705     fHomogeneous_dna->AddMaterial(fGuanine_PU, 0.07);
706     fHomogeneous_dna->AddMaterial(fCytosine_PY, 0.05);
707     fHomogeneous_dna->AddMaterial(fThymine_PY, 0.05);
708 
709     fMaterialVect.push_back(fHomogeneous_dna);
710 }
711 
712 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
713 
714 G4LogicalVolume* PhysGeoImport::CreateNucleusLogicVolume(const G4String& fileName)
715 {
716     // Here, we parse the nucleus file and build a logical Geant4 volume from it.
717 
718     G4cout<<"Start the nucleus creation..."<<G4endl;
719 
720     G4VSolid* nucleusSolid = 0;
721     G4LogicalVolume* nucleusLogic = 0;
722 
723     G4bool isNucleusBuild (false);
724 
725     // Open the world file
726     std::ifstream nucleusFile;
727     nucleusFile.open(fileName.c_str());
728 
729     // Check if the file was correctly opened
730     if(!nucleusFile.is_open())
731     {
732         // Geant4 exception
733         G4String msg = fileName+" could not be opened";
734         G4Exception("PhysGeoImport::ParseFile()", "Geo_InputFileNotOpened", FatalException, msg);
735     }
736 
737     // Define the line string variable
738     G4String line;
739 
740     // Read the file line per line
741     while(std::getline(nucleusFile, line) && !isNucleusBuild)
742     {
743         // Check the line to determine if it is empty
744         if(line.empty()) continue; // skip the line if it is empty
745 
746         // Data string stream
747         std::istringstream issLine(line);
748 
749         // String to determine the first letter/word
750         G4String firstItem;
751 
752         // Put the first letter/word within the string
753         issLine >> firstItem;
754 
755         // Check first letter to determine if the line is data or comment
756         if(firstItem=="#") continue; // skip the line if it is comment
757 
758         // Use the file
759 
760         else if(firstItem == "_Name")
761         {
762             issLine >> fNucleusName;
763         }
764 
765         else if(firstItem=="_Type")
766         {
767             issLine >> fNucleusType;
768 
769             if (fNucleusType == "Ellipsoid")
770             {
771                 G4float semiX, semiY, semiZ;
772 
773                 issLine >> semiX >> semiY >> semiZ;
774 
775                 semiX *= fFactor*nm;
776                 semiY *= fFactor*nm;
777                 semiZ *= fFactor*nm;
778 
779                 fNucleusData["SemiX"] = semiX;
780                 fNucleusData["SemiY"] = semiY;
781                 fNucleusData["SemiZ"] = semiZ;
782 
783                 nucleusSolid = new G4Ellipsoid(fNucleusName, semiX, semiY, semiZ);
784                 nucleusLogic = new G4LogicalVolume(nucleusSolid, fpWater, "nucleus_logic");
785                 G4double r3 = semiX*semiY*semiZ/m3;
786                 fNucleusVolume = 4.0*pi*r3/3.0;
787             } else if (fNucleusType == "EllipticCylinder") {
788               G4float semiX, semiY, semiZ;
789     
790                 issLine >> semiX >> semiY >> semiZ;
791     
792                 semiX *= fFactor*nm;
793                 semiY *= fFactor*nm;
794                 semiZ *= fFactor*nm;
795     
796                 fNucleusData["SemiX"] = semiX;
797                 fNucleusData["SemiY"] = semiY;
798                 fNucleusData["SemiZ"] = semiZ;
799 
800                 nucleusSolid = new G4EllipticalTube(fNucleusName, semiX, semiY, semiZ);
801                 nucleusLogic = new G4LogicalVolume(nucleusSolid, fpWater, "nucleus_logic");
802                 G4double r3 = 2*semiX*semiY*semiZ/m3; // multiplied by 2 cuz semiZ is hafl of height
803                 fNucleusVolume = pi*r3;
804             } else if (fNucleusType == "Spherical") {
805                 G4float radius;
806                 issLine >> radius;
807                 radius *= fFactor*nm;
808                 fNucleusData["SemiX"] = radius;
809                 fNucleusData["SemiY"] = radius;
810                 fNucleusData["SemiZ"] = radius;
811                 nucleusSolid = new G4Orb(fNucleusName,radius);
812                 nucleusLogic = new G4LogicalVolume(nucleusSolid, fpWater, "nucleus_logic");
813                 G4double r3 = radius*radius*radius/m3;
814                 fNucleusVolume = 4.0*pi*r3/3.0;
815             } else {
816                 G4String msg =fNucleusType+" is not registered.";
817                 G4Exception("PhysGeoImport::CreateNucleusLogicVolume()", 
818                 "Geo_InputFile_Unknown_Nucleus", FatalException, msg);
819             }
820 
821             isNucleusBuild = true;
822         }
823     }
824 
825     nucleusFile.close();
826 
827     if(!isNucleusBuild)
828     {
829         G4String msg = "Nucleus data were not found for "+fNucleusType;
830         G4Exception("PhysGeoImport::CreateNucleusLogicVolume()", 
831         "Geo_InputFile_Unknown_Nucleus", FatalException, msg);
832     }
833 
834     return nucleusLogic;
835 }
836 
837 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
838 
839 std::vector<Voxel>* PhysGeoImport::CreateVoxelsData(const G4String& fileName)
840 {
841     fTotalNbBpPlacedInGeo = 0;
842     fTotalNbHistonePlacedInGeo = 0;
843     G4cout<<"Start the voxel data generation..."<<G4endl;
844 
845     std::vector<Voxel>* voxels = new std::vector<Voxel>;
846 
847     // Clear any previously loaded world
848     voxels->clear();
849 
850     voxels->reserve(1000000);
851 
852     // Open the world file
853     std::ifstream nucleusFile;
854     nucleusFile.open(fileName.c_str());
855 
856     // Check if the file was correctly opened
857     if(!nucleusFile.is_open())
858     {
859         // Geant4 exception
860         G4String msg = fileName+" could not be opened";
861         G4Exception("PhysGeoImport::ParseFile()", "Geo_InputFileNotOpened", FatalException, msg);
862     }
863 
864     // Initialise the copy number to 0
865     G4int voxelCopyNumber = 0;
866 
867     // Define the line string variable
868     G4String line;
869 
870     auto whatchromatintype = [] (Voxel::VoxelType voxt) -> ChromatinType {
871         if (voxt == Voxel::Straight || voxt == Voxel::Left || voxt == Voxel::Right ||
872             voxt == Voxel::Up || voxt == Voxel::Down) return ChromatinType::fHeterochromatin;
873         else if (voxt == Voxel::Straight2 || voxt == Voxel::Left2 || voxt == Voxel::Right2 ||
874             voxt == Voxel::Up2 || voxt == Voxel::Down2) return ChromatinType::fEuchromatin;
875         else return ChromatinType::fUnspecified;   
876     };
877     // Read the file line per line
878     while(std::getline(nucleusFile, line) )
879     {
880         // Check the line to determine if it is empty
881         if(line.empty()) continue; // skip the line if it is empty
882 
883         // Data string stream
884         std::istringstream issLine(line);
885 
886         // String to determine the first letter/word
887         G4String firstItem;
888 
889         // Put the first letter/word within the string
890         issLine >> firstItem;
891 
892         // Check first letter to determine if the line is data or comment
893         if(firstItem=="#") continue; // skip the line if it is comment
894 
895         // Use the file
896 
897         else if(firstItem == "_pl")
898         {
899             // Set the flags
900 
901             G4String name;
902             issLine >> name;
903 
904             G4int chromoCpN;
905             issLine >> chromoCpN;
906 
907             G4int domainCpN;
908             issLine >> domainCpN;
909 
910             G4double x;
911             issLine >> x;
912             x *= fFactor*nm;
913 
914             G4double y;
915             issLine >> y;
916             y *= fFactor*nm;
917 
918             G4double z;
919             issLine >> z;
920             z *= fFactor*nm;
921             
922             G4double xx;
923             issLine >> xx;
924             
925             G4double yx;
926             issLine >> yx;
927             
928             G4double zx;
929             issLine >> zx;
930             
931             G4double xy;
932             issLine >> xy;
933             
934             G4double yy;
935             issLine >> yy;
936             
937             G4double zy;
938             issLine >> zy;
939             
940             G4double xz;
941             issLine >> xz;
942             
943             G4double yz;
944             issLine >> yz;
945             
946             G4double zz;
947             issLine >> zz;
948             
949 
950             G4RotationMatrix* rot = new G4RotationMatrix(G4ThreeVector(xx,xy,xz),
951             G4ThreeVector(yx,yy,yz),G4ThreeVector(zx,zy,zz));
952 
953             Voxel::VoxelType type = Voxel::Other;
954 
955             if(name=="voxelStraight" || name=="VoxelStraight")
956             {
957                 type = Voxel::Straight;
958             }
959             else if(name=="voxelUp" || name=="VoxelUp")
960             {
961                 type = Voxel::Up;
962             }
963             else if(name=="voxelDown" || name=="VoxelDown")
964             {
965                 type = Voxel::Down;
966             }
967             else if(name=="voxelRight" || name=="VoxelRight")
968             {
969                 type = Voxel::Right;
970             }
971             else if(name=="voxelLeft" || name=="VoxelLeft")
972             {
973                 type = Voxel::Left;
974             }
975             else if(name=="voxelStraight2" || name=="VoxelStraight2")
976             {
977                 type = Voxel::Straight2;
978             }
979             else if(name=="voxelUp2" || name=="VoxelUp2")
980             {
981                 type = Voxel::Up2;
982             }
983             else if(name=="voxelDown2" || name=="VoxelDown2")
984             {
985                 type = Voxel::Down2;
986             }
987             else if(name=="voxelRight2" || name=="VoxelRight2")
988             {
989                 type = Voxel::Right2;
990             }
991             else if(name=="voxelLeft2" || name=="VoxelLeft2")
992             {
993                 type = Voxel::Left2;
994             }
995             else
996             {
997                 G4ExceptionDescription msg ;
998                 msg << "The voxel named "<<name<<" is not registered in the parser";
999                 G4Exception("PhysGeoImport::CreateVoxelsData", "", FatalException, msg, "");
1000             }
1001 
1002             voxels->push_back(Voxel(voxelCopyNumber, chromoCpN, domainCpN, type, 
1003             G4ThreeVector(x,y,z), rot) );
1004             fTotalNbBpPlacedInGeo += fVoxelNbBpMap[name];
1005             fTotalNbHistonePlacedInGeo += fVoxelNbHistoneMap[name];
1006             voxelCopyNumber++;
1007             auto chromatintype =  whatchromatintype(type);
1008             if (fChromatinTypeCount.find(chromatintype) == fChromatinTypeCount.end()) 
1009                 fChromatinTypeCount.insert({chromatintype,1});
1010             else fChromatinTypeCount[chromatintype]++;
1011         }
1012     }
1013 
1014     nucleusFile.close();
1015 
1016     G4cout<<"End the voxel data generation"<<G4endl;
1017 
1018     return voxels;
1019 }
1020 
1021 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......