Index: trunk/examples/extended/electromagnetic/TestEm9/src/DetectorConstruction.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/DetectorConstruction.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/DetectorConstruction.cc	(revision 1230)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: DetectorConstruction.cc,v 1.9 2006/06/29 17:02:59 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: DetectorConstruction.cc,v 1.11 2008/04/07 18:09:05 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //
@@ -58,4 +58,5 @@
 #include "G4Region.hh"
 #include "G4RegionStore.hh"
+#include "G4ProductionCuts.hh"
 #include "G4PhysicalVolumeStore.hh"
 #include "G4LogicalVolumeStore.hh"
@@ -72,5 +73,5 @@
 
 DetectorConstruction::DetectorConstruction()
-  :logicC(0),logicA1(0),logicA2(0)
+  : G4VUserDetectorConstruction()
 {
   detectorMessenger = new DetectorMessenger(this);
@@ -84,5 +85,10 @@
   vertexRegion = 0;
   muonRegion   = 0;
+  logicWorld   = 0;
+  logicCal     = 0;
+  logicA1      = 0;
   DefineMaterials();
+  vertexDetectorCuts = new G4ProductionCuts();
+  muonDetectorCuts   = new G4ProductionCuts();
 }
 
@@ -90,5 +96,9 @@
 
 DetectorConstruction::~DetectorConstruction()
-{ delete detectorMessenger;}
+{ 
+  delete detectorMessenger;
+  delete vertexDetectorCuts;
+  delete muonDetectorCuts;
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -106,5 +116,5 @@
 
   G4NistManager* man = G4NistManager::Instance();
-  man->SetVerbose(1);
+  //  man->SetVerbose(1);
   worldMaterial = man->FindOrBuildMaterial("G4_AIR");
   absMaterial   = man->FindOrBuildMaterial("G4_Al");
@@ -118,17 +128,24 @@
 G4VPhysicalVolume* DetectorConstruction::ConstructVolumes()
 {
-  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
-
   // Cleanup old geometry
 
   G4GeometryManager::GetInstance()->OpenGeometry();
+
+  if(G4NistManager::Instance()->GetVerbose() > 0)
+    G4cout << *(G4Material::GetMaterialTable()) << G4endl;
+
+  if(vertexRegion) {
+    delete vertexRegion;
+    delete muonRegion;
+  }
+  vertexRegion = new G4Region("VertexDetector");
+  vertexRegion->SetProductionCuts(vertexDetectorCuts);
+
+  muonRegion   = new G4Region("MuonDetector");
+  muonRegion->SetProductionCuts(muonDetectorCuts);
+
+  G4SolidStore::GetInstance()->Clean();
+  G4LogicalVolumeStore::GetInstance()->Clean();
   G4PhysicalVolumeStore::GetInstance()->Clean();
-  G4LogicalVolumeStore::GetInstance()->Clean();
-  G4SolidStore::GetInstance()->Clean();
-
-  if(vertexRegion) delete vertexRegion;
-  if(muonRegion) delete muonRegion;
-  vertexRegion = new G4Region("VertexDetector");
-  muonRegion   = new G4Region("MuonDetector");
 
   if(vertexLength < padLength*5.0) vertexLength = padLength*5.0;
@@ -137,11 +154,13 @@
   G4double york   = 10.*cm;
 
-           worldZ = 2.*vertexLength + 3.*absLength + 0.5*(ecalLength + york) + biggap*2.;
+  worldZ = 2.*vertexLength + 3.*absLength + 0.5*(ecalLength + york) + biggap*2.;
+
   G4double worldX = ecalWidth*3.0;
   G4double vertexZ= -worldZ + vertexLength*2.0 + absLength     + biggap;
   G4double absZ2  = -worldZ + vertexLength*4.0 + absLength*3.5 + biggap;
-  G4double ecalZ  = -worldZ + vertexLength*4.0 + absLength*4.0 + ecalLength*0.5 + 2.*biggap;
+  G4double ecalZ  = -worldZ + vertexLength*4.0 + absLength*4.0 + ecalLength*0.5 
+    + 2.*biggap;
   G4double yorkZ  = -worldZ + vertexLength*4.0 + absLength*5.0 + ecalLength
-                            + york*0.5 + 3.*biggap;
+    + york*0.5 + 3.*biggap;
 
   //
@@ -149,8 +168,7 @@
   //
   G4Box* solidW = new G4Box("World",worldX,worldX,worldZ);
-  G4LogicalVolume* logicW = new G4LogicalVolume( solidW,worldMaterial,
-                                                "World");
+  logicWorld = new G4LogicalVolume( solidW,worldMaterial,"World");
   G4VPhysicalVolume* world = new G4PVPlacement(0,G4ThreeVector(),
-                                       "World",logicW,0,false,0);
+					       "World",logicWorld,0,false,0);
 
   //
@@ -158,14 +176,13 @@
   //
   G4Box* solidE = new G4Box("VolE",worldX,worldX,ecalLength*0.5 + gap);
-  G4LogicalVolume* logicE = new G4LogicalVolume( solidE,worldMaterial,
-                                                "VolE");
+  logicECal = new G4LogicalVolume( solidE,worldMaterial,"VolE");
   G4VPhysicalVolume* physE = new G4PVPlacement(0,G4ThreeVector(0.,0.,ecalZ),
-                                       "VolE",logicE,world,false,0);
+					       "VolE",logicECal,world,false,0);
 
   G4Box* solidC = new G4Box("Ecal",ecalWidth*0.5,ecalWidth*0.5,ecalLength*0.5);
-  logicC = new G4LogicalVolume( solidC,calMaterial,"Ecal");
+  logicCal = new G4LogicalVolume( solidC,calMaterial,"Ecal");
 
   G4cout << "Ecal is " << G4BestUnit(ecalLength,"Length")
-       << " of " << calMaterial->GetName() << G4endl;
+	 << " of " << calMaterial->GetName() << G4endl;
 
   // Crystals
@@ -182,5 +199,5 @@
     for (j=0; j<5; j++) {
 
-      pv = new G4PVPlacement(0,G4ThreeVector(x,y,0.),"Ecal",logicC,
+      pv = new G4PVPlacement(0,G4ThreeVector(x,y,0.),"Ecal",logicCal,
                                     physE,false,k);
       k++;
@@ -195,43 +212,41 @@
   logicA2 = new G4LogicalVolume( solidA,absMaterial,"Abs2");
   pv = new G4PVPlacement(0,G4ThreeVector(0.,0.,absZ2),
-                                       "Abs2",logicA2,world,false,0);
+			 "Abs2",logicA2,world,false,0);
 
   G4cout << "Absorber is " << G4BestUnit(absLength,"Length")
-       << " of " << absMaterial->GetName() << G4endl;
+	 << " of " << absMaterial->GetName() << G4endl;
 
   //York
 
   G4Box* solidYV = new G4Box("VolY",worldX,worldX,york*0.5+absLength);
-  G4LogicalVolume* logicYV = new G4LogicalVolume( solidYV,yorkMaterial,"VolY");
+  logicYV = new G4LogicalVolume( solidYV,yorkMaterial,"VolY");
   G4VPhysicalVolume* physYV = new G4PVPlacement(0,G4ThreeVector(0.,0.,yorkZ),
-                                       "VolY",logicYV,world,false,0);
+						"VolY",logicYV,world,false,0);
 
   G4Box* solidY = new G4Box("York",worldX,worldX,york*0.5);
-  G4LogicalVolume* logicY = new G4LogicalVolume( solidY,yorkMaterial,"York");
+  logicY = new G4LogicalVolume( solidY,yorkMaterial,"York");
   pv = new G4PVPlacement(0,G4ThreeVector(),
-                                       "York",logicY,physYV,false,0);
+			 "York",logicY,physYV,false,0);
 
   logicA3 = new G4LogicalVolume( solidA,absMaterial,"Abs3");
   logicA4 = new G4LogicalVolume( solidA,absMaterial,"Abs4");
+
   pv = new G4PVPlacement(0,G4ThreeVector(0.,0.,-(york+absLength)*0.5),
-                                       "Abs3",logicA3,physYV,false,0);
+			 "Abs3",logicA3,physYV,false,0);
   pv = new G4PVPlacement(0,G4ThreeVector(0.,0.,(york+absLength)*0.5),
-                                       "Abs4",logicA4,physYV,false,0);
+			 "Abs4",logicA4,physYV,false,0);
 
   //Vertex volume
-
   G4Box* solidVV = new G4Box("VolV",worldX,worldX,vertexLength*2.+absLength+gap);
-  G4LogicalVolume* logicVV = new G4LogicalVolume( solidVV,worldMaterial,"VolV");
+  logicVV = new G4LogicalVolume( solidVV,worldMaterial,"VolV");
   G4VPhysicalVolume* physVV = new G4PVPlacement(0,G4ThreeVector(0.,0.,vertexZ),
-                                       "VolV",logicVV,world,false,0);
+						"VolV",logicVV,world,false,0);
 
   //Absorber
-
   logicA1 = new G4LogicalVolume( solidA,absMaterial,"Abs1");
   pv = new G4PVPlacement(0,G4ThreeVector(0.,0.,vertexLength*2.-absLength*0.5),
-                                       "Abs1",logicA1,physVV,false,0);
+			 "Abs1",logicA1,physVV,false,0);
 
   //Vertex
-
   G4double vertWidth = ecalWidth/5.;
   G4int npads = (G4int)(vertWidth/padWidth);
@@ -242,8 +257,9 @@
 
   G4Box* solidVD = new G4Box("VertDet",x1,ecalWidth*0.5+gap,padLength*0.5);
-  G4LogicalVolume* logicVD = new G4LogicalVolume( solidVD,vertMaterial,"VertDet");
+  logicVD = new G4LogicalVolume( solidVD,vertMaterial,"VertDet");
+  logicVD->SetSolid(solidVD);
 
   G4Box* solidV = new G4Box("Vert",padWidth*0.5,ecalWidth*0.5,padLength*0.5);
-  G4LogicalVolume* logicV = new G4LogicalVolume( solidV,vertMaterial,"Vert");
+  logicV = new G4LogicalVolume( solidV,vertMaterial,"Vert");
 
   for (i=0; i<3; i++) {
@@ -267,21 +283,18 @@
 
   // Define region for the vertex detector
-
   vertexRegion->AddRootLogicalVolume(logicVV);
   vertexRegion->AddRootLogicalVolume(logicA3);
 
   // Define region for the muon detector
-
   muonRegion->AddRootLogicalVolume(logicYV);
 
   // color regions
-
   logicVV-> SetVisAttributes(G4VisAttributes::Invisible);
   logicV-> SetVisAttributes(G4VisAttributes::Invisible);
-  logicE-> SetVisAttributes(G4VisAttributes::Invisible);
+  logicECal-> SetVisAttributes(G4VisAttributes::Invisible);
   logicYV-> SetVisAttributes(G4VisAttributes::Invisible);
 
   G4VisAttributes* regWcolor = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3));
-  logicW->SetVisAttributes(regWcolor);
+  logicWorld->SetVisAttributes(regWcolor);
 
   G4VisAttributes* regVcolor = new G4VisAttributes(G4Colour(0., 0.3, 0.7));
@@ -289,5 +302,5 @@
 
   G4VisAttributes* regCcolor = new G4VisAttributes(G4Colour(0., 0.7, 0.3));
-  logicC->SetVisAttributes(regCcolor);
+  logicCal->SetVisAttributes(regCcolor);
 
   G4VisAttributes* regAcolor = new G4VisAttributes(G4Colour(1., 0.5, 0.5));
@@ -302,5 +315,5 @@
   // always return world
   G4cout << "### New geometry is constructed" << G4endl;
-
+ 
   return world;
 }
@@ -313,5 +326,11 @@
   G4Material* pttoMaterial = 
     G4NistManager::Instance()->FindOrBuildMaterial(mat, false);
-  if (pttoMaterial) calMaterial = pttoMaterial;
+  if (pttoMaterial) {
+    calMaterial = pttoMaterial;
+    if(logicCal) {
+      logicCal->SetMaterial(calMaterial);
+      G4RunManager::GetRunManager()->PhysicsHasBeenModified();
+    }
+  }
 }
 
@@ -323,5 +342,14 @@
   G4Material* pttoMaterial = 
     G4NistManager::Instance()->FindOrBuildMaterial(mat, false);
-  if (pttoMaterial) absMaterial = pttoMaterial;
+  if (pttoMaterial) {
+    absMaterial = pttoMaterial;
+    if(logicA1) {
+      logicA1->SetMaterial(absMaterial);
+      logicA2->SetMaterial(absMaterial);
+      logicA3->SetMaterial(absMaterial);
+      logicA4->SetMaterial(absMaterial);
+      G4RunManager::GetRunManager()->PhysicsHasBeenModified();
+    }
+  }
 }
 
@@ -330,6 +358,67 @@
 void DetectorConstruction::UpdateGeometry()
 {
+  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
   G4RunManager::GetRunManager()->DefineWorldVolume(ConstructVolumes());
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void DetectorConstruction::SetEcalLength (G4double val)   
+{
+  if(val > 0.0) {
+    ecalLength = val;
+    G4RunManager::GetRunManager()->GeometryHasBeenModified();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void DetectorConstruction::SetEcalWidth  (G4double val)   
+{
+  if(val > 0.0) {
+    ecalWidth = val;
+    G4RunManager::GetRunManager()->GeometryHasBeenModified();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void DetectorConstruction::SetVertexLength (G4double val) 
+{
+  if(val > 0.0) {
+    vertexLength = val;
+    G4RunManager::GetRunManager()->GeometryHasBeenModified();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void DetectorConstruction::SetPadLength  (G4double val)   
+{
+  if(val > 0.0) {
+    padLength = val;
+    G4RunManager::GetRunManager()->GeometryHasBeenModified();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void DetectorConstruction::SetPadWidth  (G4double val)    
+{
+  if(val > 0.0) {
+    padWidth = val;
+    G4RunManager::GetRunManager()->GeometryHasBeenModified();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void DetectorConstruction::SetAbsLength(G4double val)     
+{
+  if(val > 0.0) {
+    absLength = val;
+    G4RunManager::GetRunManager()->GeometryHasBeenModified();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm9/src/DetectorMessenger.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/DetectorMessenger.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/DetectorMessenger.cc	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: DetectorMessenger.cc,v 1.3 2006/06/29 17:03:01 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //
Index: trunk/examples/extended/electromagnetic/TestEm9/src/EmAcceptance.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/EmAcceptance.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/EmAcceptance.cc	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: EmAcceptance.cc,v 1.4 2006/06/29 17:03:03 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //
Index: trunk/examples/extended/electromagnetic/TestEm9/src/EventAction.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/EventAction.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/EventAction.cc	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: EventAction.cc,v 1.4 2006/06/29 17:03:06 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Index: trunk/examples/extended/electromagnetic/TestEm9/src/EventActionMessenger.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/EventActionMessenger.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/EventActionMessenger.cc	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: EventActionMessenger.cc,v 1.3 2006/06/29 17:03:08 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm9/src/Histo.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/Histo.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/Histo.cc	(revision 1230)
@@ -53,5 +53,5 @@
   verbose    = 0;
   histName   = "testem9";
-  histType   = "hbook";
+  histType   = "root";
   nHisto     = 0;
   defaultAct = 1;
@@ -93,6 +93,7 @@
 
   G4String nam = histName + "." + histType;
-
-  tree = tf->create(nam,histType,false,true,"uncompress");
+  G4String options  = "--noErrors export=root uncompress";
+  
+  tree = tf->create(nam,histType,false,true,options);
   delete tf;
   if(tree) {
@@ -107,6 +108,9 @@
   // Creating an 1-dimensional histograms in the root directory of the tree
   for(G4int i=0; i<nHisto; i++) {
-    if(active[i]) 
-      histo[i] = hf->createHistogram1D(ids[i], titles[i], bins[i], xmin[i], xmax[i]);
+    if(active[i]) {
+      G4String ss = ids[i];
+      if(histType == "root") ss = "h" + ids[i];
+      histo[i] = hf->createHistogram1D(ss, titles[i], bins[i], xmin[i], xmax[i]);
+    }
   }
   delete hf;
Index: trunk/examples/extended/electromagnetic/TestEm9/src/HistoManager.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/HistoManager.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/HistoManager.cc	(revision 1230)
@@ -87,5 +87,5 @@
   nBinsED= 100;
   nTuple = false;
-  nHisto = 10;
+  nHisto = 13;
 
   // initialise acceptance
@@ -97,11 +97,11 @@
 
   histo->add1D("10",
-    "Energy deposit (MeV) in central crystal",nBinsED,0.0,beamEnergy,MeV);
+    "Evis/E0 in central crystal",nBinsED,0.0,1,1.0);
 
   histo->add1D("11",
-    "Energy deposit (MeV) in 3x3",nBinsED,0.0,beamEnergy,MeV);
+    "Evis/E0 in 3x3",nBinsED,0.0,1.0,1.0);
 
   histo->add1D("12",
-    "Energy deposit (MeV) in 5x5",nBinsED,0.0,beamEnergy,MeV);
+    "Evis/E0 in 5x5",nBinsED,0.0,1.0,1.0);
 
   histo->add1D("13",
@@ -125,4 +125,13 @@
   histo->add1D("19",
     "Number of vertex hits",20,-0.5,19.5,1.0);
+
+  histo->add1D("20",
+    "E1/E9 Ratio",nBinsED,0.0,1,1.0);
+
+  histo->add1D("21",
+    "E1/25 Ratio",nBinsED,0.0,1.0,1.0);
+
+  histo->add1D("22",
+    "E9/E25 Ratio",nBinsED,0.0,1.0,1.0);
 
   if(nTuple) {
@@ -141,11 +150,14 @@
   n_gam  = 0;
   n_step = 0;
-
-  for(G4int i=0; i<nmax; i++) {
+  n_lowe = 0;
+
+  for(G4int i=0; i<6; i++) {
     stat[i] = 0;
     edep[i] = 0.0;
     erms[i] = 0.0;
-    edeptr[i] = 0.0;
-    ermstr[i] = 0.0;
+    if(i < 3) {
+      edeptr[i] = 0.0;
+      ermstr[i] = 0.0;
+    }
   }
 
@@ -164,5 +176,5 @@
 
   G4cout << "HistoManager: End of run actions are started" << G4endl;
-  G4String nam[3] = {"1x1", "3x3", "5x5"};
+  G4String nam[6] = {"1x1", "3x3", "5x5", "E1/E9 ", "E1/E25", "E9/E25"};
 
   // average
@@ -175,6 +187,6 @@
 
     // total mean
-    edep[j] *= x/beamEnergy;
-    G4double y = erms[j]*x/(beamEnergy*beamEnergy) - edep[j]*edep[j];
+    edep[j] *= x;
+    G4double y = erms[j]*x - edep[j]*edep[j];
     if(y < 0.0) y = 0.0;
     erms[j] = std::sqrt(y);
@@ -183,6 +195,6 @@
     G4double xx = G4double(stat[j]);
     if(xx > 0.0) xx = 1.0/xx;
-    edeptr[j] *= xx/beamEnergy;
-    y = ermstr[j]*xx/(beamEnergy*beamEnergy) - edeptr[j]*edeptr[j];
+    edeptr[j] *= xx;
+    y = ermstr[j]*xx - edeptr[j]*edeptr[j];
     if(y < 0.0) y = 0.0;
     ermstr[j] = std::sqrt(y);
@@ -191,5 +203,5 @@
   G4double xg = x*(G4double)n_gam;
   G4double xp = x*(G4double)n_posit;
-  G4double xs = x*(G4double)n_step;
+  G4double xs = x*n_step;
 
   G4double f = 100.*std::sqrt(beamEnergy/GeV);
@@ -201,8 +213,10 @@
   G4cout << std::setprecision(4) << "Average number of steps      " << xs << G4endl;
   
-  for(j=0; j<nmax; j++) {
+  for(j=0; j<3; j++) {
     G4double e = edeptr[j];
     G4double s = ermstr[j];
-    G4double r = s*std::sqrt(x);
+    G4double xx= G4double(stat[j]);
+    if(xx > 0.0) xx = 1.0/xx;
+    G4double r = s*std::sqrt(xx);
     G4cout << std::setprecision(4) << "Edep " << nam[j] << " =                   " << e
            << " +- " << r;
@@ -222,4 +236,20 @@
     }
   }
+  G4cout<<"===========  Ratios without trancating ==========================="<<G4endl;
+  for(j=3; j<6; j++) {
+    G4double e = edep[j];
+    G4double xx= G4double(stat[j]);
+    if(xx > 0.0) xx = 1.0/xx;
+    e *= xx;
+    G4double y = erms[j]*xx - e*e;
+    G4double r = 0.0;
+    if(y > 0.0) r = std::sqrt(y*xx);
+    G4cout << "  " << nam[j] << " =                   " << e
+           << " +- " << r;
+    G4cout << G4endl;
+  }
+  G4cout << std::setprecision(4) << "Beam Energy                  " << beamEnergy/GeV 
+	 << " GeV" << G4endl;
+  if(n_lowe > 0)          G4cout << "Number of events E/E0<0.8    " << n_lowe << G4endl; 
   G4cout<<"=================================================================="<<G4endl;
   G4cout<<G4endl;
@@ -264,5 +294,5 @@
   Evertex.clear();
   Nvertex.clear();
-  for (int i=0; i<25; i++) {
+  for (G4int i=0; i<25; i++) {
     E[i] = 0.0;
   }
@@ -275,9 +305,37 @@
   G4double e9 = 0.0;
   G4double e25= 0.0;
-  for (int i=0; i<25; i++) {
+  for (G4int i=0; i<25; i++) {
+    E[i] /= beamEnergy;
     e25 += E[i];
     if( ( 6<=i &&  8>=i) || (11<=i && 13>=i) || (16<=i && 18>=i)) e9 += E[i];
   }
-  histo->fill(0,E[12],1.0);
+
+  if(e25 < 0.8) {
+    n_lowe++;
+    G4cout << "### in the event# " << n_evt << "  E25= " << e25 << G4endl;
+  }
+
+  // compute ratios
+  G4double e0 = E[12];
+  G4double e19  = 0.0;
+  G4double e125 = 0.0;
+  G4double e925 = 0.0;
+  if(e9 > 0.0) {
+    e19 = e0/e9;
+    e125 = e0/e25;
+    e925 = e9/e25;
+    edep[3] += e19;
+    erms[3] += e19*e19;
+    edep[4] += e125;
+    erms[4] += e125*e125;
+    edep[5] += e925;
+    erms[5] += e925*e925;
+    stat[3] += 1;
+    stat[4] += 1;
+    stat[5] += 1;
+  }
+
+  // fill histo
+  histo->fill(0,e0,1.0);
   histo->fill(1,e9,1.0);
   histo->fill(2,e25,1.0);
@@ -286,9 +344,10 @@
   histo->fill(7,Eabs3,1.0);
   histo->fill(8,Eabs4,1.0);
-  float nn = (double)(Nvertex.size());
-  histo->fill(9,nn,1.0);
-
-  G4double e0 = E[12];
-
+  histo->fill(9,G4double(Nvertex.size()),1.0);
+  histo->fill(10,e19,1.0);
+  histo->fill(11,e125,1.0);
+  histo->fill(12,e925,1.0);
+
+  // compute sums
   edep[0] += e0;
   erms[0] += e0*e0;
@@ -299,15 +358,15 @@
 
   // trancated mean
-  if(limittrue[0] == DBL_MAX || std::abs(e0/beamEnergy-edeptrue[0])<rmstrue[0]*limittrue[0]) {
+  if(limittrue[0] == DBL_MAX || std::abs(e0-edeptrue[0])<rmstrue[0]*limittrue[0]) {
     stat[0] += 1;
     edeptr[0] += e0;
     ermstr[0] += e0*e0;
   }
-  if(limittrue[1] == DBL_MAX || std::abs(e9/beamEnergy-edeptrue[1])<rmstrue[1]*limittrue[1]) {
+  if(limittrue[1] == DBL_MAX || std::abs(e9-edeptrue[1])<rmstrue[1]*limittrue[1]) {
     stat[1] += 1;
     edeptr[1] += e9;
     ermstr[1] += e9*e9;
   }
-  if(limittrue[2] == DBL_MAX || std::abs(e25/beamEnergy-edeptrue[2])<rmstrue[2]*limittrue[2]) {
+  if(limittrue[2] == DBL_MAX || std::abs(e25-edeptrue[2])<rmstrue[2]*limittrue[2]) {
     stat[2] += 1;
     edeptr[2] += e25;
@@ -333,4 +392,5 @@
   if(0 == pid) {
 
+    beamEnergy = kinE;
     histo->fillTuple("TKIN", kinE/MeV);
 
Index: trunk/examples/extended/electromagnetic/TestEm9/src/HistoMessenger.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/HistoMessenger.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/HistoMessenger.cc	(revision 1230)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: HistoMessenger.cc,v 1.6 2006/06/29 17:03:16 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: HistoMessenger.cc,v 1.7 2008/08/22 14:11:53 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -60,5 +60,5 @@
   //
   G4UIparameter* ih = new G4UIparameter("ih",'i',false);
-  ih->SetGuidance("histo number : from 1 to MaxHisto");
+  ih->SetGuidance("histo number : from 0 to MaxHisto-1");
   histoCmd->SetParameter(ih);
   //
@@ -110,5 +110,6 @@
      G4String unit = unts;
      G4double vUnit = 1. ;
-     if (unit != "none") vUnit = G4UIcommand::ValueOf(unit);
+     if(unit != "none") vUnit = G4UIcommand::ValueOf(unit);
+     if(vUnit <= 0.0)  vUnit = 1.;
      histo->setHisto1D(ih,nbBins,vmin,vmax,vUnit);
    }
Index: trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmLivermore.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmLivermore.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmLivermore.cc	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: PhysListEmLivermore.cc,v 1.2 2007/01/08 16:29:42 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmPenelope.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmPenelope.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmPenelope.cc	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: PhysListEmPenelope.cc,v 1.1 2006/11/17 17:45:57 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmStandard.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmStandard.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmStandard.cc	(revision 1230)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: PhysListEmStandard.cc,v 1.8 2006/11/22 19:09:12 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: PhysListEmStandard.cc,v 1.13 2008/11/16 21:01:10 maire Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -39,5 +39,6 @@
 #include "G4PhotoElectricEffect.hh"
 
-#include "G4MultipleScattering.hh"
+#include "G4eMultipleScattering.hh"
+#include "G4hMultipleScattering.hh"
 
 #include "G4eIonisation.hh"
@@ -48,7 +49,12 @@
 #include "G4MuBremsstrahlung.hh"
 #include "G4MuPairProduction.hh"
+#include "G4hBremsstrahlung.hh"
+#include "G4hPairProduction.hh"
 
 #include "G4hIonisation.hh"
 #include "G4ionIonisation.hh"
+
+#include "G4EmProcessOptions.hh"
+#include "G4MscStepLimitType.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -83,22 +89,31 @@
     } else if (particleName == "e-") {
   
-      pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
-      pmanager->AddProcess(new G4eIonisation(),      -1, 2,2);
-      pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3,3);
+      pmanager->AddProcess(new G4eMultipleScattering, -1, 1, 1);
+      pmanager->AddProcess(new G4eIonisation(),       -1, 2, 2);
+      pmanager->AddProcess(new G4eBremsstrahlung,     -1, 3, 3);
 	    
     } else if (particleName == "e+") {
 
-      pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
-      pmanager->AddProcess(new G4eIonisation(),      -1, 2,2);
-      pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3,3);
-      pmanager->AddProcess(new G4eplusAnnihilation,   0,-1,4);
+      pmanager->AddProcess(new G4eMultipleScattering, -1, 1, 1);
+      pmanager->AddProcess(new G4eIonisation(),       -1, 2, 2);
+      pmanager->AddProcess(new G4eBremsstrahlung,     -1, 3, 3);
+      pmanager->AddProcess(new G4eplusAnnihilation,    0,-1, 4);
       
     } else if (particleName == "mu+" || 
                particleName == "mu-"    ) {
 
-      pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
-      pmanager->AddProcess(new G4MuIonisation,      -1, 2,2);
-      pmanager->AddProcess(new G4MuBremsstrahlung,  -1, 3,3);
-      pmanager->AddProcess(new G4MuPairProduction,  -1, 4,4);       
+      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+      pmanager->AddProcess(new G4MuIonisation,        -1, 2, 2);
+      pmanager->AddProcess(new G4MuBremsstrahlung,    -1, 3, 3);
+      pmanager->AddProcess(new G4MuPairProduction,    -1, 4, 4);       
+
+    } else if (particleName == "proton" || 
+               particleName == "pi-" ||  
+	       particleName == "pi+") {
+
+      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+      pmanager->AddProcess(new G4hIonisation,         -1, 2, 2);
+      pmanager->AddProcess(new G4hBremsstrahlung,     -1, 3, 3);
+      pmanager->AddProcess(new G4hPairProduction,     -1, 4, 4);       
 
     } else if (particleName == "alpha" ||
@@ -106,6 +121,6 @@
 	       particleName == "GenericIon") {
 
-      pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
-      pmanager->AddProcess(new G4ionIonisation,      -1, 2,2);
+      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+      pmanager->AddProcess(new G4ionIonisation,       -1, 2, 2);
      
     } else if ((!particle->IsShortLived()) &&
@@ -113,8 +128,39 @@
 	       (particle->GetParticleName() != "chargedgeantino")) {
 
-      pmanager->AddProcess(new G4MultipleScattering,-1,1,1);
-      pmanager->AddProcess(new G4hIonisation,       -1,2,2);
+      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
+      pmanager->AddProcess(new G4hIonisation,         -1, 2, 2);
     }
   }
+  
+  // Em options
+  //
+  // Main options and setting parameters are shown here.
+  // Several of them have default values.
+  //  
+  G4EmProcessOptions emOptions;
+  
+  //physics tables
+  //
+  emOptions.SetMinEnergy(100*eV);	//default    
+  emOptions.SetMaxEnergy(100*TeV);	//default  
+  emOptions.SetDEDXBinning(12*20);	//default=12*7  
+  emOptions.SetLambdaBinning(12*20);	//default=12*7
+  emOptions.SetSplineFlag(true);	//default  
+    
+  //coulomb scattering
+  //
+  emOptions.SetMscStepLimitation(fUseDistanceToBoundary);  //default=fUseSafety
+  emOptions.SetMscRangeFactor(0.04);	//default
+  emOptions.SetMscGeomFactor (2.5);	//default       
+  emOptions.SetSkin(3.);		//default
+          
+  //energy loss
+  //
+  emOptions.SetStepFunction(0.2, 100*um);	//default=(0.2, 1*mm)   
+  emOptions.SetLinearLossLimit(1.e-2);		//default
+   
+  //ionization
+  //
+  emOptions.SetSubCutoff(true);		//default=false  
 }
 
Index: trunk/examples/extended/electromagnetic/TestEm9/src/PhysicsList.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/PhysicsList.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/PhysicsList.cc	(revision 1230)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: PhysicsList.cc,v 1.19 2007/11/13 14:44:26 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: PhysicsList.cc,v 1.24 2008/10/16 11:48:58 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //---------------------------------------------------------------------------
@@ -46,8 +46,8 @@
 
 #include "PhysListEmStandard.hh"
-#include "PhysListEmStandardIG.hh"
 #include "G4EmStandardPhysics.hh"
 #include "G4EmStandardPhysics_option1.hh"
 #include "G4EmStandardPhysics_option2.hh"
+#include "G4EmStandardPhysics_option3.hh"
 #include "PhysListEmLivermore.hh"
 #include "PhysListEmPenelope.hh"
@@ -85,4 +85,6 @@
   cutForElectron  = defaultCutValue;
   cutForPositron  = defaultCutValue;
+  cutForVertexDetector = defaultCutValue;
+  cutForMuonDetector   = defaultCutValue;
 
   vertexDetectorCuts = 0;
@@ -96,5 +98,4 @@
   SetVerboseLevel(1);
 
-  mscStepLimit   = true;
   helIsRegisted  = false;
   bicIsRegisted  = false;
@@ -106,6 +107,6 @@
 
   // EM physics
-  emName = G4String("standard");
-  emPhysicsList = new PhysListEmStandard(emName);
+  emName = G4String("emstandard");
+  emPhysicsList = new G4EmStandardPhysics();
 }
 
@@ -173,5 +174,12 @@
       G4cout << "PhysicsList::Set " << name << " EM physics" << G4endl;
 
-  } else if (name == "standard") {
+  } else if (name == "emstandard_opt3") {
+    emName = name;
+    delete emPhysicsList;
+    emPhysicsList = new G4EmStandardPhysics_option3();
+    if (verboseLevel > 0) 
+      G4cout << "PhysicsList::Set " << name << " EM physics" << G4endl;
+
+  } else if (name == "emstandard_local") {
     emName = name;
     delete emPhysicsList;
@@ -179,11 +187,4 @@
     if (verboseLevel > 0) 
       G4cout << "PhysicsList::Set " << name << " EM physics" << G4endl;
-
-  } else if (name == "standardIG") {
-    emName = name;
-    delete emPhysicsList;
-    emPhysicsList = new PhysListEmStandardIG();
-    if (verboseLevel > 0) 
-      G4cout << "PhysicsList::Set StandardIG EM physics" << G4endl;
 
   } else if (name == "livermore") {
@@ -255,21 +256,23 @@
 void PhysicsList::SetCuts()
 {
-
   SetCutValue(cutForGamma, "gamma", "DefaultRegionForTheWorld");
   SetCutValue(cutForElectron, "e-", "DefaultRegionForTheWorld");
   SetCutValue(cutForPositron, "e+", "DefaultRegionForTheWorld");
-  G4cout << "PhysicsList: world cuts are set cutG= " << cutForGamma/mm 
-	 << " mm    cutE= " << cutForElectron/mm << " mm " << G4endl;
-
-  if( !vertexDetectorCuts ) SetVertexCut(cutForElectron);
+  //  G4cout << "PhysicsList: world cuts are set cutG= " << cutForGamma/mm 
+  //	 << " mm    cutE= " << cutForElectron/mm << " mm " << G4endl;
+
+  //G4cout << " cutV= " << cutForVertexDetector 
+  //     << " cutM= " << cutForMuonDetector<<G4endl;
+
   G4Region* region = (G4RegionStore::GetInstance())->GetRegion("VertexDetector");
-  region->SetProductionCuts(vertexDetectorCuts);
-  G4cout << "Vertex cuts are set" << G4endl;
-
-  if( !muonDetectorCuts ) SetMuonCut(cutForElectron);
+  vertexDetectorCuts = region->GetProductionCuts();
+  SetVertexCut(cutForVertexDetector);
+  //  G4cout << "Vertex cuts are set" << G4endl;
+ 
   region = (G4RegionStore::GetInstance())->GetRegion("MuonDetector");
-  region->SetProductionCuts(muonDetectorCuts);
-  G4cout << "Muon cuts are set" << G4endl;
-
+  muonDetectorCuts = region->GetProductionCuts();
+  SetMuonCut(cutForMuonDetector);
+  //G4cout << "Muon cuts are set " <<muonRegion << " " << muonDetectorCuts << G4endl;
+  
   if (verboseLevel>0) DumpCutValuesTable();
 }
@@ -303,10 +306,11 @@
 void PhysicsList::SetVertexCut(G4double cut)
 {
-  if( !vertexDetectorCuts ) vertexDetectorCuts = new G4ProductionCuts();
-
-  vertexDetectorCuts->SetProductionCut(cut, idxG4GammaCut);
-  vertexDetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
-  vertexDetectorCuts->SetProductionCut(cut, idxG4PositronCut);
-
+  cutForVertexDetector = cut;
+  
+  if( vertexDetectorCuts ) {
+    vertexDetectorCuts->SetProductionCut(cut, idxG4GammaCut);
+    vertexDetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
+    vertexDetectorCuts->SetProductionCut(cut, idxG4PositronCut);
+  } 
 }
 
@@ -315,18 +319,13 @@
 void PhysicsList::SetMuonCut(G4double cut)
 {
-  if( !muonDetectorCuts ) muonDetectorCuts = new G4ProductionCuts();
-
-  muonDetectorCuts->SetProductionCut(cut, idxG4GammaCut);
-  muonDetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
-  muonDetectorCuts->SetProductionCut(cut, idxG4PositronCut);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void PhysicsList::SetMscStepLimit(G4bool val)
-{
-  mscStepLimit = val;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
+  cutForMuonDetector   = cut;
+
+  if( muonDetectorCuts ) {
+    muonDetectorCuts->SetProductionCut(cut, idxG4GammaCut);
+    muonDetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
+    muonDetectorCuts->SetProductionCut(cut, idxG4PositronCut);
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
Index: trunk/examples/extended/electromagnetic/TestEm9/src/PhysicsListMessenger.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/PhysicsListMessenger.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/PhysicsListMessenger.cc	(revision 1230)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: PhysicsListMessenger.cc,v 1.3 2006/06/29 17:03:43 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: PhysicsListMessenger.cc,v 1.4 2008/04/03 15:07:55 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // 
@@ -143,6 +143,6 @@
    { pPhysicsList->AddPhysicsList(newValue);}
 
-  if( command == mscCmd )
-   { pPhysicsList->SetMscStepLimit(mscCmd->GetNewBoolValue(newValue));}
+  //  if( command == mscCmd )
+  // { pPhysicsList->SetMscStepLimit(mscCmd->GetNewBoolValue(newValue));}
 }
 
Index: trunk/examples/extended/electromagnetic/TestEm9/src/PrimaryGeneratorAction.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/PrimaryGeneratorAction.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/PrimaryGeneratorAction.cc	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: PrimaryGeneratorAction.cc,v 1.3 2006/06/29 17:03:46 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //
Index: trunk/examples/extended/electromagnetic/TestEm9/src/StepMax.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/StepMax.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/StepMax.cc	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: StepMax.cc,v 1.3 2006/06/29 17:03:52 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm9/src/StepMaxMessenger.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/StepMaxMessenger.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/StepMaxMessenger.cc	(revision 1230)
@@ -25,5 +25,5 @@
 //
 // $Id: StepMaxMessenger.cc,v 1.2 2006/06/29 17:03:54 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/examples/extended/electromagnetic/TestEm9/src/SteppingAction.cc
===================================================================
--- trunk/examples/extended/electromagnetic/TestEm9/src/SteppingAction.cc	(revision 812)
+++ trunk/examples/extended/electromagnetic/TestEm9/src/SteppingAction.cc	(revision 1230)
@@ -26,5 +26,5 @@
 //
 // $Id: SteppingAction.cc,v 1.2 2006/06/29 17:03:57 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //
