Ignore:
Timestamp:
Jan 8, 2010, 3:02:48 PM (14 years ago)
Author:
garnier
Message:

update to geant4.9.3

Location:
trunk/examples/extended/electromagnetic/TestEm9/src
Files:
17 edited

Legend:

Unmodified
Added
Removed

  • trunk/examples/extended/electromagnetic/TestEm9/src/DetectorConstruction.cc

    r807 r1230  
    2525//
    2626//
    27 // $Id: DetectorConstruction.cc,v 1.9 2006/06/29 17:02:59 gunter Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: DetectorConstruction.cc,v 1.11 2008/04/07 18:09:05 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    5858#include "G4Region.hh"
    5959#include "G4RegionStore.hh"
     60#include "G4ProductionCuts.hh"
    6061#include "G4PhysicalVolumeStore.hh"
    6162#include "G4LogicalVolumeStore.hh"
    … …  
    7273
    7374DetectorConstruction::DetectorConstruction()
    74   :logicC(0),logicA1(0),logicA2(0)
     75  : G4VUserDetectorConstruction()
    7576{
    7677  detectorMessenger = new DetectorMessenger(this);
    … …  
    8485  vertexRegion = 0;
    8586  muonRegion   = 0;
     87  logicWorld   = 0;
     88  logicCal     = 0;
     89  logicA1      = 0;
    8690  DefineMaterials();
     91  vertexDetectorCuts = new G4ProductionCuts();
     92  muonDetectorCuts   = new G4ProductionCuts();
    8793}
    8894
    … …  
    9096
    9197DetectorConstruction::~DetectorConstruction()
    92 { delete detectorMessenger;}
     98{
     99  delete detectorMessenger;
     100  delete vertexDetectorCuts;
     101  delete muonDetectorCuts;
     102}
    93103
    94104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    106116
    107117  G4NistManager* man = G4NistManager::Instance();
    108   man->SetVerbose(1);
     118  //  man->SetVerbose(1);
    109119  worldMaterial = man->FindOrBuildMaterial("G4_AIR");
    110120  absMaterial   = man->FindOrBuildMaterial("G4_Al");
    … …  
    118128G4VPhysicalVolume* DetectorConstruction::ConstructVolumes()
    119129{
    120   G4cout << *(G4Material::GetMaterialTable()) << G4endl;
    121 
    122130  // Cleanup old geometry
    123131
    124132  G4GeometryManager::GetInstance()->OpenGeometry();
     133
     134  if(G4NistManager::Instance()->GetVerbose() > 0)
     135    G4cout << *(G4Material::GetMaterialTable()) << G4endl;
     136
     137  if(vertexRegion) {
     138    delete vertexRegion;
     139    delete muonRegion;
     140  }
     141  vertexRegion = new G4Region("VertexDetector");
     142  vertexRegion->SetProductionCuts(vertexDetectorCuts);
     143
     144  muonRegion   = new G4Region("MuonDetector");
     145  muonRegion->SetProductionCuts(muonDetectorCuts);
     146
     147  G4SolidStore::GetInstance()->Clean();
     148  G4LogicalVolumeStore::GetInstance()->Clean();
    125149  G4PhysicalVolumeStore::GetInstance()->Clean();
    126   G4LogicalVolumeStore::GetInstance()->Clean();
    127   G4SolidStore::GetInstance()->Clean();
    128 
    129   if(vertexRegion) delete vertexRegion;
    130   if(muonRegion) delete muonRegion;
    131   vertexRegion = new G4Region("VertexDetector");
    132   muonRegion   = new G4Region("MuonDetector");
    133150
    134151  if(vertexLength < padLength*5.0) vertexLength = padLength*5.0;
    … …  
    137154  G4double york   = 10.*cm;
    138155
    139            worldZ = 2.*vertexLength + 3.*absLength + 0.5*(ecalLength + york) + biggap*2.;
     156  worldZ = 2.*vertexLength + 3.*absLength + 0.5*(ecalLength + york) + biggap*2.;
     157
    140158  G4double worldX = ecalWidth*3.0;
    141159  G4double vertexZ= -worldZ + vertexLength*2.0 + absLength     + biggap;
    142160  G4double absZ2  = -worldZ + vertexLength*4.0 + absLength*3.5 + biggap;
    143   G4double ecalZ  = -worldZ + vertexLength*4.0 + absLength*4.0 + ecalLength*0.5 + 2.*biggap;
     161  G4double ecalZ  = -worldZ + vertexLength*4.0 + absLength*4.0 + ecalLength*0.5
     162    + 2.*biggap;
    144163  G4double yorkZ  = -worldZ + vertexLength*4.0 + absLength*5.0 + ecalLength
    145                             + york*0.5 + 3.*biggap;
     164    + york*0.5 + 3.*biggap;
    146165
    147166  //
    … …  
    149168  //
    150169  G4Box* solidW = new G4Box("World",worldX,worldX,worldZ);
    151   G4LogicalVolume* logicW = new G4LogicalVolume( solidW,worldMaterial,
    152                                                 "World");
     170  logicWorld = new G4LogicalVolume( solidW,worldMaterial,"World");
    153171  G4VPhysicalVolume* world = new G4PVPlacement(0,G4ThreeVector(),
    154                                        "World",logicW,0,false,0);
     172                                               "World",logicWorld,0,false,0);
    155173
    156174  //
    … …  
    158176  //
    159177  G4Box* solidE = new G4Box("VolE",worldX,worldX,ecalLength*0.5 + gap);
    160   G4LogicalVolume* logicE = new G4LogicalVolume( solidE,worldMaterial,
    161                                                 "VolE");
     178  logicECal = new G4LogicalVolume( solidE,worldMaterial,"VolE");
    162179  G4VPhysicalVolume* physE = new G4PVPlacement(0,G4ThreeVector(0.,0.,ecalZ),
    163                                        "VolE",logicE,world,false,0);
     180                                               "VolE",logicECal,world,false,0);
    164181
    165182  G4Box* solidC = new G4Box("Ecal",ecalWidth*0.5,ecalWidth*0.5,ecalLength*0.5);
    166   logicC = new G4LogicalVolume( solidC,calMaterial,"Ecal");
     183  logicCal = new G4LogicalVolume( solidC,calMaterial,"Ecal");
    167184
    168185  G4cout << "Ecal is " << G4BestUnit(ecalLength,"Length")
    169       << " of " << calMaterial->GetName() << G4endl;
     186        << " of " << calMaterial->GetName() << G4endl;
    170187
    171188  // Crystals
    … …  
    182199    for (j=0; j<5; j++) {
    183200
    184       pv = new G4PVPlacement(0,G4ThreeVector(x,y,0.),"Ecal",logicC,
     201      pv = new G4PVPlacement(0,G4ThreeVector(x,y,0.),"Ecal",logicCal,
    185202                                    physE,false,k);
    186203      k++;
    … …  
    195212  logicA2 = new G4LogicalVolume( solidA,absMaterial,"Abs2");
    196213  pv = new G4PVPlacement(0,G4ThreeVector(0.,0.,absZ2),
    197                                       "Abs2",logicA2,world,false,0);
     214                        "Abs2",logicA2,world,false,0);
    198215
    199216  G4cout << "Absorber is " << G4BestUnit(absLength,"Length")
    200       << " of " << absMaterial->GetName() << G4endl;
     217        << " of " << absMaterial->GetName() << G4endl;
    201218
    202219  //York
    203220
    204221  G4Box* solidYV = new G4Box("VolY",worldX,worldX,york*0.5+absLength);
    205   G4LogicalVolume* logicYV = new G4LogicalVolume( solidYV,yorkMaterial,"VolY");
     222  logicYV = new G4LogicalVolume( solidYV,yorkMaterial,"VolY");
    206223  G4VPhysicalVolume* physYV = new G4PVPlacement(0,G4ThreeVector(0.,0.,yorkZ),
    207                                        "VolY",logicYV,world,false,0);
     224                                                "VolY",logicYV,world,false,0);
    208225
    209226  G4Box* solidY = new G4Box("York",worldX,worldX,york*0.5);
    210   G4LogicalVolume* logicY = new G4LogicalVolume( solidY,yorkMaterial,"York");
     227  logicY = new G4LogicalVolume( solidY,yorkMaterial,"York");
    211228  pv = new G4PVPlacement(0,G4ThreeVector(),
    212                                       "York",logicY,physYV,false,0);
     229                        "York",logicY,physYV,false,0);
    213230
    214231  logicA3 = new G4LogicalVolume( solidA,absMaterial,"Abs3");
    215232  logicA4 = new G4LogicalVolume( solidA,absMaterial,"Abs4");
     233
    216234  pv = new G4PVPlacement(0,G4ThreeVector(0.,0.,-(york+absLength)*0.5),
    217                                       "Abs3",logicA3,physYV,false,0);
     235                        "Abs3",logicA3,physYV,false,0);
    218236  pv = new G4PVPlacement(0,G4ThreeVector(0.,0.,(york+absLength)*0.5),
    219                                       "Abs4",logicA4,physYV,false,0);
     237                        "Abs4",logicA4,physYV,false,0);
    220238
    221239  //Vertex volume
    222 
    223240  G4Box* solidVV = new G4Box("VolV",worldX,worldX,vertexLength*2.+absLength+gap);
    224   G4LogicalVolume* logicVV = new G4LogicalVolume( solidVV,worldMaterial,"VolV");
     241  logicVV = new G4LogicalVolume( solidVV,worldMaterial,"VolV");
    225242  G4VPhysicalVolume* physVV = new G4PVPlacement(0,G4ThreeVector(0.,0.,vertexZ),
    226                                        "VolV",logicVV,world,false,0);
     243                                                "VolV",logicVV,world,false,0);
    227244
    228245  //Absorber
    229 
    230246  logicA1 = new G4LogicalVolume( solidA,absMaterial,"Abs1");
    231247  pv = new G4PVPlacement(0,G4ThreeVector(0.,0.,vertexLength*2.-absLength*0.5),
    232                                       "Abs1",logicA1,physVV,false,0);
     248                        "Abs1",logicA1,physVV,false,0);
    233249
    234250  //Vertex
    235 
    236251  G4double vertWidth = ecalWidth/5.;
    237252  G4int npads = (G4int)(vertWidth/padWidth);
    … …  
    242257
    243258  G4Box* solidVD = new G4Box("VertDet",x1,ecalWidth*0.5+gap,padLength*0.5);
    244   G4LogicalVolume* logicVD = new G4LogicalVolume( solidVD,vertMaterial,"VertDet");
     259  logicVD = new G4LogicalVolume( solidVD,vertMaterial,"VertDet");
     260  logicVD->SetSolid(solidVD);
    245261
    246262  G4Box* solidV = new G4Box("Vert",padWidth*0.5,ecalWidth*0.5,padLength*0.5);
    247   G4LogicalVolume* logicV = new G4LogicalVolume( solidV,vertMaterial,"Vert");
     263  logicV = new G4LogicalVolume( solidV,vertMaterial,"Vert");
    248264
    249265  for (i=0; i<3; i++) {
    … …  
    267283
    268284  // Define region for the vertex detector
    269 
    270285  vertexRegion->AddRootLogicalVolume(logicVV);
    271286  vertexRegion->AddRootLogicalVolume(logicA3);
    272287
    273288  // Define region for the muon detector
    274 
    275289  muonRegion->AddRootLogicalVolume(logicYV);
    276290
    277291  // color regions
    278 
    279292  logicVV-> SetVisAttributes(G4VisAttributes::Invisible);
    280293  logicV-> SetVisAttributes(G4VisAttributes::Invisible);
    281   logicE-> SetVisAttributes(G4VisAttributes::Invisible);
     294  logicECal-> SetVisAttributes(G4VisAttributes::Invisible);
    282295  logicYV-> SetVisAttributes(G4VisAttributes::Invisible);
    283296
    284297  G4VisAttributes* regWcolor = new G4VisAttributes(G4Colour(0.3, 0.3, 0.3));
    285   logicW->SetVisAttributes(regWcolor);
     298  logicWorld->SetVisAttributes(regWcolor);
    286299
    287300  G4VisAttributes* regVcolor = new G4VisAttributes(G4Colour(0., 0.3, 0.7));
    … …  
    289302
    290303  G4VisAttributes* regCcolor = new G4VisAttributes(G4Colour(0., 0.7, 0.3));
    291   logicC->SetVisAttributes(regCcolor);
     304  logicCal->SetVisAttributes(regCcolor);
    292305
    293306  G4VisAttributes* regAcolor = new G4VisAttributes(G4Colour(1., 0.5, 0.5));
    … …  
    302315  // always return world
    303316  G4cout << "### New geometry is constructed" << G4endl;
    304 
     317 
    305318  return world;
    306319}
    … …  
    313326  G4Material* pttoMaterial =
    314327    G4NistManager::Instance()->FindOrBuildMaterial(mat, false);
    315   if (pttoMaterial) calMaterial = pttoMaterial;
     328  if (pttoMaterial) {
     329    calMaterial = pttoMaterial;
     330    if(logicCal) {
     331      logicCal->SetMaterial(calMaterial);
     332      G4RunManager::GetRunManager()->PhysicsHasBeenModified();
     333    }
     334  }
    316335}
    317336
    … …  
    323342  G4Material* pttoMaterial =
    324343    G4NistManager::Instance()->FindOrBuildMaterial(mat, false);
    325   if (pttoMaterial) absMaterial = pttoMaterial;
     344  if (pttoMaterial) {
     345    absMaterial = pttoMaterial;
     346    if(logicA1) {
     347      logicA1->SetMaterial(absMaterial);
     348      logicA2->SetMaterial(absMaterial);
     349      logicA3->SetMaterial(absMaterial);
     350      logicA4->SetMaterial(absMaterial);
     351      G4RunManager::GetRunManager()->PhysicsHasBeenModified();
     352    }
     353  }
    326354}
    327355
    … …  
    330358void DetectorConstruction::UpdateGeometry()
    331359{
     360  G4RunManager::GetRunManager()->PhysicsHasBeenModified();
    332361  G4RunManager::GetRunManager()->DefineWorldVolume(ConstructVolumes());
    333362}
    334363
    335364//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     365
     366void DetectorConstruction::SetEcalLength (G4double val)   
     367{
     368  if(val > 0.0) {
     369    ecalLength = val;
     370    G4RunManager::GetRunManager()->GeometryHasBeenModified();
     371  }
     372}
     373
     374//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     375
     376void DetectorConstruction::SetEcalWidth  (G4double val)   
     377{
     378  if(val > 0.0) {
     379    ecalWidth = val;
     380    G4RunManager::GetRunManager()->GeometryHasBeenModified();
     381  }
     382}
     383
     384//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     385
     386void DetectorConstruction::SetVertexLength (G4double val)
     387{
     388  if(val > 0.0) {
     389    vertexLength = val;
     390    G4RunManager::GetRunManager()->GeometryHasBeenModified();
     391  }
     392}
     393
     394//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     395
     396void DetectorConstruction::SetPadLength  (G4double val)   
     397{
     398  if(val > 0.0) {
     399    padLength = val;
     400    G4RunManager::GetRunManager()->GeometryHasBeenModified();
     401  }
     402}
     403
     404//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     405
     406void DetectorConstruction::SetPadWidth  (G4double val)   
     407{
     408  if(val > 0.0) {
     409    padWidth = val;
     410    G4RunManager::GetRunManager()->GeometryHasBeenModified();
     411  }
     412}
     413
     414//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     415
     416void DetectorConstruction::SetAbsLength(G4double val)     
     417{
     418  if(val > 0.0) {
     419    absLength = val;
     420    G4RunManager::GetRunManager()->GeometryHasBeenModified();
     421  }
     422}
     423
     424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm9/src/DetectorMessenger.cc

    r807 r1230  
    2626//
    2727// $Id: DetectorMessenger.cc,v 1.3 2006/06/29 17:03:01 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//

  • trunk/examples/extended/electromagnetic/TestEm9/src/EmAcceptance.cc

    r807 r1230  
    2626//
    2727// $Id: EmAcceptance.cc,v 1.4 2006/06/29 17:03:03 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//

  • trunk/examples/extended/electromagnetic/TestEm9/src/EventAction.cc

    r807 r1230  
    2525//
    2626// $Id: EventAction.cc,v 1.4 2006/06/29 17:03:06 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

  • trunk/examples/extended/electromagnetic/TestEm9/src/EventActionMessenger.cc

    r807 r1230  
    2525//
    2626// $Id: EventActionMessenger.cc,v 1.3 2006/06/29 17:03:08 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm9/src/Histo.cc

    r807 r1230  
    5353  verbose    = 0;
    5454  histName   = "testem9";
    55   histType   = "hbook";
     55  histType   = "root";
    5656  nHisto     = 0;
    5757  defaultAct = 1;
    … …  
    9393
    9494  G4String nam = histName + "." + histType;
    95 
    96   tree = tf->create(nam,histType,false,true,"uncompress");
     95  G4String options  = "--noErrors export=root uncompress";
     96 
     97  tree = tf->create(nam,histType,false,true,options);
    9798  delete tf;
    9899  if(tree) {
    … …  
    107108  // Creating an 1-dimensional histograms in the root directory of the tree
    108109  for(G4int i=0; i<nHisto; i++) {
    109     if(active[i])
    110       histo[i] = hf->createHistogram1D(ids[i], titles[i], bins[i], xmin[i], xmax[i]);
     110    if(active[i]) {
     111      G4String ss = ids[i];
     112      if(histType == "root") ss = "h" + ids[i];
     113      histo[i] = hf->createHistogram1D(ss, titles[i], bins[i], xmin[i], xmax[i]);
     114    }
    111115  }
    112116  delete hf;

  • trunk/examples/extended/electromagnetic/TestEm9/src/HistoManager.cc

    r807 r1230  
    8787  nBinsED= 100;
    8888  nTuple = false;
    89   nHisto = 10;
     89  nHisto = 13;
    9090
    9191  // initialise acceptance
    … …  
    9797
    9898  histo->add1D("10",
    99     "Energy deposit (MeV) in central crystal",nBinsED,0.0,beamEnergy,MeV);
     99    "Evis/E0 in central crystal",nBinsED,0.0,1,1.0);
    100100
    101101  histo->add1D("11",
    102     "Energy deposit (MeV) in 3x3",nBinsED,0.0,beamEnergy,MeV);
     102    "Evis/E0 in 3x3",nBinsED,0.0,1.0,1.0);
    103103
    104104  histo->add1D("12",
    105     "Energy deposit (MeV) in 5x5",nBinsED,0.0,beamEnergy,MeV);
     105    "Evis/E0 in 5x5",nBinsED,0.0,1.0,1.0);
    106106
    107107  histo->add1D("13",
    … …  
    125125  histo->add1D("19",
    126126    "Number of vertex hits",20,-0.5,19.5,1.0);
     127
     128  histo->add1D("20",
     129    "E1/E9 Ratio",nBinsED,0.0,1,1.0);
     130
     131  histo->add1D("21",
     132    "E1/25 Ratio",nBinsED,0.0,1.0,1.0);
     133
     134  histo->add1D("22",
     135    "E9/E25 Ratio",nBinsED,0.0,1.0,1.0);
    127136
    128137  if(nTuple) {
    … …  
    141150  n_gam  = 0;
    142151  n_step = 0;
    143 
    144   for(G4int i=0; i<nmax; i++) {
     152  n_lowe = 0;
     153
     154  for(G4int i=0; i<6; i++) {
    145155    stat[i] = 0;
    146156    edep[i] = 0.0;
    147157    erms[i] = 0.0;
    148     edeptr[i] = 0.0;
    149     ermstr[i] = 0.0;
     158    if(i < 3) {
     159      edeptr[i] = 0.0;
     160      ermstr[i] = 0.0;
     161    }
    150162  }
    151163
    … …  
    164176
    165177  G4cout << "HistoManager: End of run actions are started" << G4endl;
    166   G4String nam[3] = {"1x1", "3x3", "5x5"};
     178  G4String nam[6] = {"1x1", "3x3", "5x5", "E1/E9 ", "E1/E25", "E9/E25"};
    167179
    168180  // average
    … …  
    175187
    176188    // total mean
    177     edep[j] *= x/beamEnergy;
    178     G4double y = erms[j]*x/(beamEnergy*beamEnergy) - edep[j]*edep[j];
     189    edep[j] *= x;
     190    G4double y = erms[j]*x - edep[j]*edep[j];
    179191    if(y < 0.0) y = 0.0;
    180192    erms[j] = std::sqrt(y);
    … …  
    183195    G4double xx = G4double(stat[j]);
    184196    if(xx > 0.0) xx = 1.0/xx;
    185     edeptr[j] *= xx/beamEnergy;
    186     y = ermstr[j]*xx/(beamEnergy*beamEnergy) - edeptr[j]*edeptr[j];
     197    edeptr[j] *= xx;
     198    y = ermstr[j]*xx - edeptr[j]*edeptr[j];
    187199    if(y < 0.0) y = 0.0;
    188200    ermstr[j] = std::sqrt(y);
    … …  
    191203  G4double xg = x*(G4double)n_gam;
    192204  G4double xp = x*(G4double)n_posit;
    193   G4double xs = x*(G4double)n_step;
     205  G4double xs = x*n_step;
    194206
    195207  G4double f = 100.*std::sqrt(beamEnergy/GeV);
    … …  
    201213  G4cout << std::setprecision(4) << "Average number of steps      " << xs << G4endl;
    202214 
    203   for(j=0; j<nmax; j++) {
     215  for(j=0; j<3; j++) {
    204216    G4double e = edeptr[j];
    205217    G4double s = ermstr[j];
    206     G4double r = s*std::sqrt(x);
     218    G4double xx= G4double(stat[j]);
     219    if(xx > 0.0) xx = 1.0/xx;
     220    G4double r = s*std::sqrt(xx);
    207221    G4cout << std::setprecision(4) << "Edep " << nam[j] << " =                   " << e
    208222           << " +- " << r;
    … …  
    222236    }
    223237  }
     238  G4cout<<"===========  Ratios without trancating ==========================="<<G4endl;
     239  for(j=3; j<6; j++) {
     240    G4double e = edep[j];
     241    G4double xx= G4double(stat[j]);
     242    if(xx > 0.0) xx = 1.0/xx;
     243    e *= xx;
     244    G4double y = erms[j]*xx - e*e;
     245    G4double r = 0.0;
     246    if(y > 0.0) r = std::sqrt(y*xx);
     247    G4cout << "  " << nam[j] << " =                   " << e
     248           << " +- " << r;
     249    G4cout << G4endl;
     250  }
     251  G4cout << std::setprecision(4) << "Beam Energy                  " << beamEnergy/GeV
     252         << " GeV" << G4endl;
     253  if(n_lowe > 0)          G4cout << "Number of events E/E0<0.8    " << n_lowe << G4endl;
    224254  G4cout<<"=================================================================="<<G4endl;
    225255  G4cout<<G4endl;
    … …  
    264294  Evertex.clear();
    265295  Nvertex.clear();
    266   for (int i=0; i<25; i++) {
     296  for (G4int i=0; i<25; i++) {
    267297    E[i] = 0.0;
    268298  }
    … …  
    275305  G4double e9 = 0.0;
    276306  G4double e25= 0.0;
    277   for (int i=0; i<25; i++) {
     307  for (G4int i=0; i<25; i++) {
     308    E[i] /= beamEnergy;
    278309    e25 += E[i];
    279310    if( ( 6<=i &&  8>=i) || (11<=i && 13>=i) || (16<=i && 18>=i)) e9 += E[i];
    280311  }
    281   histo->fill(0,E[12],1.0);
     312
     313  if(e25 < 0.8) {
     314    n_lowe++;
     315    G4cout << "### in the event# " << n_evt << "  E25= " << e25 << G4endl;
     316  }
     317
     318  // compute ratios
     319  G4double e0 = E[12];
     320  G4double e19  = 0.0;
     321  G4double e125 = 0.0;
     322  G4double e925 = 0.0;
     323  if(e9 > 0.0) {
     324    e19 = e0/e9;
     325    e125 = e0/e25;
     326    e925 = e9/e25;
     327    edep[3] += e19;
     328    erms[3] += e19*e19;
     329    edep[4] += e125;
     330    erms[4] += e125*e125;
     331    edep[5] += e925;
     332    erms[5] += e925*e925;
     333    stat[3] += 1;
     334    stat[4] += 1;
     335    stat[5] += 1;
     336  }
     337
     338  // fill histo
     339  histo->fill(0,e0,1.0);
    282340  histo->fill(1,e9,1.0);
    283341  histo->fill(2,e25,1.0);
    … …  
    286344  histo->fill(7,Eabs3,1.0);
    287345  histo->fill(8,Eabs4,1.0);
    288   float nn = (double)(Nvertex.size());
    289   histo->fill(9,nn,1.0);
    290 
    291   G4double e0 = E[12];
    292 
     346  histo->fill(9,G4double(Nvertex.size()),1.0);
     347  histo->fill(10,e19,1.0);
     348  histo->fill(11,e125,1.0);
     349  histo->fill(12,e925,1.0);
     350
     351  // compute sums
    293352  edep[0] += e0;
    294353  erms[0] += e0*e0;
    … …  
    299358
    300359  // trancated mean
    301   if(limittrue[0] == DBL_MAX || std::abs(e0/beamEnergy-edeptrue[0])<rmstrue[0]*limittrue[0]) {
     360  if(limittrue[0] == DBL_MAX || std::abs(e0-edeptrue[0])<rmstrue[0]*limittrue[0]) {
    302361    stat[0] += 1;
    303362    edeptr[0] += e0;
    304363    ermstr[0] += e0*e0;
    305364  }
    306   if(limittrue[1] == DBL_MAX || std::abs(e9/beamEnergy-edeptrue[1])<rmstrue[1]*limittrue[1]) {
     365  if(limittrue[1] == DBL_MAX || std::abs(e9-edeptrue[1])<rmstrue[1]*limittrue[1]) {
    307366    stat[1] += 1;
    308367    edeptr[1] += e9;
    309368    ermstr[1] += e9*e9;
    310369  }
    311   if(limittrue[2] == DBL_MAX || std::abs(e25/beamEnergy-edeptrue[2])<rmstrue[2]*limittrue[2]) {
     370  if(limittrue[2] == DBL_MAX || std::abs(e25-edeptrue[2])<rmstrue[2]*limittrue[2]) {
    312371    stat[2] += 1;
    313372    edeptr[2] += e25;
    … …  
    333392  if(0 == pid) {
    334393
     394    beamEnergy = kinE;
    335395    histo->fillTuple("TKIN", kinE/MeV);
    336396

  • trunk/examples/extended/electromagnetic/TestEm9/src/HistoMessenger.cc

    r807 r1230  
    2525//
    2626//
    27 // $Id: HistoMessenger.cc,v 1.6 2006/06/29 17:03:16 gunter Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: HistoMessenger.cc,v 1.7 2008/08/22 14:11:53 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    6060  //
    6161  G4UIparameter* ih = new G4UIparameter("ih",'i',false);
    62   ih->SetGuidance("histo number : from 1 to MaxHisto");
     62  ih->SetGuidance("histo number : from 0 to MaxHisto-1");
    6363  histoCmd->SetParameter(ih);
    6464  //
    … …  
    110110     G4String unit = unts;
    111111     G4double vUnit = 1. ;
    112      if (unit != "none") vUnit = G4UIcommand::ValueOf(unit);
     112     if(unit != "none") vUnit = G4UIcommand::ValueOf(unit);
     113     if(vUnit <= 0.0)  vUnit = 1.;
    113114     histo->setHisto1D(ih,nbBins,vmin,vmax,vUnit);
    114115   }

  • trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmLivermore.cc

    r807 r1230  
    2626//
    2727// $Id: PhysListEmLivermore.cc,v 1.2 2007/01/08 16:29:42 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmPenelope.cc

    r807 r1230  
    2626//
    2727// $Id: PhysListEmPenelope.cc,v 1.1 2006/11/17 17:45:57 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm9/src/PhysListEmStandard.cc

    r807 r1230  
    2525//
    2626//
    27 // $Id: PhysListEmStandard.cc,v 1.8 2006/11/22 19:09:12 vnivanch Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: PhysListEmStandard.cc,v 1.13 2008/11/16 21:01:10 maire Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    3939#include "G4PhotoElectricEffect.hh"
    4040
    41 #include "G4MultipleScattering.hh"
     41#include "G4eMultipleScattering.hh"
     42#include "G4hMultipleScattering.hh"
    4243
    4344#include "G4eIonisation.hh"
    … …  
    4849#include "G4MuBremsstrahlung.hh"
    4950#include "G4MuPairProduction.hh"
     51#include "G4hBremsstrahlung.hh"
     52#include "G4hPairProduction.hh"
    5053
    5154#include "G4hIonisation.hh"
    5255#include "G4ionIonisation.hh"
     56
     57#include "G4EmProcessOptions.hh"
     58#include "G4MscStepLimitType.hh"
    5359
    5460//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    8389    } else if (particleName == "e-") {
    8490 
    85       pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
    86       pmanager->AddProcess(new G4eIonisation(),      -1, 2,2);
    87       pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3,3);
     91      pmanager->AddProcess(new G4eMultipleScattering, -1, 1, 1);
     92      pmanager->AddProcess(new G4eIonisation(),       -1, 2, 2);
     93      pmanager->AddProcess(new G4eBremsstrahlung,     -1, 3, 3);
    8894           
    8995    } else if (particleName == "e+") {
    9096
    91       pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
    92       pmanager->AddProcess(new G4eIonisation(),      -1, 2,2);
    93       pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3,3);
    94       pmanager->AddProcess(new G4eplusAnnihilation,   0,-1,4);
     97      pmanager->AddProcess(new G4eMultipleScattering, -1, 1, 1);
     98      pmanager->AddProcess(new G4eIonisation(),       -1, 2, 2);
     99      pmanager->AddProcess(new G4eBremsstrahlung,     -1, 3, 3);
     100      pmanager->AddProcess(new G4eplusAnnihilation,    0,-1, 4);
    95101     
    96102    } else if (particleName == "mu+" ||
    97103               particleName == "mu-"    ) {
    98104
    99       pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
    100       pmanager->AddProcess(new G4MuIonisation,      -1, 2,2);
    101       pmanager->AddProcess(new G4MuBremsstrahlung,  -1, 3,3);
    102       pmanager->AddProcess(new G4MuPairProduction,  -1, 4,4);       
     105      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
     106      pmanager->AddProcess(new G4MuIonisation,        -1, 2, 2);
     107      pmanager->AddProcess(new G4MuBremsstrahlung,    -1, 3, 3);
     108      pmanager->AddProcess(new G4MuPairProduction,    -1, 4, 4);       
     109
     110    } else if (particleName == "proton" ||
     111               particleName == "pi-" || 
     112               particleName == "pi+") {
     113
     114      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
     115      pmanager->AddProcess(new G4hIonisation,         -1, 2, 2);
     116      pmanager->AddProcess(new G4hBremsstrahlung,     -1, 3, 3);
     117      pmanager->AddProcess(new G4hPairProduction,     -1, 4, 4);       
    103118
    104119    } else if (particleName == "alpha" ||
    … …  
    106121               particleName == "GenericIon") {
    107122
    108       pmanager->AddProcess(new G4MultipleScattering, -1, 1,1);
    109       pmanager->AddProcess(new G4ionIonisation,      -1, 2,2);
     123      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
     124      pmanager->AddProcess(new G4ionIonisation,       -1, 2, 2);
    110125     
    111126    } else if ((!particle->IsShortLived()) &&
    … …  
    113128               (particle->GetParticleName() != "chargedgeantino")) {
    114129
    115       pmanager->AddProcess(new G4MultipleScattering,-1,1,1);
    116       pmanager->AddProcess(new G4hIonisation,       -1,2,2);
     130      pmanager->AddProcess(new G4hMultipleScattering, -1, 1, 1);
     131      pmanager->AddProcess(new G4hIonisation,         -1, 2, 2);
    117132    }
    118133  }
     134 
     135  // Em options
     136  //
     137  // Main options and setting parameters are shown here.
     138  // Several of them have default values.
     139  // 
     140  G4EmProcessOptions emOptions;
     141 
     142  //physics tables
     143  //
     144  emOptions.SetMinEnergy(100*eV);       //default   
     145  emOptions.SetMaxEnergy(100*TeV);      //default 
     146  emOptions.SetDEDXBinning(12*20);      //default=12*7 
     147  emOptions.SetLambdaBinning(12*20);    //default=12*7
     148  emOptions.SetSplineFlag(true);        //default 
     149   
     150  //coulomb scattering
     151  //
     152  emOptions.SetMscStepLimitation(fUseDistanceToBoundary);  //default=fUseSafety
     153  emOptions.SetMscRangeFactor(0.04);    //default
     154  emOptions.SetMscGeomFactor (2.5);     //default       
     155  emOptions.SetSkin(3.);                //default
     156         
     157  //energy loss
     158  //
     159  emOptions.SetStepFunction(0.2, 100*um);       //default=(0.2, 1*mm)   
     160  emOptions.SetLinearLossLimit(1.e-2);          //default
     161   
     162  //ionization
     163  //
     164  emOptions.SetSubCutoff(true);         //default=false 
    119165}
    120166

  • trunk/examples/extended/electromagnetic/TestEm9/src/PhysicsList.cc

    r807 r1230  
    2424// ********************************************************************
    2525//
    26 // $Id: PhysicsList.cc,v 1.19 2007/11/13 14:44:26 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-01-patch-02 $
     26// $Id: PhysicsList.cc,v 1.24 2008/10/16 11:48:58 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//---------------------------------------------------------------------------
    … …  
    4646
    4747#include "PhysListEmStandard.hh"
    48 #include "PhysListEmStandardIG.hh"
    4948#include "G4EmStandardPhysics.hh"
    5049#include "G4EmStandardPhysics_option1.hh"
    5150#include "G4EmStandardPhysics_option2.hh"
     51#include "G4EmStandardPhysics_option3.hh"
    5252#include "PhysListEmLivermore.hh"
    5353#include "PhysListEmPenelope.hh"
    … …  
    8585  cutForElectron  = defaultCutValue;
    8686  cutForPositron  = defaultCutValue;
     87  cutForVertexDetector = defaultCutValue;
     88  cutForMuonDetector   = defaultCutValue;
    8789
    8890  vertexDetectorCuts = 0;
    … …  
    9698  SetVerboseLevel(1);
    9799
    98   mscStepLimit   = true;
    99100  helIsRegisted  = false;
    100101  bicIsRegisted  = false;
    … …  
    106107
    107108  // EM physics
    108   emName = G4String("standard");
    109   emPhysicsList = new PhysListEmStandard(emName);
     109  emName = G4String("emstandard");
     110  emPhysicsList = new G4EmStandardPhysics();
    110111}
    111112
    … …  
    173174      G4cout << "PhysicsList::Set " << name << " EM physics" << G4endl;
    174175
    175   } else if (name == "standard") {
     176  } else if (name == "emstandard_opt3") {
     177    emName = name;
     178    delete emPhysicsList;
     179    emPhysicsList = new G4EmStandardPhysics_option3();
     180    if (verboseLevel > 0)
     181      G4cout << "PhysicsList::Set " << name << " EM physics" << G4endl;
     182
     183  } else if (name == "emstandard_local") {
    176184    emName = name;
    177185    delete emPhysicsList;
    … …  
    179187    if (verboseLevel > 0)
    180188      G4cout << "PhysicsList::Set " << name << " EM physics" << G4endl;
    181 
    182   } else if (name == "standardIG") {
    183     emName = name;
    184     delete emPhysicsList;
    185     emPhysicsList = new PhysListEmStandardIG();
    186     if (verboseLevel > 0)
    187       G4cout << "PhysicsList::Set StandardIG EM physics" << G4endl;
    188189
    189190  } else if (name == "livermore") {
    … …  
    255256void PhysicsList::SetCuts()
    256257{
    257 
    258258  SetCutValue(cutForGamma, "gamma", "DefaultRegionForTheWorld");
    259259  SetCutValue(cutForElectron, "e-", "DefaultRegionForTheWorld");
    260260  SetCutValue(cutForPositron, "e+", "DefaultRegionForTheWorld");
    261   G4cout << "PhysicsList: world cuts are set cutG= " << cutForGamma/mm
    262          << " mm    cutE= " << cutForElectron/mm << " mm " << G4endl;
    263 
    264   if( !vertexDetectorCuts ) SetVertexCut(cutForElectron);
     261  //  G4cout << "PhysicsList: world cuts are set cutG= " << cutForGamma/mm
     262  //     << " mm    cutE= " << cutForElectron/mm << " mm " << G4endl;
     263
     264  //G4cout << " cutV= " << cutForVertexDetector
     265  //     << " cutM= " << cutForMuonDetector<<G4endl;
     266
    265267  G4Region* region = (G4RegionStore::GetInstance())->GetRegion("VertexDetector");
    266   region->SetProductionCuts(vertexDetectorCuts);
    267   G4cout << "Vertex cuts are set" << G4endl;
    268 
    269   if( !muonDetectorCuts ) SetMuonCut(cutForElectron);
     268  vertexDetectorCuts = region->GetProductionCuts();
     269  SetVertexCut(cutForVertexDetector);
     270  //  G4cout << "Vertex cuts are set" << G4endl;
     271 
    270272  region = (G4RegionStore::GetInstance())->GetRegion("MuonDetector");
    271   region->SetProductionCuts(muonDetectorCuts);
    272   G4cout << "Muon cuts are set" << G4endl;
    273 
     273  muonDetectorCuts = region->GetProductionCuts();
     274  SetMuonCut(cutForMuonDetector);
     275  //G4cout << "Muon cuts are set " <<muonRegion << " " << muonDetectorCuts << G4endl;
     276 
    274277  if (verboseLevel>0) DumpCutValuesTable();
    275278}
    … …  
    303306void PhysicsList::SetVertexCut(G4double cut)
    304307{
    305   if( !vertexDetectorCuts ) vertexDetectorCuts = new G4ProductionCuts();
    306 
    307   vertexDetectorCuts->SetProductionCut(cut, idxG4GammaCut);
    308   vertexDetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
    309   vertexDetectorCuts->SetProductionCut(cut, idxG4PositronCut);
    310 
     308  cutForVertexDetector = cut;
     309 
     310  if( vertexDetectorCuts ) {
     311    vertexDetectorCuts->SetProductionCut(cut, idxG4GammaCut);
     312    vertexDetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
     313    vertexDetectorCuts->SetProductionCut(cut, idxG4PositronCut);
     314  }
    311315}
    312316
    … …  
    315319void PhysicsList::SetMuonCut(G4double cut)
    316320{
    317   if( !muonDetectorCuts ) muonDetectorCuts = new G4ProductionCuts();
    318 
    319   muonDetectorCuts->SetProductionCut(cut, idxG4GammaCut);
    320   muonDetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
    321   muonDetectorCuts->SetProductionCut(cut, idxG4PositronCut);
    322 }
    323 
    324 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    325 
    326 void PhysicsList::SetMscStepLimit(G4bool val)
    327 {
    328   mscStepLimit = val;
    329 }
    330 
    331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    332 
     321  cutForMuonDetector   = cut;
     322
     323  if( muonDetectorCuts ) {
     324    muonDetectorCuts->SetProductionCut(cut, idxG4GammaCut);
     325    muonDetectorCuts->SetProductionCut(cut, idxG4ElectronCut);
     326    muonDetectorCuts->SetProductionCut(cut, idxG4PositronCut);
     327  }
     328}
     329
     330//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     331

  • trunk/examples/extended/electromagnetic/TestEm9/src/PhysicsListMessenger.cc

    r807 r1230  
    2525//
    2626//
    27 // $Id: PhysicsListMessenger.cc,v 1.3 2006/06/29 17:03:43 gunter Exp $
    28 // GEANT4 tag $Name: $
     27// $Id: PhysicsListMessenger.cc,v 1.4 2008/04/03 15:07:55 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
    … …  
    143143   { pPhysicsList->AddPhysicsList(newValue);}
    144144
    145   if( command == mscCmd )
    146   { pPhysicsList->SetMscStepLimit(mscCmd->GetNewBoolValue(newValue));}
     145  //  if( command == mscCmd )
     146  // { pPhysicsList->SetMscStepLimit(mscCmd->GetNewBoolValue(newValue));}
    147147}
    148148

  • trunk/examples/extended/electromagnetic/TestEm9/src/PrimaryGeneratorAction.cc

    r807 r1230  
    2626//
    2727// $Id: PrimaryGeneratorAction.cc,v 1.3 2006/06/29 17:03:46 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//

  • trunk/examples/extended/electromagnetic/TestEm9/src/StepMax.cc

    r807 r1230  
    2525//
    2626// $Id: StepMax.cc,v 1.3 2006/06/29 17:03:52 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm9/src/StepMaxMessenger.cc

    r807 r1230  
    2525//
    2626// $Id: StepMaxMessenger.cc,v 1.2 2006/06/29 17:03:54 gunter Exp $
    27 // GEANT4 tag $Name: $
     27// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2828//
    2929//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

  • trunk/examples/extended/electromagnetic/TestEm9/src/SteppingAction.cc

    r807 r1230  
    2626//
    2727// $Id: SteppingAction.cc,v 1.2 2006/06/29 17:03:57 gunter Exp $
    28 // GEANT4 tag $Name: $
     28// GEANT4 tag $Name: geant4-09-03-cand-01 $
    2929//
    3030//
Note: See TracChangeset for help on using the changeset viewer.