Changeset 1337 for trunk/examples/extended/medical/DICOM/src

Timestamp:

Sep 30, 2010, 2:47:17 PM (14 years ago)

Author:

garnier

Message:

tag geant4.9.4 beta 1 + modifs locales

Location:

trunk/examples/extended/medical/DICOM/src

Files:

• DicomDetectorConstruction.cc (modified) (4 diffs)
• DicomEventAction.cc (modified) (2 diffs)
• DicomHandler.cc (modified) (12 diffs)
• DicomNestedPhantomParameterisation.cc (modified) (4 diffs)
• DicomPatientZSliceHeader.cc (modified) (1 diff)
• DicomPhantomParameterisationColour.cc (modified) (1 diff)
• DicomPhysicsList.cc (modified) (2 diffs)
• DicomPrimaryGeneratorAction.cc (modified) (1 diff)

trunk/examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc

@@ -314 +314 @@
   G4double densityDiff = 0.1; 
   std::map<G4int,G4double>  fDensityDiffs; // to be able to use a different densityDiff for each material
-  for( size_t ii = 0; ii < fOriginalMaterials.size(); ii++ ){
+  for( unsigned int ii = 0; ii < fOriginalMaterials.size(); ii++ ){
     fDensityDiffs[ii] = densityDiff; //currently all materials with same difference
   }
@@ -327 +327 @@
   //--- If first slice, initiliaze fMateIDs
   if( fZSliceHeaders.size() == 1 ) {
+    //fMateIDs = new unsigned int[fNoFiles*nVoxels];
     fMateIDs = new size_t[fNoFiles*nVoxels];
-  }
-
-  size_t mateID;
+
+  }
+
+  unsigned int mateID;
   G4int voxelCopyNo = (fZSliceHeaders.size()-1)*nVoxels; // number of voxels from previously read slices
   for( G4int ii = 0; ii < nVoxels; ii++, voxelCopyNo++ ){
@@ -353 +355 @@
     
     //-- Look if a material with this name is already created (because a previous voxel was already in this density bin)
-    size_t im;
+    unsigned int im;
     for( im = 0; im < fMaterials.size(); im++ ){
       if( fMaterials[im]->GetName() == newMateName ) {
@@ -377 +379 @@
   //----- Images must have the same dimension ... 
   fZSliceHeaderMerged = new DicomPatientZSliceHeader( *fZSliceHeaders[0] );
-  for( size_t ii = 1; ii < fZSliceHeaders.size(); ii++ ) {
+  for( unsigned int ii = 1; ii < fZSliceHeaders.size(); ii++ ) {
     *fZSliceHeaderMerged += *fZSliceHeaders[ii];
   };

trunk/examples/extended/medical/DICOM/src/DicomEventAction.cc

@@ -42 +42 @@
 #include "DicomEventAction.hh"
 #include "G4Event.hh"
-#include "G4TrajectoryContainer.hh"
-#include "G4Trajectory.hh"
-#include "G4VVisManager.hh"
 
 DicomEventAction::DicomEventAction():drawFlag("all")
@@ -55 +52 @@
 {
   G4cout << "EV: " << evt->GetEventID() << G4endl;
- }
-
-void DicomEventAction::EndOfEventAction(const G4Event* evt)
-{
-  // extract the trajectories and draw them
-
-  if (G4VVisManager::GetConcreteInstance())
-    {
-      G4TrajectoryContainer * trajectoryContainer = evt->GetTrajectoryContainer();
-      G4int nTrajectories = 0;
-      if ( trajectoryContainer )
-        nTrajectories = trajectoryContainer->entries();
-
-      for ( G4int i = 0;i < nTrajectories; i++ )
-        {
-          G4Trajectory* trj = (G4Trajectory*)((*(evt->GetTrajectoryContainer()))[i]);
-          if (drawFlag == "all")
-            trj->DrawTrajectory(50);
-          else if ((drawFlag == "charged")&&(trj->GetCharge() != 0.))
-            trj->DrawTrajectory(50);
-          else if ((drawFlag == "neutral")&&(trj->GetCharge() == 0.))
-            trj->DrawTrajectory(50);
-        }
-    }
 }
 
+void DicomEventAction::EndOfEventAction(const G4Event*)
+{
+}

trunk/examples/extended/medical/DICOM/src/DicomHandler.cc

@@ -61 +61 @@
 
 DicomHandler::DicomHandler()
-    : DATABUFFSIZE(8192), LINEBUFFSIZE(128), FILENAMESIZE(512),
+    : DATABUFFSIZE(8192), LINEBUFFSIZE(5020), FILENAMESIZE(512),
       compression(0), nFiles(0), rows(0), columns(0),
       bitAllocated(0), maxPixelValue(0), minPixelValue(0),
@@ -98 +98 @@
     short readElementId;  // identify a particular type information
     short elementLength2; // deal with element length in 2 bytes
+    //unsigned int elementLength4; // deal with element length in 4 bytes
     G4int elementLength4; // deal with element length in 4 bytes
 
     char * data = new char[DATABUFFSIZE];
 
+ 
     // Read information up to the pixel data
     while(true) {
@@ -117 +119 @@
         // Creating a tag to be identified afterward
         G4int tagDictionary = readGroupId*0x10000 + readElementId;
-
-
-        // VR or element length
+        
+        // beginning of the pixels
+            if(tagDictionary == 0x7FE00010) break;
+
+      // VR or element length
         std::fread(buffer,2,1,dicom);
         GetValue(buffer, elementLength2);
-
-        // If value representation (VR) is OB, OW, SQ, UN, 
+          
+        // If value representation (VR) is OB, OW, SQ, UN, added OF and UT 
         //the next length is 32 bits
         if((elementLength2 == 0x424f ||  // "OB"
             elementLength2 == 0x574f ||  // "OW"
-            elementLength2 == 0x5153 ||  // "SQ"
+            elementLength2 == 0x464f ||  // "OF"
+            elementLength2 == 0x5455 ||  // "UT"
+            elementLength2 == 0x5153 || //  "SQ"
             elementLength2 == 0x4e55) && // "UN"
            !implicitEndian ) {           // explicit VR
@@ -136 +142 @@
             std::fread(buffer, 4, 1, dicom);
             GetValue(buffer, elementLength4);
-
-            // beginning of the pixels
-            if(tagDictionary == 0x7FE00010) break;
-
+            
+            if(elementLength2 == 0x5153)
+            {
+             if(elementLength4 == 0xFFFFFFFF)
+              read_undefined_nested( dicom );
+             else{
+              if(read_defined_nested( dicom, elementLength4 )==0){
+               G4cerr << "Reading defined nested functin failed!" << G4endl;
+               exit(-10);              }
+              }
+            } else  { 
             // Reading the information with data
             std::fread(data, elementLength4,1,dicom);
-
-
-        } else { // length is 16 bits :
-
-            if(!implicitEndian || readGroupId == 2) {
-                // element length (2 bytes)
-                std::fread(buffer, 2, 1, dicom);
-                GetValue(buffer, elementLength2);
-                elementLength4 = elementLength2;
-
-            } else { 
-                // element length (4 bytes)
-                if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
-                    G4cerr << "[DicomHandler] fseek failed" << G4endl;
-                    exit(-10);
-                }
-                std::fread(buffer, 4, 1, dicom);
-                GetValue(buffer, elementLength4);
-            }
-
-            // beginning of the pixels
-            if(tagDictionary == 0x7FE00010) break;
-
-            std::fread(data, elementLength4, 1, dicom);
+            }  
+
+                
+        }  else { 
+
+                if(!implicitEndian || readGroupId == 2) {  //  explicit with VR different than previous ones
+                
+                  //G4cout << "Reading  DICOM files with Explicit VR"<< G4endl;
+                  // element length (2 bytes)
+                  std::fread(buffer, 2, 1, dicom);
+                  GetValue(buffer, elementLength2);
+                  elementLength4 = elementLength2;
+                  
+                  
+                  
+                  std::fread(data, elementLength4, 1, dicom);
+                
+                } else {                                 // Implicit VR
+
+                  //G4cout << "Reading  DICOM files with Implicit VR"<< G4endl;
+   
+                  // element length (4 bytes)
+                  if(std::fseek(dicom, -2, SEEK_CUR) != 0) {
+                      G4cerr << "[DicomHandler] fseek failed" << G4endl;
+                      exit(-10);}
+
+                  std::fread(buffer, 4, 1, dicom);
+                  GetValue(buffer, elementLength4);
+
+                  //G4cout <<  std::hex<< elementLength4 << G4endl;
+              
+                  if(elementLength4 == 0xFFFFFFFF) 
+                   read_undefined_nested(dicom);
+                  else {
+                  
+                  
+                  std::fread(data, elementLength4, 1, dicom);
+                
+                 } 
+                      
+               } 
         }
 
@@ -179 +209 @@
     G4String fnameG4DCM = G4String(filename2) + ".g4dcm";
     foutG4DCM.open(fnameG4DCM);
-    G4cout << " opened fnameG4DCM file " << fnameG4DCM << G4endl;
+    G4cout << "### Writing of " << fnameG4DCM << " ### " << G4endl;
 
     foutG4DCM << fMaterialIndices.size() << G4endl;
     //--- Write materials
-    size_t ii = 0;
-    std::map<G4double,G4String>::const_iterator ite;
+    unsigned int ii = 0;
+    std::map<G4float,G4String>::const_iterator ite;
     for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ite++, ii++ ){
       foutG4DCM << ii << " " << (*ite).second << G4endl;
@@ -194 +224 @@
     foutG4DCM << -pixelSpacingY*columns/2 << " " << pixelSpacingY*columns/2 << G4endl;
     foutG4DCM << sliceLocation-sliceThickness/2. << " " << sliceLocation+sliceThickness/2. << G4endl;
+    //    foutG4DCM << compression << G4endl;
     
     ReadData( dicom, filename2 );
@@ -431 +462 @@
 void DicomHandler::ReadMaterialIndices( std::ifstream& finData)
 {
-  size_t nMate;
+  unsigned int nMate;
   G4String mateName;
-  G4double densityMax;
+  G4float densityMax;
   finData >> nMate;
   if( finData.eof() ) return;
 
   G4cout << " ReadMaterialIndices " << nMate << G4endl;
-  for( size_t ii = 0; ii < nMate; ii++ ){
+  for( unsigned int ii = 0; ii < nMate; ii++ ){
     finData >> mateName >> densityMax;
     fMaterialIndices[densityMax] = mateName;
@@ -446 +477 @@
 }
 
-size_t DicomHandler::GetMaterialIndex( G4double density )
+unsigned int DicomHandler::GetMaterialIndex( G4float density )
 {
-  size_t mateID;
-  std::map<G4double,G4String>::reverse_iterator ite;
+  std::map<G4float,G4String>::reverse_iterator ite;
   G4int ii = fMaterialIndices.size();
   for( ite = fMaterialIndices.rbegin(); ite != fMaterialIndices.rend(); ite++, ii-- ) {
@@ -508 +538 @@
 
     // Creation of .g4 files wich contains averaged density data
-
     char * nameProcessed = new char[FILENAMESIZE];
-    FILE* processed;
-
-    std::sprintf(nameProcessed,"%s.g4",filename2);
-    processed = std::fopen(nameProcessed,"w+b");
+    FILE* fileOut;
+
+    std::sprintf(nameProcessed,"%s.g4dcmb",filename2);
+    fileOut = std::fopen(nameProcessed,"w+b");
     std::printf("### Writing of %s ###\n",nameProcessed);
 
-    std::fwrite(&rows, 2, 1, processed);
-    std::fwrite(&columns, 2, 1, processed);
-    std::fwrite(&pixelSpacingX, 8, 1, processed);
-    std::fwrite(&pixelSpacingY, 8, 1, processed);
-    std::fwrite(&sliceThickness, 8, 1, processed);
-    std::fwrite(&sliceLocation, 8, 1, processed);
-    std::fwrite(&compression, 2, 1, processed);
+    unsigned int nMate = fMaterialIndices.size();
+    std::fwrite(&nMate, sizeof(unsigned int), 1, fileOut);
+    //--- Write materials
+    std::map<G4float,G4String>::const_iterator ite;
+    for( ite = fMaterialIndices.begin(); ite != fMaterialIndices.end(); ite++ ){
+      G4String mateName = (*ite).second;
+      for( G4int ii = (*ite).second.length(); ii < 40; ii++ ) {
+        mateName += " ";
+      }  //mateName = const_cast<char*>(((*ite).second).c_str());
+
+      const char* mateNameC = mateName.c_str();
+      std::fwrite(mateNameC, sizeof(char),40, fileOut);
+    }
+
+    unsigned int rowsC = rows/compression;
+    unsigned int columnsC = columns/compression;
+    unsigned int planesC = 1;
+    G4float pixelLocationXM = -pixelSpacingX*rows/2.;
+    G4float pixelLocationXP = pixelSpacingX*rows/2.;
+    G4float pixelLocationYM = -pixelSpacingY*rows/2.;
+    G4float pixelLocationYP = pixelSpacingY*rows/2.;
+    G4float sliceLocationZM = sliceLocation-sliceThickness/2.;
+    G4float sliceLocationZP = sliceLocation+sliceThickness/2.;
+    //--- Write number of voxels (assume only one voxel in Z)
+    std::fwrite(&rowsC, sizeof(unsigned int), 1, fileOut);
+    std::fwrite(&columnsC, sizeof(unsigned int), 1, fileOut);
+    std::fwrite(&planesC, sizeof(unsigned int), 1, fileOut);
+    //--- Write minimum and maximum extensions
+    std::fwrite(&pixelLocationXM, sizeof(G4float), 1, fileOut);
+    std::fwrite(&pixelLocationXP, sizeof(G4float), 1, fileOut);
+    std::fwrite(&pixelLocationYM, sizeof(G4float), 1, fileOut);
+    std::fwrite(&pixelLocationYP, sizeof(G4float), 1, fileOut);
+    std::fwrite(&sliceLocationZM, sizeof(G4float), 1, fileOut);
+    std::fwrite(&sliceLocationZP, sizeof(G4float), 1, fileOut);
+    //    std::fwrite(&compression, sizeof(unsigned int), 1, fileOut);
 
     std::printf("%8i   %8i\n",rows,columns);
@@ -532 +589 @@
     G4int compSize = compression;
     G4int mean;
-    G4double density;
+    G4float density;
     G4bool overflow = false;
     G4int cpt=1;
 
+    //----- Write index of material for each pixel
     if(compSize == 1) { // no compression: each pixel has a density value)
-        for( G4int ww = 0; ww < rows; ww++) {
-            for( G4int xx = 0; xx < columns; xx++) {
-                mean = tab[ww][xx];
-                density = Pixel2density(mean);
-                std::fwrite(&density, sizeof(G4double), 1, processed);
-            }
-        }
+      for( G4int ww = 0; ww < rows; ww++) {
+        for( G4int xx = 0; xx < columns; xx++) {
+          mean = tab[ww][xx];
+          density = Pixel2density(mean);
+          unsigned int mateID = GetMaterialIndex( density );
+          std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
+        }
+      }
 
     } else {
-        // density value is the average of a square region of
-        // compression*compression pixels
-        for(G4int ww = 0; ww < rows ;ww += compSize ) {
-            for(G4int xx = 0; xx < columns ;xx +=compSize ) {
-                overflow = false;
-                mean = 0;
-                for(int sumx = 0; sumx < compSize; sumx++) {
-                    for(int sumy = 0; sumy < compSize; sumy++) {
-                        if(ww+sumy >= rows || xx+sumx >= columns) overflow = true;
-                        mean += tab[ww+sumy][xx+sumx];
-                    }
-                    if(overflow) break;
-                }
-                mean /= compSize*compSize;
-                cpt = 1;
-
-                if(!overflow) {
-                    G4double density = Pixel2density(mean);
-                    std::fwrite(&density, sizeof(G4double), 1, processed);
-                }
+      // density value is the average of a square region of
+      // compression*compression pixels
+      for(G4int ww = 0; ww < rows ;ww += compSize ) {
+        for(G4int xx = 0; xx < columns ;xx +=compSize ) {
+          overflow = false;
+          mean = 0;
+          for(int sumx = 0; sumx < compSize; sumx++) {
+            for(int sumy = 0; sumy < compSize; sumy++) {
+              if(ww+sumy >= rows || xx+sumx >= columns) overflow = true;
+              mean += tab[ww+sumy][xx+sumx];
             }
-
-        }
-    }
-    std::fclose(processed);
-
+            if(overflow) break;
+          }
+          mean /= compSize*compSize;
+          cpt = 1;
+          
+          if(!overflow) {
+            density = Pixel2density(mean);
+            unsigned int mateID = GetMaterialIndex( density );
+            std::fwrite(&mateID, sizeof(unsigned int), 1, fileOut);
+          }
+        }
+        
+      }
+    }
+
+    //----- Write density for each pixel
+    if(compSize == 1) { // no compression: each pixel has a density value)
+      for( G4int ww = 0; ww < rows; ww++) {
+        for( G4int xx = 0; xx < columns; xx++) {
+          mean = tab[ww][xx];
+          density = Pixel2density(mean);
+          std::fwrite(&density, sizeof(G4float), 1, fileOut);
+        }
+      }
+      
+    } else {
+      // density value is the average of a square region of
+      // compression*compression pixels
+      for(G4int ww = 0; ww < rows ;ww += compSize ) {
+        for(G4int xx = 0; xx < columns ;xx +=compSize ) {
+          overflow = false;
+          mean = 0;
+          for(int sumx = 0; sumx < compSize; sumx++) {
+            for(int sumy = 0; sumy < compSize; sumy++) {
+              if(ww+sumy >= rows || xx+sumx >= columns) overflow = true;
+              mean += tab[ww+sumy][xx+sumx];
+            }
+            if(overflow) break;
+          }
+          mean /= compSize*compSize;
+          cpt = 1;
+          
+          if(!overflow) {
+            density = Pixel2density(mean);
+            std::fwrite(&density, sizeof(G4float), 1, fileOut);
+          }
+        }
+        
+      }
+    }
+    
+    std::fclose(fileOut);
+    
     delete [] nameProcessed;
 
@@ -595 +692 @@
 */
 
-G4double DicomHandler::Pixel2density(G4int pixel)
+G4float DicomHandler::Pixel2density(G4int pixel)
 {
-    G4double density = -1.;
+    G4float density = -1.;
     G4int nbrequali = 0;
     G4double deltaCT = 0;
@@ -752 +849 @@
     _rval = *(Type *)_val;
 }
+
+G4int DicomHandler::read_defined_nested(FILE * nested,G4int SQ_Length)
+{
+  //      VARIABLES
+  unsigned short item_GroupNumber;
+  unsigned short item_ElementNumber;
+  G4int item_Length;
+  G4int items_array_length=0;
+  char * buffer= new char[LINEBUFFSIZE];
+  
+
+  while(items_array_length < SQ_Length)
+  {
+   std::fread(buffer, 2, 1, nested);
+   GetValue(buffer, item_GroupNumber);
+   
+   std::fread(buffer, 2, 1, nested);
+   GetValue(buffer, item_ElementNumber);
+   
+   std::fread(buffer, 4, 1, nested);
+   GetValue(buffer, item_Length);
+   
+   std::fread(buffer, item_Length, 1, nested);
+   
+   items_array_length= items_array_length+8+item_Length;
+  }
+ 
+  delete []buffer;
+  
+  if(SQ_Length>items_array_length)
+   return 0;
+  else
+   return 1;
+
+}
+
+void DicomHandler::read_undefined_nested(FILE * nested)
+{
+  //      VARIABLES
+  unsigned short item_GroupNumber;
+  unsigned short item_ElementNumber;
+  G4int item_Length;
+  char * buffer= new char[LINEBUFFSIZE];
+  
+
+  do
+  {
+   std::fread(buffer, 2, 1, nested);
+   GetValue(buffer, item_GroupNumber);
+   
+   std::fread(buffer, 2, 1, nested);
+   GetValue(buffer, item_ElementNumber);
+   
+   std::fread(buffer, 4, 1, nested);
+   GetValue(buffer, item_Length);
+   
+   if(item_Length!=0xffffffff)
+    std::fread(buffer, item_Length, 1, nested);
+   else
+    read_undefined_item(nested);
+   
+   
+  }while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE0DD || item_Length!=0);
+
+  delete [] buffer;
+
+}
+
+void DicomHandler::read_undefined_item(FILE * nested)
+{
+  //      VARIABLES
+ unsigned short item_GroupNumber;
+ unsigned short item_ElementNumber;
+ G4int item_Length;
+ char *buffer= new char[LINEBUFFSIZE];
+ 
+ do
+ {
+  std::fread(buffer, 2, 1, nested);
+  GetValue(buffer, item_GroupNumber);
+   
+  std::fread(buffer, 2, 1, nested);
+  GetValue(buffer, item_ElementNumber);
+   
+  std::fread(buffer, 4, 1, nested);
+  GetValue(buffer, item_Length);
+
+
+  if(item_Length!=0)
+   std::fread(buffer,item_Length,1,nested);
+
+ }
+ while(item_GroupNumber!=0xFFFE || item_ElementNumber!=0xE00D || item_Length!=0);
+ 
+ delete []buffer;
+ 
+}

trunk/examples/extended/medical/DICOM/src/DicomNestedPhantomParameterisation.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: DicomNestedPhantomParameterisation.cc,v 1.5 2009/01/27 10:44:58 gcosmo Exp $
-// GEANT4 tag $Name: geant4-09-03-cand-01 $
+// $Id: DicomNestedPhantomParameterisation.cc,v 1.6 2009/11/10 18:48:46 arce Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-01 $
 //
 // --------------------------------------------------------------------
@@ -53 +53 @@
 
 void DicomNestedPhantomParameterisation::
-SetNoVoxel( size_t nx, size_t ny, size_t nz )
+SetNoVoxel( unsigned int nx, unsigned int ny, unsigned int nz )
 {
   fnX = nx;
@@ -80 +80 @@
   G4int copyNo = ix + fnX*iy + fnX*fnY*iz;
 
-  size_t matIndex = GetMaterialIndex(copyNo);
+  unsigned int matIndex = GetMaterialIndex(copyNo);
 
   return fMaterials[ matIndex ];
@@ -86 +86 @@
 
 //------------------------------------------------------------------
-size_t DicomNestedPhantomParameterisation::
-GetMaterialIndex( size_t copyNo ) const
+unsigned int DicomNestedPhantomParameterisation::
+GetMaterialIndex( unsigned int copyNo ) const
 {
   return *(fMaterialIndices+copyNo);

trunk/examples/extended/medical/DICOM/src/DicomPatientZSliceHeader.cc

@@ -148 +148 @@
     G4Exception("");
   }
-  for( size_t ii = 0; ii < fMaterialNames.size(); ii++ ) {
+  for( unsigned int ii = 0; ii < fMaterialNames.size(); ii++ ) {
     if( fMaterialNames[ii] != fMaterialNames2[ii] ) {
       G4cerr << "DicomPatientZSliceHeader error adding two slice headers: !!! Different material number " << ii << " : " << fMaterialNames[ii] << " =? " << fMaterialNames2[ii] << G4endl;

trunk/examples/extended/medical/DICOM/src/DicomPhantomParameterisationColour.cc

@@ -77 +77 @@
   if( physVol ) {
     G4String mateName = mate->GetName();
-    size_t iuu = mateName.find("__");
+    unsigned int iuu = mateName.find("__");
     if( iuu != std::string::npos ) {
       mateName = mateName.substr( 0, iuu );

trunk/examples/extended/medical/DICOM/src/DicomPhysicsList.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// The code was written by :
-//      *Louis Archambault louis.archambault@phy.ulaval.ca,
-//      *Luc Beaulieu beaulieu@phy.ulaval.ca
-//      +Vincent Hubert-Tremblay at tigre.2@sympatico.ca
-//
-//
-// *Centre Hospitalier Universitaire de Quebec (CHUQ),
-// Hotel-Dieu de Quebec, departement de Radio-oncologie
-// 11 cote du palais. Quebec, QC, Canada, G1R 2J6
-// tel (418) 525-4444 #6720
-// fax (418) 691 5268
-//
-// + Université Laval, Québec (QC) Canada
-//
-// History: 30.11.07  P.Arce default cut changed to 1 mm
-//*******************************************************
-
-#include "DicomPhysicsList.hh"
+// -------------------------------------------------------------------
+// $Id: DicomPhysicsList.cc,v 1.11 2009/10/26 11:20:49 chauvie Exp $
+// -------------------------------------------------------------------
 
 #include "G4ParticleDefinition.hh"
-#include "G4ParticleWithCuts.hh"
 #include "G4ProcessManager.hh"
 #include "G4ParticleTypes.hh"
-#include "G4ParticleTable.hh"
-#include "G4Material.hh"
-#include "G4UnitsTable.hh"
-#include "G4ios.hh"
+#include "G4StepLimiter.hh"
+#include "G4BaryonConstructor.hh"                     
+#include "G4IonConstructor.hh"   
+#include "G4MesonConstructor.hh"         
+
+#include "DicomPhysicsList.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
 DicomPhysicsList::DicomPhysicsList():  G4VUserPhysicsList()
 {
-  defaultCutValue = 1.*mm;
-  cutForGamma     = 1.*mm;
+  defaultCutValue = 0.01*micrometer;
+  cutForGamma     = defaultCutValue;
   cutForElectron  = defaultCutValue;
   cutForPositron  = defaultCutValue;
-
-  SetVerboseLevel(0);
-}
+  SetVerboseLevel(1);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
 DicomPhysicsList::~DicomPhysicsList()
 {}
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void DicomPhysicsList::ConstructParticle()
 {
-  // In this method, static member functions should be called
-  // for all particles which you want to use.
-  // This ensures that objects of these particle types will be
-  // created in the program.
-
   ConstructBosons();
   ConstructLeptons();
-}
+  ConstructBaryons();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
 void DicomPhysicsList::ConstructBosons()
-{
+{ 
   // gamma
   G4Gamma::GammaDefinition();
 
-}
+  // optical photon
+  G4OpticalPhoton::OpticalPhotonDefinition();
+}
+ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
 void DicomPhysicsList::ConstructLeptons()
@@ -90 +82 @@
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void DicomPhysicsList::ConstructBaryons()
+{
+  //  baryons
+  G4BaryonConstructor bConstructor;
+  bConstructor.ConstructParticle();
+
+  G4IonConstructor iConstructor;
+  iConstructor.ConstructParticle();
+
+  G4MesonConstructor mConstructor;
+  mConstructor.ConstructParticle();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void DicomPhysicsList::ConstructProcess()
 {
   AddTransportation();
   ConstructEM();
-}
-
-#include "G4MultipleScattering.hh"
+  ConstructHad();
+  ConstructGeneral();
+}
+
+// *** Processes and models
+
 // gamma
-#include "G4LowEnergyRayleigh.hh"
-#include "G4LowEnergyPhotoElectric.hh"
-#include "G4LowEnergyCompton.hh"
-#include "G4LowEnergyGammaConversion.hh"
+
+#include "G4PhotoElectricEffect.hh"
+#include "G4LivermorePhotoElectricModel.hh"
+
+#include "G4ComptonScattering.hh"
+#include "G4LivermoreComptonModel.hh"
+
+#include "G4GammaConversion.hh"
+#include "G4LivermoreGammaConversionModel.hh"
+
+#include "G4RayleighScattering.hh" 
+#include "G4LivermoreRayleighModel.hh"
+
 // e-
-#include "G4LowEnergyIonisation.hh"
-#include "G4LowEnergyBremsstrahlung.hh"
+
+#include "G4eMultipleScattering.hh"
+#include "G4UniversalFluctuation.hh"
+
+#include "G4eIonisation.hh"
+#include "G4LivermoreIonisationModel.hh"
+
+#include "G4eBremsstrahlung.hh"
+#include "G4LivermoreBremsstrahlungModel.hh"
+
 // e+
-#include "G4eIonisation.hh"
-#include "G4eBremsstrahlung.hh"
+
 #include "G4eplusAnnihilation.hh"
 
+// mu
+
+#include "G4MuIonisation.hh"
+#include "G4MuBremsstrahlung.hh"
+#include "G4MuPairProduction.hh"
+
+// hadrons
+
+#include "G4hMultipleScattering.hh"
+#include "G4MscStepLimitType.hh"
+
+#include "G4hBremsstrahlung.hh"
+#include "G4hPairProduction.hh"
+
+#include "G4hIonisation.hh"
+#include "G4ionIonisation.hh"
+#include "G4IonParametrisedLossModel.hh"
+
+//
+
+#include "G4LossTableManager.hh"
+#include "G4EmProcessOptions.hh"
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 void DicomPhysicsList::ConstructEM()
 {
+  theParticleIterator->reset();
+
+  while( (*theParticleIterator)() ){
+
+    G4ParticleDefinition* particle = theParticleIterator->value();
+    G4ProcessManager* pmanager = particle->GetProcessManager();
+    G4String particleName = particle->GetParticleName();
+
+    if (particleName == "gamma") {
+
+
+      G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
+      G4LivermorePhotoElectricModel* theLivermorePhotoElectricModel = new G4LivermorePhotoElectricModel();
+      thePhotoElectricEffect->AddEmModel(0, theLivermorePhotoElectricModel);
+      pmanager->AddDiscreteProcess(thePhotoElectricEffect);
+
+      G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
+      G4LivermoreComptonModel* theLivermoreComptonModel = new G4LivermoreComptonModel();
+      theComptonScattering->AddEmModel(0, theLivermoreComptonModel);
+      pmanager->AddDiscreteProcess(theComptonScattering);
+
+      G4GammaConversion* theGammaConversion = new G4GammaConversion();
+      G4LivermoreGammaConversionModel* theLivermoreGammaConversionModel = new G4LivermoreGammaConversionModel();
+      theGammaConversion->AddEmModel(0, theLivermoreGammaConversionModel);
+      pmanager->AddDiscreteProcess(theGammaConversion);
+
+      G4RayleighScattering* theRayleigh = new G4RayleighScattering();
+      G4LivermoreRayleighModel* theRayleighModel = new G4LivermoreRayleighModel();
+      theRayleigh->AddEmModel(0, theRayleighModel);
+      pmanager->AddDiscreteProcess(theRayleigh);
+      
+      pmanager->AddProcess(new G4StepLimiter(), -1, -1, 5);
+      
+    } else if (particleName == "e-") {
+     
+      G4eMultipleScattering* msc = new G4eMultipleScattering();
+      msc->SetStepLimitType(fUseDistanceToBoundary);
+      pmanager->AddProcess(msc,                   -1, 1, 1);
+      
+      // Ionisation
+      G4eIonisation* eIoni = new G4eIonisation();
+      eIoni->AddEmModel(0, new G4LivermoreIonisationModel(), new G4UniversalFluctuation() );
+      eIoni->SetStepFunction(0.2, 100*um); //     
+      pmanager->AddProcess(eIoni,                 -1, 2, 2);
+      
+      // Bremsstrahlung
+      G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
+      eBrem->AddEmModel(0, new G4LivermoreBremsstrahlungModel());
+      pmanager->AddProcess(eBrem,                 -1,-3, 3);
+
+      pmanager->AddProcess(new G4StepLimiter(), -1, -1, 4);
+      
+    } else if (particleName == "e+") {
+
+      // Identical to G4EmStandardPhysics_option3
+    
+      G4eMultipleScattering* msc = new G4eMultipleScattering();
+      msc->SetStepLimitType(fUseDistanceToBoundary);
+      pmanager->AddProcess(msc,                   -1, 1, 1);
+
+      G4eIonisation* eIoni = new G4eIonisation();
+      eIoni->SetStepFunction(0.2, 100*um);      
+      pmanager->AddProcess(eIoni,                 -1, 2, 2);
+      
+      pmanager->AddProcess(new G4eBremsstrahlung, -1,-3, 3);
+      
+      pmanager->AddProcess(new G4eplusAnnihilation,0,-1, 4);
+      
+      pmanager->AddProcess(new G4StepLimiter(), -1, -1, 5);
+
+    } else if (particleName == "GenericIon") {
+
+      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+
+      G4ionIonisation* ionIoni = new G4ionIonisation();
+      ionIoni->SetEmModel(new G4IonParametrisedLossModel());
+      ionIoni->SetStepFunction(0.1, 20*um);
+      pmanager->AddProcess(ionIoni,                   -1, 2, 2);
+
+      pmanager->AddProcess(new G4StepLimiter(), -1, -1, 3);
+
+    } else if (particleName == "alpha" ||
+               particleName == "He3" ) {
+
+      // Identical to G4EmStandardPhysics_option3
+      
+      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+
+      G4ionIonisation* ionIoni = new G4ionIonisation();
+      ionIoni->SetStepFunction(0.1, 20*um);
+      pmanager->AddProcess(ionIoni,                   -1, 2, 2);
+
+      pmanager->AddProcess(new G4StepLimiter(), -1, -1, 3);
+    }
+
+   //
+
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4HadronElasticProcess.hh"
+#include "G4LElastic.hh"
+
+#include "G4AlphaInelasticProcess.hh"
+#include "G4BinaryLightIonReaction.hh"
+#include "G4TripathiCrossSection.hh"
+#include "G4IonsShenCrossSection.hh"
+#include "G4LEAlphaInelastic.hh"
+
+void DicomPhysicsList::ConstructHad()
+{
+
+  G4HadronElasticProcess * theElasticProcess = new G4HadronElasticProcess;
+  theElasticProcess->RegisterMe( new G4LElastic() );
+
   theParticleIterator->reset();
   while( (*theParticleIterator)() )
-    {
-      G4ParticleDefinition* particle = theParticleIterator->value();
-      G4ProcessManager* pmanager = particle->GetProcessManager();
-      G4String particleName = particle->GetParticleName();
-
-      //processes
-      lowePhot = new  G4LowEnergyPhotoElectric("LowEnPhotoElec");
-      loweIon  = new G4LowEnergyIonisation("LowEnergyIoni");
-      loweBrem = new G4LowEnergyBremsstrahlung("LowEnBrem");
-
-      if (particleName == "gamma")
-        {
-          //gamma
-          pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh);
-          pmanager->AddDiscreteProcess(lowePhot);
-          pmanager->AddDiscreteProcess(new G4LowEnergyCompton);
-          pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion);
-
-        }
-      else if (particleName == "e-")
-        {
-          //electron
-          pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
-          pmanager->AddProcess(loweIon,     -1, 2,2);
-          pmanager->AddProcess(loweBrem,    -1,-1,3);
-
-        }
-      else if (particleName == "e+")
-        {
-          //positron
-          pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
-          pmanager->AddProcess(new G4eIonisation,        -1, 2,2);
-          pmanager->AddProcess(new G4eBremsstrahlung,    -1,-1,3);
-          pmanager->AddProcess(new G4eplusAnnihilation,   0,-1,4);
-
-        }
-    }
-}
-#include "G4Decay.hh"
+  {
+    G4ParticleDefinition* particle = theParticleIterator->value();
+    G4ProcessManager* pManager = particle->GetProcessManager();
+
+    if (particle->GetParticleName() == "alpha") 
+       { 
+
+          // INELASTIC SCATTERING
+          // Binary Cascade
+          G4BinaryLightIonReaction* theBC = new G4BinaryLightIonReaction();
+          theBC -> SetMinEnergy(80.*MeV);
+          theBC -> SetMaxEnergy(40.*GeV);
+  
+          // TRIPATHI CROSS SECTION
+          // Implementation of formulas in analogy to NASA technical paper 3621 by 
+          // Tripathi, et al. Cross-sections for ion ion scattering
+          G4TripathiCrossSection* TripathiCrossSection = new G4TripathiCrossSection;
+  
+          // IONS SHEN CROSS SECTION
+          // Implementation of formulas 
+          // Shen et al. Nuc. Phys. A 491 130 (1989) 
+          // Total Reaction Cross Section for Heavy-Ion Collisions
+          G4IonsShenCrossSection* aShen = new G4IonsShenCrossSection;
+  
+          // Final state production model for Alpha inelastic scattering below 20 GeV
+          G4LEAlphaInelastic* theAIModel = new G4LEAlphaInelastic;
+          theAIModel -> SetMaxEnergy(100.*MeV);
+
+          G4AlphaInelasticProcess * theIPalpha = new G4AlphaInelasticProcess;             
+          theIPalpha->AddDataSet(TripathiCrossSection);
+          theIPalpha->AddDataSet(aShen);
+
+          // Register the Alpha Inelastic and Binary Cascade Model
+          theIPalpha->RegisterMe(theAIModel);
+          theIPalpha->RegisterMe(theBC);
+          
+          // Activate the alpha inelastic scattering using the alpha inelastic and binary cascade model
+          pManager -> AddDiscreteProcess(theIPalpha);
+          
+          // Activate the Hadron Elastic Process
+          pManager -> AddDiscreteProcess(theElasticProcess); 
+            
+       }
+  }
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 void DicomPhysicsList::ConstructGeneral()
-{
-  // Add Decay Process
-  G4Decay* theDecayProcess = new G4Decay();
-  theParticleIterator->reset();
-  while( (*theParticleIterator)() )
-    {
-      G4ParticleDefinition* particle = theParticleIterator->value();
-      G4ProcessManager* pmanager = particle->GetProcessManager();
-      if (theDecayProcess->IsApplicable(*particle))
-        {
-          pmanager ->AddProcess(theDecayProcess);
-          // set ordering for PostStepDoIt and AtRestDoIt
-          pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
-          pmanager ->SetProcessOrdering(theDecayProcess, idxAtRest);
-        }
-    }
-}
+{ }
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
 void DicomPhysicsList::SetCuts()
 {
-  if (verboseLevel >0)
-    {
-      G4cout << "DicomPhysicsList::SetCuts:";
-      G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
-    }
-
+  if (verboseLevel >0){
+    G4cout << "DicomPhysicsList::SetCuts:";
+    G4cout << "CutLength : " << G4BestUnit(defaultCutValue,"Length") << G4endl;
+  }  
+  
   // set cut values for gamma at first and for e- second and next for e+,
-  // because some processes for e+/e- need cut values for gamma
+  // because some processes for e+/e- need cut values for gamma 
   SetCutValue(cutForGamma, "gamma");
   SetCutValue(cutForElectron, "e-");
   SetCutValue(cutForPositron, "e+");
-
-
-  if (verboseLevel>0)
-    DumpCutValuesTable();
-}
-
-void DicomPhysicsList::SetGammaCut(G4double val)
-{
-  ResetCuts();
-  cutForGamma = val;
-}
-
-void DicomPhysicsList::SetElectronCut(G4double val)
-{
-  //  ResetCuts();
-  cutForElectron = val;
-}
-
-void DicomPhysicsList::SetPositronCut(G4double val)
-{
-  //  ResetCuts();
-  cutForPositron = val;
-}
-
+  
+  if (verboseLevel>0) DumpCutValuesTable();
+
+}
+

trunk/examples/extended/medical/DICOM/src/DicomPrimaryGeneratorAction.cc

@@ -74 +74 @@
   particleGun->SetParticleEnergy(5.*MeV);
   //put it at SAD = 1 m on xy plane of central slice
-  particleGun->SetParticlePosition(G4ThreeVector(0.,-99.9*cm,-27.));
+  particleGun->SetParticlePosition(G4ThreeVector(0.,-99.9*cm,87.5*mm));
   //particleGun->SetParticlePosition(G4ThreeVector(0.,0.,-22.));
   particleGun->GeneratePrimaryVertex(anEvent);