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 ]

Diff markup

Differences between /examples/advanced/dna/dsbandrepair/src/PhysGeoImport.cc (Version 11.3.0) and /examples/advanced/dna/dsbandrepair/src/PhysGeoImport.cc (Version 10.2)


  1 //                                                  1 
  2 // *******************************************    
  3 // * License and Disclaimer                       
  4 // *                                              
  5 // * The  Geant4 software  is  copyright of th    
  6 // * the Geant4 Collaboration.  It is provided    
  7 // * conditions of the Geant4 Software License    
  8 // * LICENSE and available at  http://cern.ch/    
  9 // * include a list of copyright holders.         
 10 // *                                              
 11 // * Neither the authors of this software syst    
 12 // * institutes,nor the agencies providing fin    
 13 // * work  make  any representation or  warran    
 14 // * regarding  this  software system or assum    
 15 // * use.  Please see the license in the file     
 16 // * for the full disclaimer and the limitatio    
 17 // *                                              
 18 // * This  code  implementation is the result     
 19 // * technical work of the GEANT4 collaboratio    
 20 // * By using,  copying,  modifying or  distri    
 21 // * any work based  on the software)  you  ag    
 22 // * use  in  resulting  scientific  publicati    
 23 // * acceptance of all terms of the Geant4 Sof    
 24 // *******************************************    
 25 //                                                
 26 //                                                
 27 /// \file PhysGeoImport.cc                        
 28 /// \brief Implementation of the PhysGeoImport    
 29                                                   
 30 #include "PhysGeoImport.hh"                       
 31 #include "Randomize.hh"                           
 32 #include "G4Sphere.hh"                            
 33 #include "G4PhysicalConstants.hh"                 
 34                                                   
 35 //....oooOO0OOooo........oooOO0OOooo........oo    
 36                                                   
 37 PhysGeoImport::PhysGeoImport() :                  
 38     fFactor(1.),                                  
 39     fGeoName("")                                  
 40 {                                                 
 41     DefineMaterial();                             
 42 }                                                 
 43                                                   
 44 //....oooOO0OOooo........oooOO0OOooo........oo    
 45                                                   
 46 PhysGeoImport::PhysGeoImport(G4bool isVisu) :     
 47     fIsVisu(isVisu),                              
 48     fFactor(1.),                                  
 49     fGeoName("")                                  
 50 {                                                 
 51     if (fIsVisu) G4cout<<" For Checking Visual    
 52     DefineMaterial();                             
 53 }                                                 
 54                                                   
 55 //....oooOO0OOooo........oooOO0OOooo........oo    
 56                                                   
 57 G4LogicalVolume* PhysGeoImport::CreateLogicVol    
 58 {                                                 
 59     // The idea is to return a pointer to the     
 60     // Create a temporary material map to assi    
 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<    
 66     {                                             
 67         G4String matName("");                     
 68                                                   
 69         if(fMaterialVect[i]->GetName()=="adeni    
 70         else if(fMaterialVect[i]->GetName()=="    
 71         else if(fMaterialVect[i]->GetName()=="    
 72         else if(fMaterialVect[i]->GetName()=="    
 73         else if(fMaterialVect[i]->GetName()=="    
 74         else if(fMaterialVect[i]->GetName()=="    
 75         else if(fMaterialVect[i]->GetName()=="    
 76         else if(fMaterialVect[i]->GetName()=="    
 77         else if(fMaterialVect[i]->GetName()=="    
 78         else                                      
 79         {                                         
 80             matName = fMaterialVect[i]->GetNam    
 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,    
 95                                                   
 96     G4String boxNameLogic = fGeoName+"_logic";    
 97     G4LogicalVolume* box_logic = new G4Logical    
 98                                                   
 99     // sort the molecules                         
100     std::sort(fMolecules.begin(), fMolecules.e    
101                                                   
102     // Loop on all the parsed molecules           
103     for(G4int i=0, ie=fMolecules.size(); i<ie;    
104     {                                             
105         // Retrieve general molecule informati    
106         //                                        
107         G4String name = fMolecules[i].fName;      
108         G4String materialName = "water";// fMo    
109         G4double radius = fMolecules[i].fRadiu    
110         G4double waterRadius = fMolecules[i].f    
111         G4ThreeVector moleculePosition = fMole    
112         G4int copyNum = fMolecules[i].fCopyNum    
113                                                   
114         // Water hydration shell volume part      
115                                                   
116         G4Orb* moleculeWater_solid = 0;           
117         G4VSolid* moleculeWaterCut_solid = 0;     
118         G4LogicalVolume* moleculeWater_logic =    
119                                                   
120         // If water radius != 0 then we have a    
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+"_w    
129             moleculeWater_solid = new G4Orb(na    
130             moleculeWaterCut_solid = CreateCut    
131                                                   
132             G4String nameWaterLogic = name+"_w    
133             moleculeWater_logic = new G4Logica    
134                                                   
135             G4String nameWaterPhys = name+"_wa    
136             new G4PVPlacement(0, moleculePosit    
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].fNa    
146         molecule_solid = new G4Orb(nameSolid,     
147         moleculeCut_solid = CreateCutSolid(mol    
148                                                   
149         G4String nameLogic = name+"_logic";       
150         molecule_logic = new G4LogicalVolume(m    
151                                                   
152         // If there was a water hydration shel    
153         // placed within it and, thus, its rel    
154         if(waterRadius > (0 + tol)*nm)            
155         {                                         
156             G4ThreeVector position(0.,0.,0.);     
157             G4String namePhys = name+"_phys";     
158             new G4PVPlacement(0, position, mol    
159         }                                         
160         // If not, coordinates are those of th    
161         else                                      
162         {                                         
163             G4ThreeVector position = moleculeP    
164             G4String namePhys = name+"_phys";     
165             new G4PVPlacement(0, position, mol    
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........oo    
178                                                   
179 G4String PhysGeoImport::ParseFile(const G4Stri    
180 {                                                 
181     G4cout<<"Start parsing of "<<fileName<<G4e    
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    
196         G4Exception("PhysGeoImport::ParseFile(    
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 i    
209         if(line.empty()) continue; // skip the    
210                                                   
211         // Data string stream                     
212         std::istringstream issLine(line);         
213                                                   
214         // String to determine the first lette    
215         G4String firstItem;                       
216                                                   
217         // Put the first letter/word within th    
218         issLine >> firstItem;                     
219                                                   
220         // Check first letter to determine if     
221         if(firstItem=="#") continue; // skip t    
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] = waterRadiu    
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 != preHistoneCp    
279                     NoHistone ++;                 
280                     preHistoneCpNo = copyNumbe    
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    
302                                                   
303             fMolecules.push_back(molecule);       
304         }                                         
305         else                                      
306         {                                         
307             // Geant4 exception                   
308             G4String msg = firstItem+" is not     
309             G4Exception("PhysGeoImport::ParseF    
310         }                                         
311     }                                             
312                                                   
313     // Close the file once the reading is done    
314     file.close();                                 
315                                                   
316     G4cout<<"End parsing of "<<fileName<<G4end    
317     fVoxelNbHistoneMap.insert({fGeoName,NoHist    
318     fVoxelNbBpMap.insert({fGeoName,NoBp});        
319     return fGeoName;                              
320 }                                                 
321                                                   
322 //....oooOO0OOooo........oooOO0OOooo........oo    
323                                                   
324 G4VSolid* PhysGeoImport::CreateCutSolid(G4Orb     
325                                                   
326                                                   
327                                                   
328 {                                                 
329     // The idea behing this method is to cut o    
330     // If a reference and a target volumes are    
331     // The reference is already selected when     
332                                                   
333     // Use the tiny space to differentiate the    
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 = mol    
346     else radiusRef = molRef.fRadiusWater;         
347                                                   
348     // Reference volume position                  
349     G4ThreeVector posRef = molRef.fPosition;      
350     // Before looping on all the volumes we ch    
351     //reference volume overlaps with its conta    
352     if(std::abs(posRef.x() ) + radiusRef > fSi    
353             || std::abs(posRef.y() ) + radiusR    
354             || std::abs(posRef.z() ) + radiusR    
355     {                                             
356         // If we enter here, then the referenc    
357         // Box used to cut                        
358         G4Box* solidBox = new G4Box("solid_box    
359         G4ThreeVector posBox;                     
360                                                   
361         // Create a dummy rotation matrix         
362         G4RotationMatrix *rotMat = new G4Rotat    
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 )     
369         {                                         
370            posBox = -posRef +  G4ThreeVector(0    
371                                                   
372             // If the volume is cutted for the    
373             if(!isCutted)                         
374             {                                     
375                 solidCut = new G4SubtractionSo    
376                 isCutted = true;                  
377             }                                     
378             // For the other times                
379             else solidCut = new G4SubtractionS    
380         }                                         
381                                                   
382         // Down                                   
383         if(std::abs( posRef.y() - radiusRef )     
384         {                                         
385             posBox = -posRef + G4ThreeVector(0    
386                                                   
387             // If the volume is cutted for the    
388             if(!isCutted)                         
389             {                                     
390                 solidCut = new G4SubtractionSo    
391             isCutted = true;                      
392         }                                         
393             // For the other times                
394             else solidCut = new G4SubtractionS    
395         }                                         
396                                                   
397         // Left                                   
398         if(std::abs( posRef.x() + radiusRef )     
399         {                                         
400             posBox = -posRef + G4ThreeVector(f    
401                                                   
402             // If the volume is cutted for the    
403             if(!isCutted)                         
404             {                                     
405                 solidCut = new G4SubtractionSo    
406                 isCutted = true;                  
407             }                                     
408             // For the other times                
409             else solidCut = new G4SubtractionS    
410         }                                         
411                                                   
412         // Right                                  
413         if(std::abs( posRef.x() - radiusRef )     
414         {                                         
415             posBox = -posRef + G4ThreeVector(-    
416                                                   
417             // If the volume is cutted for the    
418             if(!isCutted)                         
419             {                                     
420                 solidCut = new G4SubtractionSo    
421                 isCutted = true;                  
422             }                                     
423             // For the other times                
424             else solidCut = new G4SubtractionS    
425         }                                         
426                                                   
427         // Forward                                
428         if(std::abs( posRef.z() + radiusRef )     
429         {                                         
430             posBox = -posRef + G4ThreeVector(0    
431                                                   
432             // If the volume is cutted for the    
433             if(!isCutted)                         
434             {                                     
435                 solidCut = new G4SubtractionSo    
436                 isCutted = true;                  
437             }                                     
438             // For the other times                
439             else solidCut = new G4SubtractionS    
440         }                                         
441                                                   
442         // Backward                               
443         if(std::abs( posRef.z() - radiusRef )     
444         {                                         
445             posBox = -posRef + G4ThreeVector(0    
446                                                   
447             // If the volume is cutted for the    
448             if(!isCutted)                         
449             {                                     
450                 solidCut = new G4SubtractionSo    
451                 isCutted = true;                  
452             }                                     
453             // For the other times                
454             else solidCut = new G4SubtractionS    
455         }                                         
456     }                                             
457                                                   
458     // Look the other volumes of the voxel        
459     // Loop on all the target volumes (other v    
460     for(G4int i=0, ie=molList.size(); i<ie; ++    
461     {                                             
462         G4ThreeVector posTar = molList[i].fPos    
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 info    
477         G4double radiusTar;                       
478         if(molList[i].fRadiusWater==0) radiusT    
479         else radiusTar = molList[i].fRadiusWat    
480                                                   
481         // Compute the distance reference-targ    
482         G4double distance = std::abs( (posRef     
483                                                   
484         // Use the distance to check if the cu    
485         // This can only happen once per loop.    
486         if(distance==0 && !isOurVol)              
487         {                                         
488             // Target volume is also reference    
489                                                   
490             // Set the flag                       
491             isOurVol = true;                      
492                                                   
493             // Next iteration                     
494             continue;                             
495         }                                         
496         // If the condition is correct more th    
497         else if(distance == 0 && isOurVol)        
498         {                                         
499             G4ExceptionDescription desmsg;        
500             desmsg << "DetectorConstruction::C    
501             G4Exception("ChemGeoImport::BuildG    
502         }                                         
503                                                   
504         // If the volumes are differents then     
505         // close enough to overlap and, thus,     
506         else if(distance <= radiusRef+radiusTa    
507         {                                         
508             // Volumes are close enough, there    
509                                                   
510             // Box used to cut                    
511             G4Box* solidBox = new G4Box("solid    
512                                                   
513             // This part is tricky.               
514             // The goal is to calculate the po    
515                                                   
516             // diff vector to from ref to tar     
517             G4ThreeVector diff = posTar - posR    
518                                                   
519             // Find the intersection point and    
520             G4double d = (std::pow(radiusRef,2    
521                         (2*distance) + solidBo    
522                                                   
523             // If we are in another volume we     
524             //differentiate without ambiguitie    
525             // (may not be necessary)             
526             if(in) d -= 2*tinySpace;              
527                                                   
528             // Position of the box in order to    
529             // "* ( diff/diff.getR() )" is nec    
530             G4ThreeVector pos = d *( diff/diff    
531                                                   
532             // Get the rotation angles because    
533             // to give the right "cut angle".     
534             G4double phi = std::acos(pos.getZ(    
535             G4double theta = std::acos( pos.ge    
536                                                   
537             if(pos.getY()<0) theta = -theta;      
538                                                   
539             G4ThreeVector rotAxisForPhi(1*fFac    
540             rotAxisForPhi.rotateZ(theta+pi/2);    
541                                                   
542             // Create the rotation matrix         
543             G4RotationMatrix *rotMat = new G4R    
544             rotMat->rotate(-phi, rotAxisForPhi    
545                                                   
546             // Rotate it again                    
547             G4ThreeVector rotZAxis(0.,0.,1*fFa    
548             rotMat->rotate(theta, rotZAxis);      
549                                                   
550             // If the volume is cutted for the    
551             if(!isCutted) solidCut = new G4Sub    
552             solidBox, rotMat, pos);               
553                                                   
554             // For the other times                
555             else solidCut = new G4SubtractionS    
556                                                   
557             // Set the cut flag                   
558             isCutted = true;                      
559         }                                         
560     }                                             
561                                                   
562     // If there was at least one cut then we r    
563     if(isCutted) return solidCut;                 
564                                                   
565     // Otherwise, we return the original volum    
566     else return solidOrbRef;                      
567                                                   
568 }                                                 
569                                                   
570 //....oooOO0OOooo........oooOO0OOooo........oo    
571                                                   
572 void PhysGeoImport::DefineMaterial()              
573 {                                                 
574     G4NistManager * man = G4NistManager::Insta    
575     fpWater = man->FindOrBuildMaterial("G4_WAT    
576     fMaterialVect.push_back(fpWater);             
577     fVacuum = man->FindOrBuildMaterial("G4_Gal    
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="Carbo    
586     a = 1.00794*g/mole;                           
587     G4Element* elH  = new G4Element(name="Hydr    
588     a = 15.9994*g/mole;                           
589     G4Element* elO  = new G4Element(name="Oxyg    
590     a = 14.00674*g/mole;                          
591     G4Element* elN  = new G4Element(name="Nitr    
592     a = 30.973762*g/mole;                         
593     G4Element* elP  = new G4Element(name="Phos    
594                                                   
595     density = 1.*g/cm3;                           
596                                                   
597     fTHF = new G4Material("THF", density, nCom    
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, nCompo    
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, nCompo    
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", d    
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.0    
638     G4double VBackbone = 1.823912/20. * std::p    
639     G4double densityBaTHF = (MBackbone/(g/mole    
640                                                   
641     fDeoxyribose = new G4Material("backbone_TH    
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    
648     VBackbone = 1.19114/20. * std::pow(1.e-9,     
649     G4double densityBaTMP = (MBackbone/(g/mole    
650                                                   
651     fPhosphate = new G4Material("backbone_TMP"    
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.00    
658     G4double VBase = 1.8527205/20. * std::pow(    
659     G4double densityCy = (MCytosine/(g/mole) *    
660                                                   
661     fCytosine_PY = new G4Material("cytosine_PY    
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 +    
669     VBase = 1.8527205/20. * std::pow(1.e-9, 3.    
670     densityCy = (MThy/(g/mole) * 1. )/(6.022e2    
671                                                   
672     fThymine_PY = new G4Material("thymine_PY",    
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.007    
680     VBase = 1.8527205/20. * std::pow(1.e-9, 3.    
681     G4double densityAd = (MAdenine/(g/mole) *     
682                                                   
683     fAdenine_PU = new G4Material("adenine_PU",    
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.*1    
690     VBase = 1.8527205/20. * std::pow(1.e-9, 3.    
691     densityAd = (MAdenine/(g/mole) * 1. )/(6.0    
692                                                   
693     fGuanine_PU = new G4Material("guanine_PU",    
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("homogen    
701     fHomogeneous_dna->AddMaterial(fpWater, 0.3    
702     fHomogeneous_dna->AddMaterial(fDeoxyribose    
703     fHomogeneous_dna->AddMaterial(fPhosphate,     
704     fHomogeneous_dna->AddMaterial(fAdenine_PU,    
705     fHomogeneous_dna->AddMaterial(fGuanine_PU,    
706     fHomogeneous_dna->AddMaterial(fCytosine_PY    
707     fHomogeneous_dna->AddMaterial(fThymine_PY,    
708                                                   
709     fMaterialVect.push_back(fHomogeneous_dna);    
710 }                                                 
711                                                   
712 //....oooOO0OOooo........oooOO0OOooo........oo    
713                                                   
714 G4LogicalVolume* PhysGeoImport::CreateNucleusL    
715 {                                                 
716     // Here, we parse the nucleus file and bui    
717                                                   
718     G4cout<<"Start the nucleus creation..."<<G    
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    
734         G4Exception("PhysGeoImport::ParseFile(    
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) && !    
742     {                                             
743         // Check the line to determine if it i    
744         if(line.empty()) continue; // skip the    
745                                                   
746         // Data string stream                     
747         std::istringstream issLine(line);         
748                                                   
749         // String to determine the first lette    
750         G4String firstItem;                       
751                                                   
752         // Put the first letter/word within th    
753         issLine >> firstItem;                     
754                                                   
755         // Check first letter to determine if     
756         if(firstItem=="#") continue; // skip t    
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 >> s    
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    
784                 nucleusLogic = new G4LogicalVo    
785                 G4double r3 = semiX*semiY*semi    
786                 fNucleusVolume = 4.0*pi*r3/3.0    
787             } else if (fNucleusType == "Ellipt    
788               G4float semiX, semiY, semiZ;        
789                                                   
790                 issLine >> semiX >> semiY >> s    
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 G4Elliptica    
801                 nucleusLogic = new G4LogicalVo    
802                 G4double r3 = 2*semiX*semiY*se    
803                 fNucleusVolume = pi*r3;           
804             } else if (fNucleusType == "Spheri    
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(fNucl    
812                 nucleusLogic = new G4LogicalVo    
813                 G4double r3 = radius*radius*ra    
814                 fNucleusVolume = 4.0*pi*r3/3.0    
815             } else {                              
816                 G4String msg =fNucleusType+" i    
817                 G4Exception("PhysGeoImport::Cr    
818                 "Geo_InputFile_Unknown_Nucleus    
819             }                                     
820                                                   
821             isNucleusBuild = true;                
822         }                                         
823     }                                             
824                                                   
825     nucleusFile.close();                          
826                                                   
827     if(!isNucleusBuild)                           
828     {                                             
829         G4String msg = "Nucleus data were not     
830         G4Exception("PhysGeoImport::CreateNucl    
831         "Geo_InputFile_Unknown_Nucleus", Fatal    
832     }                                             
833                                                   
834     return nucleusLogic;                          
835 }                                                 
836                                                   
837 //....oooOO0OOooo........oooOO0OOooo........oo    
838                                                   
839 std::vector<Voxel>* PhysGeoImport::CreateVoxel    
840 {                                                 
841     fTotalNbBpPlacedInGeo = 0;                    
842     fTotalNbHistonePlacedInGeo = 0;               
843     G4cout<<"Start the voxel data generation..    
844                                                   
845     std::vector<Voxel>* voxels = new std::vect    
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    
861         G4Exception("PhysGeoImport::ParseFile(    
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::VoxelT    
871         if (voxt == Voxel::Straight || voxt ==    
872             voxt == Voxel::Up || voxt == Voxel    
873         else if (voxt == Voxel::Straight2 || v    
874             voxt == Voxel::Up2 || voxt == Voxe    
875         else return ChromatinType::fUnspecifie    
876     };                                            
877     // Read the file line per line                
878     while(std::getline(nucleusFile, line) )       
879     {                                             
880         // Check the line to determine if it i    
881         if(line.empty()) continue; // skip the    
882                                                   
883         // Data string stream                     
884         std::istringstream issLine(line);         
885                                                   
886         // String to determine the first lette    
887         G4String firstItem;                       
888                                                   
889         // Put the first letter/word within th    
890         issLine >> firstItem;                     
891                                                   
892         // Check first letter to determine if     
893         if(firstItem=="#") continue; // skip t    
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 G4Rota    
951             G4ThreeVector(yx,yy,yz),G4ThreeVec    
952                                                   
953             Voxel::VoxelType type = Voxel::Oth    
954                                                   
955             if(name=="voxelStraight" || name==    
956             {                                     
957                 type = Voxel::Straight;           
958             }                                     
959             else if(name=="voxelUp" || name=="    
960             {                                     
961                 type = Voxel::Up;                 
962             }                                     
963             else if(name=="voxelDown" || name=    
964             {                                     
965                 type = Voxel::Down;               
966             }                                     
967             else if(name=="voxelRight" || name    
968             {                                     
969                 type = Voxel::Right;              
970             }                                     
971             else if(name=="voxelLeft" || name=    
972             {                                     
973                 type = Voxel::Left;               
974             }                                     
975             else if(name=="voxelStraight2" ||     
976             {                                     
977                 type = Voxel::Straight2;          
978             }                                     
979             else if(name=="voxelUp2" || name==    
980             {                                     
981                 type = Voxel::Up2;                
982             }                                     
983             else if(name=="voxelDown2" || name    
984             {                                     
985                 type = Voxel::Down2;              
986             }                                     
987             else if(name=="voxelRight2" || nam    
988             {                                     
989                 type = Voxel::Right2;             
990             }                                     
991             else if(name=="voxelLeft2" || name    
992             {                                     
993                 type = Voxel::Left2;              
994             }                                     
995             else                                  
996             {                                     
997                 G4ExceptionDescription msg ;      
998                 msg << "The voxel named "<<nam    
999                 G4Exception("PhysGeoImport::Cr    
1000             }                                    
1001                                                  
1002             voxels->push_back(Voxel(voxelCopy    
1003             G4ThreeVector(x,y,z), rot) );        
1004             fTotalNbBpPlacedInGeo += fVoxelNb    
1005             fTotalNbHistonePlacedInGeo += fVo    
1006             voxelCopyNumber++;                   
1007             auto chromatintype =  whatchromat    
1008             if (fChromatinTypeCount.find(chro    
1009                 fChromatinTypeCount.insert({c    
1010             else fChromatinTypeCount[chromati    
1011         }                                        
1012     }                                            
1013                                                  
1014     nucleusFile.close();                         
1015                                                  
1016     G4cout<<"End the voxel data generation"<<    
1017                                                  
1018     return voxels;                               
1019 }                                                
1020                                                  
1021 //....oooOO0OOooo........oooOO0OOooo........o