Index: trunk/examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc
===================================================================
--- trunk/examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc	(revision 1230)
+++ trunk/examples/extended/medical/DICOM/src/DicomDetectorConstruction.cc	(revision 1337)
@@ -314,5 +314,5 @@
   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,8 +327,10 @@
   //--- 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,5 +355,5 @@
     
     //-- 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,5 +379,5 @@
   //----- 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];
   };
Index: trunk/examples/extended/medical/DICOM/src/DicomEventAction.cc
===================================================================
--- trunk/examples/extended/medical/DICOM/src/DicomEventAction.cc	(revision 1230)
+++ trunk/examples/extended/medical/DICOM/src/DicomEventAction.cc	(revision 1337)
@@ -42,7 +42,4 @@
 #include "DicomEventAction.hh"
 #include "G4Event.hh"
-#include "G4TrajectoryContainer.hh"
-#include "G4Trajectory.hh"
-#include "G4VVisManager.hh"
 
 DicomEventAction::DicomEventAction():drawFlag("all")
@@ -55,28 +52,7 @@
 {
   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*)
+{
+}
Index: trunk/examples/extended/medical/DICOM/src/DicomHandler.cc
===================================================================
--- trunk/examples/extended/medical/DICOM/src/DicomHandler.cc	(revision 1230)
+++ trunk/examples/extended/medical/DICOM/src/DicomHandler.cc	(revision 1337)
@@ -61,5 +61,5 @@
 
 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,8 +98,10 @@
     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,15 +119,19 @@
 	// 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,34 +142,58 @@
 	    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,10 +209,10 @@
     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,4 +224,5 @@
     foutG4DCM << -pixelSpacingY*columns/2 << " " << pixelSpacingY*columns/2 << G4endl;
     foutG4DCM << sliceLocation-sliceThickness/2. << " " << sliceLocation+sliceThickness/2. << G4endl;
+    //    foutG4DCM << compression << G4endl;
     
     ReadData( dicom, filename2 );
@@ -431,12 +462,12 @@
 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,8 +477,7 @@
 }
 
-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,19 +538,46 @@
 
     // 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,44 +589,84 @@
     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,7 +692,7 @@
 */
 
-G4double DicomHandler::Pixel2density(G4int pixel)
+G4float DicomHandler::Pixel2density(G4int pixel)
 {
-    G4double density = -1.;
+    G4float density = -1.;
     G4int nbrequali = 0;
     G4double deltaCT = 0;
@@ -752,2 +849,99 @@
     _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;
+ 
+}
Index: trunk/examples/extended/medical/DICOM/src/DicomNestedPhantomParameterisation.cc
===================================================================
--- trunk/examples/extended/medical/DICOM/src/DicomNestedPhantomParameterisation.cc	(revision 1230)
+++ trunk/examples/extended/medical/DICOM/src/DicomNestedPhantomParameterisation.cc	(revision 1337)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $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,5 +53,5 @@
 
 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,5 +80,5 @@
   G4int copyNo = ix + fnX*iy + fnX*fnY*iz;
 
-  size_t matIndex = GetMaterialIndex(copyNo);
+  unsigned int matIndex = GetMaterialIndex(copyNo);
 
   return fMaterials[ matIndex ];
@@ -86,6 +86,6 @@
 
 //------------------------------------------------------------------
-size_t DicomNestedPhantomParameterisation::
-GetMaterialIndex( size_t copyNo ) const
+unsigned int DicomNestedPhantomParameterisation::
+GetMaterialIndex( unsigned int copyNo ) const
 {
   return *(fMaterialIndices+copyNo);
Index: trunk/examples/extended/medical/DICOM/src/DicomPatientZSliceHeader.cc
===================================================================
--- trunk/examples/extended/medical/DICOM/src/DicomPatientZSliceHeader.cc	(revision 1230)
+++ trunk/examples/extended/medical/DICOM/src/DicomPatientZSliceHeader.cc	(revision 1337)
@@ -148,5 +148,5 @@
     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;
Index: trunk/examples/extended/medical/DICOM/src/DicomPhantomParameterisationColour.cc
===================================================================
--- trunk/examples/extended/medical/DICOM/src/DicomPhantomParameterisationColour.cc	(revision 1230)
+++ trunk/examples/extended/medical/DICOM/src/DicomPhantomParameterisationColour.cc	(revision 1337)
@@ -77,5 +77,5 @@
   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 );
Index: trunk/examples/extended/medical/DICOM/src/DicomPhysicsList.cc
===================================================================
--- trunk/examples/extended/medical/DICOM/src/DicomPhysicsList.cc	(revision 1230)
+++ trunk/examples/extended/medical/DICOM/src/DicomPhysicsList.cc	(revision 1337)
@@ -24,62 +24,54 @@
 // ********************************************************************
 //
-// 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,121 +82,276 @@
 }
 
+//....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();
+
+}
+
Index: trunk/examples/extended/medical/DICOM/src/DicomPrimaryGeneratorAction.cc
===================================================================
--- trunk/examples/extended/medical/DICOM/src/DicomPrimaryGeneratorAction.cc	(revision 1230)
+++ trunk/examples/extended/medical/DICOM/src/DicomPrimaryGeneratorAction.cc	(revision 1337)
@@ -74,5 +74,5 @@
   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);
