Project

General

Profile

Feature #1497

Require to set production threshold per volume

Added by Rogers, Chris over 9 years ago. Updated about 9 years ago.

Status:
Closed
Priority:
High
Assignee:
Category:
Simulation
Target version:
Start date:
25 June 2014
Due date:
% Done:

100%

Estimated time:
Workflow:
New Issue

Description

Currently Kinetic energy threshold is set per volume using G4KinMin on the volume. We should

  • Check this works okay (I am worried there may be a per track threshold in MAUSTrackingAction or something which uses the same parameter)
  • Add a similar production threshold value

Files

install.log (912 KB) install.log Nugent, John , 08 August 2014 11:01
install.log (10.5 MB) install.log Nugent, John , 09 August 2014 22:13
install.log (10.6 MB) install.log Nugent, John , 01 September 2014 13:29

Related issues

Related to MAUS - Bug #1492: Cpp unit test failure in trunkOpenRajaram, Durga23 June 2014

Actions
Related to MAUS - Bug #1548: PIDVarATestSetUp fails in trunkOpen01 September 2014

Actions
#1

Updated by Rogers, Chris over 9 years ago

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.

#2

Updated by Nugent, John about 9 years ago

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?

#3

Updated by Nugent, John about 9 years ago

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

#4

Updated by Rogers, Chris about 9 years ago

  • Target version set to Future MAUS release

I don't think this will work. The productionThreshold is done on a per-process basis in MAUSPhysicsList (but all set to the same value). See, for example void MAUSPhysicsList::SetEnergyLoss(eloss eLossModel) in src/common_cpp/Simulation/MAUSPhysicsList.cc. I don't quite see how the G4ProductionCuts as you describe interacts with this stuff...

#5

Updated by Nugent, John about 9 years ago

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]);
  }
}
#6

Updated by Rogers, Chris about 9 years ago

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

Cheers

#7

Updated by Nugent, John about 9 years ago

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

Cheers

#8

Updated by Rogers, Chris about 9 years ago

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.

#9

Updated by Nugent, John about 9 years ago

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

#10

Updated by Nugent, John about 9 years ago

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

#11

Updated by Rogers, Chris about 9 years ago

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...

#12

Updated by Nugent, John about 9 years ago

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.

#13

Updated by Rogers, Chris about 9 years ago

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

The behaviours are as follows:

Region defined in MiceModules? Thresholds defined in datacards Result
Yes Yes Production thresholds set per particle as described above
Yes Yes but not for all particles Not defined particles acquire some Geant4 default
Yes No Acquires some Geant4 default
No No Acquires datacards production_threshold

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

#14

Updated by Rogers, Chris about 9 years ago

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:

Region defined in MiceModules? Thresholds defined in datacards Result
Yes Yes using pdg pid Production thresholds set per particle as described above
Yes Yes using G4 particle name Production thresholds set per particle as described above
Yes Yes but not for all particles Not defined particles acquire productionThreshold default
Yes Yes but "default" also specified Not defined particles acquire "default"
Yes Yes but negative Not defined particles acquire Geant4 default
Yes No Acquires production threshold
No No Acquires datacards production_threshold
No Yes Ignored

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

#15

Updated by Rogers, Chris about 9 years ago

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

#16

Updated by Rogers, Chris about 9 years ago

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

#17

Updated by Rogers, Chris about 9 years ago

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.

#18

Updated by Rogers, Chris about 9 years ago

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

#19

Updated by Rogers, Chris about 9 years ago

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?
#20

Updated by Rogers, Chris about 9 years ago

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

#21

Updated by Rogers, Chris about 9 years ago

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

#22

Updated by Nugent, John about 9 years ago

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.

#23

Updated by Rogers, Chris about 9 years ago

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.

#24

Updated by Nugent, John about 9 years ago

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

#25

Updated by Rogers, Chris about 9 years ago

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

#26

Updated by Nugent, John about 9 years ago

Here is the log.

#27

Updated by Rogers, Chris about 9 years ago

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.

#28

Updated by Rajaram, Durga about 9 years ago

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

#29

Updated by Nugent, John about 9 years ago

Hi Durga,

I just pulled down lp:maus and built on my glasgow machine. It failed on the same test as before, I have attached the entire install log.

Cheers,
John

#30

Updated by Rajaram, Durga about 9 years ago

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

#31

Updated by Rogers, Chris about 9 years ago

John, did you ever get this working?

Cheers,
Chris

#32

Updated by Nugent, John about 9 years ago

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.

#33

Updated by Nugent, John about 9 years ago

Apologies for resurrecting this issue, I was not able to make any time for this during the school. I have built lp:maus/merge on my linux machine and the install fails with the attached log. If either Chris or Durga have any time to look into this it would be a great help.

Cheers,
John

#34

Updated by Rogers, Chris about 9 years ago

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.

#35

Updated by Nugent, John about 9 years ago

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

#36

Updated by Rogers, Chris about 9 years ago

  • Status changed from Open to Closed
  • % Done changed from 0 to 100

Kewl, issue is closed.

Also available in: Atom PDF