New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

step_length_max not propagated to Geant4 #327

Closed
petricm opened this Issue Mar 1, 2018 · 21 comments

Comments

Projects
None yet
6 participants
@petricm
Member

petricm commented Mar 1, 2018

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 petricm added the bug label Mar 1, 2018

@petricm

This comment has been minimized.

Member

petricm commented Mar 1, 2018

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.

@petricm petricm changed the title from step_length_max not propagated to G4 to step_length_max not propagated to Geant4 Mar 1, 2018

@gaede

This comment has been minimized.

Contributor

gaede commented Mar 1, 2018

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

This comment has been minimized.

Contributor

rete commented Mar 1, 2018

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

This comment has been minimized.

Contributor

gaede commented Mar 1, 2018

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

This comment has been minimized.

Contributor

MarkusFrankATcernch commented Mar 2, 2018

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

This comment has been minimized.

Contributor

gaede commented Mar 2, 2018

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()

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

This comment has been minimized.

Contributor

MarkusFrankATcernch commented Mar 2, 2018

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

This comment has been minimized.

Contributor

gaede commented Mar 2, 2018

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

This comment has been minimized.

nalipour commented Mar 2, 2018

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

This comment has been minimized.

Contributor

MarkusFrankATcernch commented Mar 2, 2018

@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

This comment has been minimized.

Contributor

shaojunlu commented Mar 5, 2018

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

This comment has been minimized.

Contributor

MarkusFrankATcernch commented Mar 6, 2018

@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

This comment has been minimized.

Contributor

MarkusFrankATcernch commented Mar 6, 2018

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

This comment has been minimized.

Contributor

MarkusFrankATcernch commented Mar 7, 2018

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

This comment has been minimized.

Contributor

gaede commented Mar 8, 2018

@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

This comment has been minimized.

Contributor

MarkusFrankATcernch commented Mar 9, 2018

@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

This comment has been minimized.

Contributor

shaojunlu commented Mar 9, 2018

@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:
steplength_head-2018-03-08
The default step_length_max=5.0mm cut and this step_length_max=1.0mm cut in one plots.
steplength_head-2018-03-08_pdg11
@MarkusFrankATcernch I would like to conclude that your implementation works fine as we expected.
Please merge it.

@MarkusFrankATcernch

This comment has been minimized.

Contributor

MarkusFrankATcernch commented Mar 12, 2018

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

This comment has been minimized.

Member

petricm commented Mar 12, 2018

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

This comment has been minimized.

Contributor

shaojunlu commented Mar 12, 2018

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

This comment has been minimized.

nalipour commented Mar 14, 2018

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment