step_length_max not propagated to Geant4

[petricm]

The step_length_max that is set in <limits> is not propagated to G4.
Individually corroborated by @shaojunlu and @nalipour that no step limit is observed when looking at step length.

[petricm]

From the Geant4 manual

A G4UserLimits object must be instantiated for the duration of whatever logical volume or region to which it is assigned. It is the responsibility of the user's code to delete the object after the assigned volume(s)/region(s) have been deleted.

In addition to instantiating G4UserLimits and setting it to logical volume or region, the user has to assign the following process(es) to particle types he/she wants to affect. If none of these processes is assigned, that kind of particle is not affected by G4UserLimits.

It seems that DD4hep does not do the second step.

[gaede]

Here is what I believe is the relevant code used in Mokka:

void PhysicsListUserLimits::Enable(void)
{
  // cf. Geant 4 HyperNews, Forum "Physics List", Message 129
  // http://geant4-hn.slac.stanford.edu:5090/HyperNews/public/get/phys-list/129.html

  G4UserSpecialCuts *specialCuts = new G4UserSpecialCuts;
  G4StepLimiter     *stepLimiter = new G4StepLimiter;
  TPCStepLimiterLowPt     *tpcStepLimiterLowPt = new TPCStepLimiterLowPt;

  G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
  G4ParticleTable::G4PTblDicIterator *particleIterator = particleTable->GetIterator();
  // make sure you have called "G4RunManager::Initialize()" before

  particleIterator->reset();
  while ((*particleIterator)()) {
  // iterate through all known particles

    G4ParticleDefinition *particleDefinition = particleIterator->value();
    G4ProcessManager *processManager = particleDefinition->GetProcessManager();

    if (processManager && !particleDefinition->IsShortLived() && particleDefinition->GetPDGCharge() != 0) {
    // the process manager should exist, but we don't need to limit short-lived particles or neutrals

      processManager->AddDiscreteProcess(stepLimiter);
      processManager->AddDiscreteProcess(specialCuts);
      // these transportation-related processes are always applicable

      // if TPCStepLimiterLowPt set in steering file
      if( Control::TPCLowPtStepLimit ) {
        processManager->AddDiscreteProcess(tpcStepLimiterLowPt);
      }
    }
  }
}

[rete]

Also, in this scheme, the limits apply to all particles entering in the volume where the limits are defined. This is how we currently define the limits in one of our ILD detector :

 <limits>
    <limitset name="cal_limits">
      <limit name="step_length_max" particles="*" value="cal_steplimit_val" unit="cal_steplimit_unit" />
    </limitset>
    <limitset name="TPC_limits">
      <limit name="step_length_max" particles="*" value="tpc_steplimit_val" unit="tpc_steplimit_unit" />
    </limitset>
    <limitset name="Tracker_limits">
      <limit name="step_length_max" particles="*" value="tracker_steplimit_val" unit="tracker_steplimit_unit" />
    </limitset>
</limits>

I don't understand how the limits could apply to a specific set of particles in a given volume but not to other particles.
For example, if I do this :

 <limits>
    <limitset name="cal_limits">
      <limit name="step_length_max" particles="gamma" value="cal_steplimit_val" unit="cal_steplimit_unit" />
    </limitset>
    <limitset name="TPC_limits">
      <limit name="step_length_max" particles="e-" value="tpc_steplimit_val" unit="tpc_steplimit_unit" />
    </limitset>
</limits>

the step limiter process will be assigned to the electron particle but by design Geant4 will apply the step limits to the electron also where the cal_limits are applied because the process is globally assigned to the particle. Am I clear ? ...

[gaede]

From what I can see, the particles in the limitset are never really used.
Maybe we can drop these global limits completely - as they should be defined in the volumes/regions of the detector ?
I think @MarkusFrankATcernch needs to have a look and say where and when the above code needs to be applied.

[MarkusFrankATcernch]

Frank,
I think your code snippet is somehow incomplete.
There are no values propagated to the G4StepLimiter nor to the G4UserSpecialCuts.
There must be another fragment, where these things are actually set.
Could you have a second look please ?

[gaede]

Markus, the way this works is that one has to add the StepLimiter to all charged long lived particles. The actual limits are then taken from the volume or the region (so probably nothing to be done there in DD4hep).
I just found the following example code in examples/basic/B5/exampleB5.cc :

  auto physicsList = new FTFP_BERT;
  physicsList->RegisterPhysics(new G4StepLimiterPhysics());
  runManager->SetUserInitialization(physicsList);

The G4StepLimiterPhysics includes the same code as above for Mokka in ConstructProcess()

• see here:
  http://www.apc.univ-paris7.fr/~franco/g4doxy4.10/html/class_g4_step_limiter_physics.html

An extra complication is the TPCStepLimiterLowPt we used to have optionally in Mokka. This we eventually also need in DD4hep - I can look into this, once we have thew basic user step limits enabled and working.

[MarkusFrankATcernch]

Yes. I found that as well.
It is even easier (and already implemented from python):

Now build the physics list:

  phys = simple.setupPhysics('QGSP_BERT')
  ph = DDG4.PhysicsList(kernel,'Geant4PhysicsList/Myphysics')
  #....
  #Add this:
  ph.addPhysicsConstructor('G4StepLimiterPhysics')
  #....
  phys.add(ph)

Would it be possible, that @shaojunlu tries this out?

[gaede]

Ok, this works if we add it to DDSim - is this already the fix or should this also be added to DDG4 as default behavior ?

[nalipour]

Hi,
For FCC, DDSim is not used. So it would be useful to have it in DDG4 otherwise I can not set this parameter.

Thank you,
Nilou

[MarkusFrankATcernch]

@nalipour : You will have to instruct Gaussino to do this for you.
At this level DD4hep cannot help you, because DD4hep does neither have a handle to your run manager nor to your physics list. Hence somewhere you have to apply Frank's solution:

auto physicsList = new FTFP_BERT;
physicsList->RegisterPhysics(new G4StepLimiterPhysics());
runManager->SetUserInitialization(physicsList);

[shaojunlu]

It works if we add it to DDSim in this way.
- is this already the fix or should this also be added to DDG4 as default behaviour ?

[MarkusFrankATcernch]

@shaojunlu I do not know ...
This should become only the default behaviour if this is ALWAYS wanted.
But why is then not in Geant4?

Otherwise the behavior should be implemented in ddsim.

[MarkusFrankATcernch]

There is still the issue mentioned by @rete from above.
I must say I am not competent enough to answer this. If this behavior can be achieved in Geant4,
then I shall implement it, if someone can give me the corresponding recipe. Otherwise, I think, we only document that the "particles" attribute in the Limit object is unused.

[MarkusFrankATcernch]

I have committed a dd4hep version in
https://github.com/MarkusFrankATcernch/DD4hep
which should set particle limits exactly the way the limit-sets say:

• multiple limits depending on the particle type.
• if the particle type is not explicitly found, defaults kick in:
• Limit.particles == "*" and if not present, the G4UserLimits default value

Could someone please test this?

• obviously the step limited physics needs to be added:
  ph.addPhysicsConstructor('G4StepLimiterPhysics')

It would be good if we could have an explicit test checking such a behavior.

[gaede]

@shaojunlu I am currently installing a HEAD release with this version of DD4hep at:

/afs/desy.de/project/ilcsoft/sw/x86_64_gcc49_sl6/HEAD-2018-03-08

Once this is done, could you please run a short test with your tools that check the step size for a setting that specifies different additional limits to say electrons and muons ?

@MarkusFrankATcernch could you provide example xml snippets that set limits for e+- and mu+- only ?

[MarkusFrankATcernch]

@gaede: For me such a limitset definition works (or does not give errors):

  <limits>
    <limitset name="minitel_limits">
      <limit name="step_length_max" particles="e[+-]" value="1.0" unit="mm" />
      <limit name="step_length_max" particles="mu[+-]" value="3.0" unit="mm" />
      <limit name="step_length_max" particles="*" value="5.0" unit="mm" />
    </limitset>
  </limits>

[shaojunlu]

@gaede Your three step_length_max have been used in our ILD_l5_v02 simulation.
This is the output from SimCalorimeterHits MCParticles contribution from Hybird HCAL:

The default step_length_max=5.0mm cut and this step_length_max=1.0mm cut in one plots.

@MarkusFrankATcernch I would like to conclude that your implementation works fine as we expected.
Please merge it.

[MarkusFrankATcernch]

Hi @shaojunlu ,

1. Thanks a lot for doing the tests!
2. Could you please send me the code to produce these histograms?
   I will try to somehow build a standalone example and try to incorporate it in the standard tests.
3. Do we need something like this also for G4Region objects ?

Cheers,

Markus

[petricm]

As for your question for Regions, we need the same behavior as Geant

User limits assigned to logical volume do not propagate to daughter volumes, while User limits assigned to region propagate to daughter volumes unless daughters belong to another region

[shaojunlu]

Dear @MarkusFrankATcernch

Recently, ILCsoft LCIO has added "SimCalorimeterHit::getLengthCont() (step length)".
Please have a look at "SimCalorimeterHit" implementation.
https://github.com/iLCSoft/LCIO/tree/master/src/cpp/src/IMPL
https://github.com/iLCSoft/LCIO/blob/master/src/cpp/src/IMPL/SimCalorimeterHitImpl.cc#L117

After the simulation as README.md documented in ILDConfig/StandardConfig/production.
https://github.com/iLCSoft/ILDConfig/blob/master/StandardConfig/production/README.md#2-run-the-lcgeoddsim-simulation-the-3-ttbar-example

SimCalorimeterHit collections from the calorimeter detectors will include the "step length" information.

    const LCCollection* col =  evt->getCollection( _simCaloHitCollectionName ) ; 
    if( col ){
      streamlog_out(DEBUG4) << " We examine collection " << _simCaloHitCollectionNames<< " with " << col->getNumberOfElements() << " hits " << std::endl;
	
      for( int j = 0, jN = col->getNumberOfElements() ; j<jN ; ++j ) {
	SimCalorimeterHit* simHit = (SimCalorimeterHit*) col->getElementAt( j ) ;

	for( int k = 0, kN = simHit->getNMCContributions() ; k<kN ; ++k ) {
	  int NMCC_PDG = simHit->getPDGCont(k) ;
	  float NMCC_stepLength = simHit->getLengthCont(k) ;

	  // NMCC_PDG is this contribution PDG code.
	  // NMCC_stepLength is this contribution step length.
	}

      }

    }

   // evt is the LCEvent which is defined in LCIO
   // _simCaloHitCollectionName is one calorimeter collection name to be checked within this event.

[nalipour]

Hi,

Thank you for fixing the issue. Is it possible to make a tag so that we can use it in FCCSW?

Thanks,
Best regards,
Nilou