kinetic_energy_threshold is called from src/common_cpp/Simulation/MAUSStackingAction and src/common_cpp/Simulation/DetectorConstruction.

MAUSStackingAction kills the track if at new track time if kinetic energy is less than kinetic_energy_threshold

    if (aTrack->GetKineticEnergy() < _kinetic_energy_threshold) {
        std::string ke = STLUtils::ToString(aTrack->GetKineticEnergy());
        std::string thresh = STLUtils::ToString(_kinetic_energy_threshold);
        MAUSGeant4Manager::GetInstance()->GetTracking()->SetKillReason
                (aTrack, "Kinetic Energy "+ke+" less than threshold "+thresh);
        return fKill;
    }

DetectorConstruction sets thresholds using G4UserLimits(stepMax, trackMax, timeMax, kinMin) where kinMin defaults to kinetic_energy_threshold but can be overloaded using G4KinMin on the track. Geant4 documentation reports that tracks with kinetic energy less than kinMin will be killed.

Hi Chris,

I have run simulations with G4KinMin set to both 0.001 and 0.0001 in the KL volume. The output is the same in both cases so it may be that the cut is being applied by MAUSStackingAction. I looked at pulling the value from the MICEModule for MAUSStackingAction in the same way as for src/common_cpp/Simulation/DetectorConstruction but was not able to implement this. Do you have any suggestions for how to go about setting the threshold in MAUSStackingAction on a volume-by-volume basis?

Additionally in the previous round of simulations I also lowered the thresholds production_threshold & kinetic_cutoff. Are these used by MAUS?

Hi Chris,

I have defined a new region KLregion which encompasses the KL volume in the geometry. I then set the production_cuts for this region independently of the values for the world region. The region is defined in DetectorConstruction.cc and the cuts are implemented in MAUSPhysics.cc as:

cutDoubleKL   =  _productionThresholdKL;
      // Production thresholds for detector regions
      G4Region* region;
      G4ProductionCuts* cuts;
      G4String volName;
      G4VPhysicalVolume* volume;

      regName = "KLregion";
      volName = "DummyPV";
      region = G4RegionStore::GetInstance()->GetRegion(regName);
      volume = G4PhysicalVolumeStore::GetInstance()->GetVolume(volName);
      cuts = new G4ProductionCuts;
      region->SetWorld(volume);
      cuts->SetProductionCut(0.0005*mm); // same cuts for gamma, e- and e+
      region->SetProductionCuts(cuts);

However I am not convinced that the geometry is implementing this region correctly as the output remains unchanged. Will these values be overwritten by the default values when MAUS simulates the geometry? All of the changes I have made are in my branch lp:~j-nugent-1/maus/develG4beamline

The region which I defined above is defined in void MAUSPhysicsList::SetEnergyLoss(eLoss eLossModel) and implemented if the model is energyStraggling.

void MAUSPhysicsList::SetEnergyLoss(eloss eLossModel) {
  double       cutDouble = _productionThreshold;
  double       cutDoubleKL = _productionThresholdKL;
  G4String regName;
  std::stringstream cutregion;
  std::stringstream cutStream;
  std::stringstream cutStreamRegion;
  std::string  elossActive = "activate";
  std::string  flucActive  = "true";
  switch (eLossModel) {
    case energyStraggling:
            {      elossActive = "activate";
      flucActive  = "true";
      cutDouble   =  _productionThreshold;
      cutDoubleKL   =  _productionThresholdKL;
      // Production thresholds for detector regions
      G4Region* region;
      G4ProductionCuts* cuts;
      G4String volName;
      G4VPhysicalVolume* volume;

      regName = "KLregion";
      volName = "DummyPV";
      region = G4RegionStore::GetInstance()->GetRegion(regName);
      volume = G4PhysicalVolumeStore::GetInstance()->GetVolume(volName);
      cuts = new G4ProductionCuts;
      region->SetWorld(volume);
      cuts->SetProductionCut(0.0005*mm); // same cuts for gamma, e- and e+
      region->SetProductionCuts(cuts);
//      cutStreamRegion << cutDoubleKL;
//      cutregion << regName; 
//      uiCommand.push_back("/run/setCutForRegion "+cutregion.str()+" "+cutStreamRegion.str());
      }
      break;
    case dedx:
      elossActive = "activate";
      flucActive  = "false";
      cutDouble   = 1e12*mm;
      break;
    case no_eloss:
      elossActive = "inactivate";
      flucActive  = "false";
      cutDouble   = 1e1*mm;
      break;
  }
  cutStream << cutDouble;
  cutStreamRegion << cutDoubleKL;
  cutregion << regName;
  std::vector<std::string> uiCommand;
  uiCommand.push_back("/run/setCut "+cutStream.str());
//  uiCommand.push_back("/run/setCutForRegion "+cutStreamRegion.str());
//  uiCommand.push_back("/run/setCutForRegion "+cutregion.str()+" "+cutStreamRegion.str());
  uiCommand.push_back("/process/eLoss/fluct "+flucActive);
  for (int i = 0; i < _nELossNames; i++)
    uiCommand.push_back("/process/"+elossActive+" "+_eLossNames[i]);

  for (size_t i = 0; i < uiCommand.size(); i++) {
    UIApplyCommand(uiCommand[i]);
  }
}

Humm let me look. Did you push your code anywhere?

Cheers

Everything should be on my branch on launchpad lp:~j-nugent-1/maus/develG4beamline.

Cheers

Okay, I pushed some code to

bzr+ssh://bazaar.launchpad.net/~chris-rogers/maus/1497/

That looks like it is doing something reasonable. Example/tests are in

tests/integration/test_simulation/test_physics_model_brief/test_fine_grained_production_threshold.py
tests/integration/test_simulation/test_physics_model_brief/fine_grained_production_threshold_config.py
tests/integration/test_simulation/test_physics_model_brief/FineGrainedProductionThresholdTest.dat

Should enable setting up G4Regions through the MiceModules by adding a

PropertyString Region ARegionName

and setting production threshold against these regions using fine_grained_production_threshold datacard like

fine_grained_production_threshold = {
# makes electrons, no gammas
"ARegionName":{
    "-11":0.01,
    "11":0.01,
    "22":1e12,
    "PDG PID Index":<Production_threshold>
},
}

I need to do a little more testing, documentation, merge with trunk - but it looks like it is doing something reasonable.

Thanks Chris I'll merge in the changes and try running with it.

Hi Chris,

PropertyString Region ARegionName

Does the region have to be given a position and volume or by defining in the MICE module does it automatically share the properties of the module? Also should I remove the definition of a region from DetectorConstruction.cc & MAUSPhysicsList.cc as well?

Cheers,
John

Does the region have to be given a position and volume or by defining in the MICE module does it automatically share the properties of the module?

It shares the properties of the module. E.g. look at the example, there is no special set up, I just add the Region property and then the MAUS code tells Geant4 that the logical volume defined by the MICEModule is part of the G4Region with the given name.

Also should I remove the definition of a region from DetectorConstruction.cc & MAUSPhysicsList.cc as well?

I don't know what you mean. In the branch I referred you to, the Region set up is handled in DetectorConstruction by method AddToRegion and the ProductionThreshold is put in by method SetProductionThresholdByVolume (called from SetEnergyLoss). You need this code to make it work.

Nb: the problem I am tackling - looks like volumes which are not in a region are not getting the correct default settings once I implemented the Region stuff...

Thanks for the help, I've run a few test jobs and it appears to be working. About you last comment the only other region which is currently defined is the world region does this mean that the world region is set to whatever production_threshold is used in the KL region? I would like to run some longer jobs but I'm concerned about the memory usage blowing up.

does this mean that the world region is set to whatever production_threshold is used in the KL region?

The behaviours are as follows:

I think that we want to never take the Geant4 defaults, say unless production_threshold datacard is negative.

Some more work. I added a local store of defined regions to DetectorConstruction in order to catch "Regions defined in MiceModules but not in datacards"; I added option to set a per-region default; I added the capability to use G4 name string instead of pid; I added default to G4 by setting a negative threshold. My test cases now look like:

Need to:

• automate the tests
• check I didn't break anything else
• document
• merge with trunk

Committed to bzr+ssh://bazaar.launchpad.net/~chris-rogers/maus/1497/ as r760

Unit tests and style tests all pass on my machine, try building on the test server...

Please note I am getting some test fails at integration level - e.g. nasty segmentation faults. I will have a look tomorrow... test results are visible if you click around on:

http://test.mice.rl.ac.uk/job/MAUS_rogers

I think it is now passing all tests. I also added a blob on documentation. I have had some balls up with the test server, but hopefully it should build and pass tests today sometime so that I can make a merge with the trunk.

Now merged with the trunk... I think assuming trunk tests okay this closes the issue?

By email from John

I have run a few test jobs now however for the same starting conditions changing the fine_grained_production_threshold does nothing to change the output. I have pushed the code I am using to my branch on Launchpad. Inspecting the debugging output I see that the thresholds are indeed set according to the datacard however it returns the warning message:

-------- WWWW ------- G4Exception-START -------- WWWW -------
*** G4Exception : GeomMgt1001
      issued by : G4RegionStore::GetRegion()
Region NOT found in store !
        Region KL_region NOT found in store !
        Returning NULL pointer.
*** This is just a warning message. ***
-------- WWWW -------- G4Exception-END --------- WWWW -------

I define the KL_region in the $MAUS_ROOT_DIR/src/legacy/FILES/Models/Modules/KL/KL.dat as:

Module KL
{
  Volume Box
  Dimensions 932 927.2 40.2 mm
  PropertyDouble G4StepMax 5.0 mm
  PropertyDouble G4KinMin 0.0001 MeV
  PropertyString Region KL_region    # This line defines the KL_region
  Module KL/KLFoilModule0.dat
  {
    Position 0.0 -397.5 0 mm
    Rotation 0.0 0.0 0.0 degree
  }
  Module KL/KLFoilModule1.dat
  etc. etc.

is there a problem with this definition?

I made up a geometry with two KLs in it, both with regions defined. One KL had production threshold set to 1e12, the other had production threshold set to -1. The PIDs produced by a muon track are as follows:

Pids for region KL1 [-13]
Pids for region KL2 [-13, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 22, 11]

So I believe the code is working. I have added this to the integration test, I will merge it later today if all goes well. You can see the test implementation here:

http://bazaar.launchpad.net/~chris-rogers/maus/1539/view/head:/tests/integration/test_simulation/test_physics_model_brief/test_production_threshold.py
http://bazaar.launchpad.net/~chris-rogers/maus/1539/view/head:/tests/integration/test_simulation/test_physics_model_brief/production_threshold_config.py
http://bazaar.launchpad.net/~chris-rogers/maus/1539/view/head:/tests/integration/test_simulation/test_physics_model_brief/ProductionThresholdTest.dat

John, can you confirm you are happy with this? It is all now in the trunk...

Hi Chris,

unfortunately I haven't been able to get it to work for any of the simulations that I have run. Regardless of whether I set the thresholds high or low the output is identical. I am attempting to run your tests on my build but the tests fail on an earlier test. As far as I can see my DetectorConstruction.cc and MAUSPhysics.cc are identical to the one on your branch, the region is defined in the geometry in the same way and the thresholds are set in the data card. I am not sure what else I could be missing.

Ah. I can see that e.g. the test server is running the tests and they are passing. If they fail on your machine it is quite likely you have made a mistake somewhere. Can you check out a fresh copy of the trunk and try that? You can always do

./install_build_test /PATH/TO/some_other_maus

and it will use the third party libraries from some_other_maus (e.g. geant4, root, etc). So the build time should be short. If trunk passes tests, then you know your branch is broken; if trunk fails tests then try working with a fresh install of third party stuff; if things are still failing, then we will have to think again.

I pulled down a copy of lp:maus and built with my current third-party libraries. However it is still reporting the same test failure

ok
Try saving a few standard events ... ok
Try saving a standard header; check we can only save job header once ... ok
test_OutputCppRoot: Try saving a run header ... ok
test_geometries ... ok
Check we can import libxml2 ... ok
Check we can import libxslt ... ok

======================================================================
FAIL: test_run_test_cpp_unit
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/afs/phas.gla.ac.uk/user/j/jnugent/private/PhD_Year_1/Coding/maus/maus/build/test_cpp_unit.py", line 23, in test_run_test_cpp_unit
    self.assertEqual(cpp_test, 0, 'Failed cpp tests with '+str(cpp_test))
AssertionError: Failed cpp tests with 256

----------------------------------------------------------------------
Ran 525 tests in 71.033s

FAILED (SKIP=12, failures=1)
FAIL Failed unit tests. Fatal error - aborting

I will create an entirely new installation and see what happens there but obviously it will take a few hours before I know whether it is working or not.

Cheers,
John

Ack. Can you put the full test log up? Cheers.

Thanks. I guess it's too late now, but you can ignore that test failure - it is something obscure to do with the global PID. It is already an open issue under #1492.

The trunk now has a fix for #1492
John, let me know if you still see an error with the latest trunk

Hi John -- if you pulled in lp:maus, the stable release, it won't have the fix.
The fix is in the trunk -- lp:maus/merge

John, did you ever get this working?

Cheers,
Chris

Sorry Chris, I am attending the neutrino summer school at St. Andrews this week so I haven't had a chance to try it. If I get any time over the next few days I'll rebuild MAUS and let you know the result.

This is an issue somewhere in global particle identification code, probably it is the test that is broken and not the code. I think it is safe to ignore this test failure and continue using the code.

Just to close this issue, the code I have from the merge branch builds and runs the tests without complaint and sets the threshold values correctly for each volume.

Cheers,
John