search

AIDASoft / DD4hep

Detector Description Toolkit for High Energy Physics
http://dd4hep.cern.ch
GNU Lesser General Public License v3.0
53 stars 100 forks source link

Program never stops #1027

Closed chekanov closed 1 year ago

chekanov commented 1 year ago

Hi,

I'm not sure is this is a bug or not, but I'm posting this problem here. I'm trying to design a simple ECAL with PbWO4 crystals, using the CLIC calorimeter geometry. If I just use PbWO4 without optical properties, everything works. When optical properties are added to PbWO4 , the program runs forever.

Here is GIT to reproduce this problem: https://github.com/chekanov/SimpleECAL

Note that this example is similar to SingleDualCrystal (used by Sarah Eno), which works fine for single crystals.

This infinite loop problem can be fixed by commenting out

scint = PhysicsList(kernel, 'Geant4ScintillationPhysics/ScintillationPhys')

but keeping Geant4CerenkovPhysics/CerenkovPhys and Geant4OpticalPhotonPhysics/OpticalGammaPhys in the steering card SCEPCALsteering.py

MarkusFrankATcernch commented 1 year ago

I have built the example. How do I start the simulation ?

chekanov commented 1 year ago

You can simply run a script from "compact":

> A_RUN

that makes a ROOT file or interactive method:

> A_INTERACTIVE

Then run:

/control/execute vis.mac
/run/beamOn 1
MarkusFrankATcernch commented 1 year ago

Is the A_INTERACTIVE in the repository supposed to reproduce the problem? With the current version I was happily running:

DRcalo deposit 0 position (43,52,0) string /ECalBarrel_0/stave_6/layer4_3/slice3_2 volume id 14636883489332340
DRcalo deposit 0 position (28,47,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 13229444181269588
DRcalo deposit 0 position (28,48,0) string /ECalBarrel_0/stave_6/layer4_3/slice3_2 volume id 13510919157980276
DRcalo deposit 0 position (26,44,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 12385010661203028
DRcalo deposit 0.069689 position (25,42,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 11822056412814420
DRcalo deposit 1.00183 position (25,41,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 11540581436103764
DRcalo deposit 1.91461 position (25,43,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 12103531389525076
DRcalo deposit 0 position (28,45,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 12666494227848276
DRcalo deposit 1.71739 position (29,44,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 12385023546104916
DRcalo deposit 0 position (29,45,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 12666498522815572
DRcalo deposit 0 position (30,45,0) string /ECalBarrel_0/stave_6/layer4_3/slice2_1 volume id 12666502817782868
GenerationInit   WARN  +++ Finished run 1 after 5 events (6 events in total)
0 events have been kept for refreshing and/or reviewing.
  "/vis/reviewKeptEvents" to review them one by one.
  "/vis/enable", then "/vis/viewer/flush" or "/vis/viewer/rebuild" to see them accumulated.
Idle> exit
Graphics systems deleted.
Visualization Manager deleting...

Geant4Kernel     INFO  ++ Terminate Geant4 and delete associated actions.
RootOutput       INFO  +++ Closing ROOT ourpur file testSCEPCAL.root
DDSim            INFO DDSim            INFO  Total Time:   18.87 s (User), 1.10 s (System)
[Thread 0x7fd9e78e2700 (LWP 228062) exited]
[Inferior 1 (process 228004) exited normally]
(gdb) q

... and thanks for letting me spot the typo: Closing ROOT ourpur file.... If such typos are spotted everyone should somehow bring it up.

chekanov commented 1 year ago

Yes, I posted a working version with PbWO4 defined without the optical properties. Please change PbWO4 to E_PbWO4 on the line 42 of compact/ECalBarrel_DualCrystal.xml and run this example again. Then this iterative example will run forever. I suggest to have a second terminal to kill the program otherwise you may need a hard reboot. Note that we tested E_PbWO4 using a single crystal geometry and such an example works. The problem appears when we make a full ECAL with the crystals organized as PolyhedraRegular.

MarkusFrankATcernch commented 1 year ago

Do you then get output like the following forever:

DRcalo deposit 2.60848e-06 position (13,168,0) string /ECalBarrel_0/stave_8/layer1_0/slice2_1 volume id 47287851926168660
DRcalo deposit 2.60848e-06 position (40,114,0) string /ECalBarrel_0/stave_8/layer1_0/slice2_1 volume id 32088319147910228
DRcalo deposit 2.60848e-06 position (34,102,0) string /ECalBarrel_0/stave_8/layer1_0/slice1_0 volume id 28710593657578548
DRcalo deposit 2.72093e-06 position (-10,-55,0) string /ECalBarrel_0/stave_8/layer1_0/slice2_1 volume id -15199691687844780
DRcalo deposit 2.63858e-06 position (35,138,0) string /ECalBarrel_0/stave_8/layer1_0/slice2_1 volume id 38843697114129492
DRcalo deposit 2.63858e-06 position (14,175,0) string /ECalBarrel_0/stave_8/layer1_0/slice3_2 volume id 49258181058110580
DRcalo deposit 2.65e-06 position (-39,51,0) string /ECalBarrel_0/stave_8/layer1_0/slice2_1 volume id 14636531289433172
DRcalo deposit 2.71359e-06 position (0,-36,0) string /ECalBarrel_0/stave_7/layer1_0/slice2_1 volume id -10133099157381036
DRcalo deposit 2.71359e-06 position (41,-49,0) string /ECalBarrel_0/stave_8/layer2_1/slice1_0 volume id -13792097756765132
DRcalo deposit 2.65877e-06 position (37,-3,0) string /ECalBarrel_0/stave_9/layer1_0/slice3_2 volume id -844266012137356
DRcalo deposit 2.62378e-06 position (-2,105,0) string /ECalBarrel_0/stave_8/layer1_0/slice2_1 volume id 29836338945598548
DRcalo deposit 2.83043e-06 position (5,-103,0) string /ECalBarrel_0/stave_8/layer1_0/slice3_2 volume id -28991901122157452
DRcalo deposit 2.44291e-06 position (26,61,0) string /ECalBarrel_0/stave_8/layer1_0/slice2_1 volume id 17170085252703316
DRcalo deposit 2.44291e-06 position (17,127,0) string /ECalBarrel_0/stave_8/layer1_0/slice1_0 volume id 35747395060900916
DRcalo deposit 2.72603e-06 position (1,-32,0) string /ECalBarrel_0/stave_8/layer1_0/slice1_0 volume id -9007194955570124
^C
Thread 1 "python3.9" received signal SIGINT, Interrupt.
0x00007f79aa723dc1 in G4String::operator== (this=0xd10b8a0, str=...)
    at /cvmfs/sft.cern.ch/lcg/views/LCG_101/x86_64-centos7-gcc11-dbg/include/Geant4/G4String.icc:119
119 }

For me this kind of output continues forever.... Is this what you mean ?

chekanov commented 1 year ago

Yes, exactly. It prints forever without stop.

MarkusFrankATcernch commented 1 year ago

Hmmm. This does not look immediately like a DD4hep problem, since the Geant4 kernel keeps on tracking without ever stopping. If I remove in the material description:

  <material name="E_PbWO4">
    <D value="8.28" unit="g/cm3" />
    <fraction n="0.45532661" ref="Pb"/>
    <fraction n="0.40403397" ref="W"/>
    <fraction n="0.14063942" ref="O"/>
    <property name="RINDEX"        ref="RINDEX__PbWO4"/>
      <property name="ABSLENGTH" ref="AbsLen_PS"/>
<!--
      <property name="FASTCOMPONENT" ref="scintFast_PS"/>
-->
      <constant name="SCINTILLATIONYIELD" value="13.9/keV"/>
      <constant name="FASTTIMECONSTANT" value="2.8*ns"/>
      <constant name="RESOLUTIONSCALE" value="1."/>
  </material>

the FASTCOMPONENT, then the simulation terminates, but creates many more and much smaller deposits than without any property. I think what is needed here is a better understanding what these properties really mean for Geant4. DD4hep is here only the mediator. I checked and according to my understanding are values are passed correctly to Geant4.

I don't know if this entry here sheds any light on the topic: https://geant4-forum.web.cern.ch/t/optical-properties-of-the-scintillator-material/3611

chekanov commented 1 year ago

It it possible the the light is reflected on the wall of PolyhedraRegular geometry and travels around in the phi? Maybe there is a way to insulate sides of the crystals ? I'm not sure how this can be done since I can only see how to insulate 4 crystals in the depth, but not on the sides.. A similar example for a module made from several such crystals works fine (example by Sarah Eno https://github.com/saraheno/DualTestBeam). This problem happens when we arrange the crystals using PolyhedraRegular , i.e. ECAL barrel, circular in Phi,

MarkusFrankATcernch commented 1 year ago

Hmm. Then the problem should disappear if only one stave would be constructed. This would eliminate this bit.

chekanov commented 1 year ago

Thanks, I can check this idea. Do you have any suggestion of how to modify DD4HEP part https://github.com/chekanov/SimpleECAL/blob/main/compact/ECalBarrel_DualCrystal.xml to make one stave (or, alternatively, 1/2 of the ECAL in phi)?

MarkusFrankATcernch commented 1 year ago

You will have to add a line if ( i == <X> ) break; at line 156 of EcalBarrel_geo.cpp, where <X> is the number of staves you actually want.

changing "nsides" in xml alone will not do the job.

saraheno commented 1 year ago

So Sergei: my attempt to insulate the sides of the crystals is being discussed in this thread: #1014 dd4hep does have the possibility of putting optical surfaces between optical elements, but somehow it is not quite working. I am disucssing with Andre and with Hans Wenzel if this is a problem with dd4hep or geant4. @hanswenzel

hanswenzel commented 1 year ago

Hi Sarah et al.

the provided gdml file is not valid as far as Geant4 (I was using the latest version of Geant4) is concerned. Reading it in from CaTS I get:

------- EEEE ------- G4Exception-START -------- EEEE ------- G4Exception : ReadError issued by : G4GDMLReadDefine::getMatrix() Matrix 'SCINTILLATIONYIELD_0x2dfbf80' was not found! Fatal Exception core dump Track information is not available at this moment Step information is not available at this moment

-------- EEEE -------- G4Exception-END --------- EEEE -------

that makes sense since SCINTILLATIONYIELD is defined as a constant

while the property requires a matrix:

  <property name="SCINTILLATIONYIELD" ref="SCINTILLATIONYIELD_0x2dfbf80"/>

What is the Geant4 version used by dd4hep?

Hans

--

Dr. Hans Wenzel

Physics and Detector Simulation

Fermi National Accelerator Laboratory

P.O. Box 500, MS 223

Batavia, Illinois 60510

USA

630 840 6034 office

331 250 9342 mobile

home.fnal.gov/~wenzelhttp://home.fnal.gov/~wenzel

@.**@.>

From: Sarah Eno @.> Sent: Tuesday, December 6, 2022 10:54 AM To: AIDASoft/DD4hep @.> Cc: Hans-Joachim Wenzel @.>; Mention @.> Subject: Re: [AIDASoft/DD4hep] Program never stops (Issue #1027)

So Sergei: my attempt to insulate the sides of the crystals is being discussed in this thread: #1014https://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_AIDASoft_DD4hep_issues_1014&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=bIL908M-35QimpL2FpB7Bw&m=k5-tvqiIsvJTn_DIcyCC0IVVP9GCHli3M96TLdNiq_R0U0qw3ZikNZ_KX481q8aF&s=SWajEp7Ns7wkHjleDgDA1i6qLrNwjBoS5Hw7p8GYt9Q&e= dd4hep does have the possibility of putting optical surfaces between optical elements, but somehow it is not quite working. I am disucssing with Andre and with Hans Wenzel if this is a problem with dd4hep or geant4. @hanswenzelhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_hanswenzel&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=bIL908M-35QimpL2FpB7Bw&m=k5-tvqiIsvJTn_DIcyCC0IVVP9GCHli3M96TLdNiq_R0U0qw3ZikNZ_KX481q8aF&s=1LlEBVef89yWN84KROKUFeanMswNwRnamjap3V9pOpo&e=

— Reply to this email directly, view it on GitHubhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_AIDASoft_DD4hep_issues_1027-23issuecomment-2D1339673091&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=bIL908M-35QimpL2FpB7Bw&m=k5-tvqiIsvJTn_DIcyCC0IVVP9GCHli3M96TLdNiq_R0U0qw3ZikNZ_KX481q8aF&s=OCnrbB7qn9oFI_4hbATkdd1YY-d8Sd-ibL7Re-oF2BE&e=, or unsubscribehttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AB4O4OZT3ZXUTGFH2VLUWELWL5VTZANCNFSM6AAAAAASQAEGLM&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=bIL908M-35QimpL2FpB7Bw&m=k5-tvqiIsvJTn_DIcyCC0IVVP9GCHli3M96TLdNiq_R0U0qw3ZikNZ_KX481q8aF&s=RVMPr2MWBmxK2XajSB0bQdqag_TuWzLw2VqnzfXC4UY&e=. You are receiving this because you were mentioned.Message ID: @.***>

saraheno commented 1 year ago

@hanswenzel

I initialize my area doing

source /cvmfs/sft.cern.ch/lcg/views/LCG_101/x86_64-centos7-gcc11-opt/setup.sh

I get this spew when I run

Geant4 version Name: geant4-10-07-patch-02 [MT] (11-June-2021) Copyright : Geant4 Collaboration References : NIM A 506 (2003), 250-303 : IEEE-TNS 53 (2006), 270-278 : NIM A 835 (2016), 186-225 WWW : http://geant4.org/

hanswenzel commented 1 year ago

Hi Sarah,

that's definitely recent enough. So there is definitely something wrong in the compact to gdml conversion. It's not only the SCINTILLATIONYIELD but other properties like e.g. the FASTTIMECONSTANT that should be defined as matrices not constants.

Another problem in the gdml file is:

<bordersurface name="/world/DRCrystal#HallCrys" surfaceproperty="/world/DRCrystal#mirrorSurface">
  <physvolref ref="tower_0" ref="DRCrystal_0"/>
</bordersurface>

a bordersurface needs 2 physical volumes (compared to a skin surface that only needs 1 logical volume) e.g:

So I don't know dd4hep but the gdml file that is produced and should represent the Geant4 geometry definitely is not valid. I modified the gdml file and attached it to this mail as well as a screenshot of CaTS reading in this file. I will shoot some optical photons into this geometry to see if the surface works properly but I think Geant4 is just fine. CaTS might be a better solution for single crystals or test beam set ups 🙂 .

Hans

--

Dr. Hans Wenzel

Physics and Detector Simulation

Fermi National Accelerator Laboratory

P.O. Box 500, MS 223

Batavia, Illinois 60510

USA

630 840 6034 office

331 250 9342 mobile

home.fnal.gov/~wenzelhttp://home.fnal.gov/~wenzel

@.**@.>

From: Sarah Eno @.> Sent: Tuesday, December 6, 2022 11:52 AM To: AIDASoft/DD4hep @.> Cc: Hans-Joachim Wenzel @.>; Mention @.> Subject: Re: [AIDASoft/DD4hep] Program never stops (Issue #1027)

@hanswenzelhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_hanswenzel&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=bIL908M-35QimpL2FpB7Bw&m=R_0yitd5r1Ai3FMJjxKM08AUW-ine1zoWF0HX4K6SKimR1-2ghkAYVx35rVLVQY8&s=ZBMHAHDO1nnceeM0N3jPSd78Ix7RfA0lYFPXQz1P1VE&e=

I initialize my area doing

source /cvmfs/sft.cern.ch/lcg/views/LCG_101/x86_64-centos7-gcc11-opt/setup.sh

I get this spew when I run

Geant4 version Name: geant4-10-07-patch-02 [MT] (11-June-2021) Copyright : Geant4 Collaboration References : NIM A 506 (2003), 250-303 : IEEE-TNS 53 (2006), 270-278 : NIM A 835 (2016), 186-225 WWW : http://geant4.org/https://urldefense.proofpoint.com/v2/url?u=http-3A__geant4.org_&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=bIL908M-35QimpL2FpB7Bw&m=R_0yitd5r1Ai3FMJjxKM08AUW-ine1zoWF0HX4K6SKimR1-2ghkAYVx35rVLVQY8&s=TvNllodRmuQnZKgQriSJeukVZEz9VMGZBt_CjbIiVgg&e=

— Reply to this email directly, view it on GitHubhttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_AIDASoft_DD4hep_issues_1027-23issuecomment-2D1339753769&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=bIL908M-35QimpL2FpB7Bw&m=R_0yitd5r1Ai3FMJjxKM08AUW-ine1zoWF0HX4K6SKimR1-2ghkAYVx35rVLVQY8&s=F-guYGovHSMPYELeKNrUG8f7kLVV-wH0NFJw2djnRGI&e=, or unsubscribehttps://urldefense.proofpoint.com/v2/url?u=https-3A__github.com_notifications_unsubscribe-2Dauth_AB4O4OYEZJY2LFNFNA5KOU3WL54OXANCNFSM6AAAAAASQAEGLM&d=DwMCaQ&c=gRgGjJ3BkIsb5y6s49QqsA&r=bIL908M-35QimpL2FpB7Bw&m=R_0yitd5r1Ai3FMJjxKM08AUW-ine1zoWF0HX4K6SKimR1-2ghkAYVx35rVLVQY8&s=VnIKisRP61YNlJypVlUapqdZsr0R6sfJt-Q7EzjqBns&e=. You are receiving this because you were mentioned.Message ID: @.***>

MarkusFrankATcernch commented 1 year ago

In Geant4 geant4-11.00.b01 according to geant4-11.00.b01/source/materials/include/G4MaterialPropertiesIndex.hh

kSCINTILLATIONYIELD is part of the enum G4MaterialConstPropertyIndex. Does this mean it can still be a matrix ?

saraheno commented 1 year ago

@hanswenzel yes, we have a team at Virginia working on CaTS. But eventually we need to build a full detector that will work with the CLIC/ILC/CEPC software, so it has to be in dd4hep. We cannot build an entire detector from scratch. We want to use other's trackers etc. So I'm doing the single bars etc to understand better the code, working with Sergei to eventually get a full detector in dd4hep. while the Virginia team uses CaTS for our cosmic ray tests that are starting now. (anyway, right now dd4hep is too slow... I hope somehow we can get your accelerator code into it or we may only be able to put in parameterizations of the optics, not use the GEANT4 optics)

@MarkusFrankATcernch @andresailer The place where I attach the surface is at https://github.com/saraheno/SingleDualCrystal/blob/6589de4b77b686399b7236cf3d33de8f708833f1/src/DRCrystal_geo.cpp#L213

i think I am attaching it between pv and env_phv. is that not right?

MarkusFrankATcernch commented 1 year ago

@saraheno If the gdml is not correct: how did you produce it?

The surface stuff in ROOT is rather recent (and the surface gdml writing probably not extensively debugged), but the command from the Geant4 prompt invokes the Geant4 GDML writer. This could be an explanation which the GDML is inconsistent.

saraheno commented 1 year ago

I used this command: geoConverter -compact2gdml -input DRSingleCrystal.xml -output DRSingleCrystal.gdml

hanswenzel commented 1 year ago

Hi I just checked by dumping the geometry that I constructed with my valid gdml file directly from geant4 memory to a new gdml file and the produced gdml file just looks fine. The surfaces are defined as skinsurfaces and the relevant properties are defined as matrices. So I don't think the problem is Geant4

Hans

saraheno commented 1 year ago

@hanswenzel when you run muons through this, making cherenkov photons, do any escape through the surface? or are they all trapped?

hanswenzel commented 1 year ago

Hi Sarah

Haven't tried that yet . I will try and let you know after lunch. I just shot single optical photons. But after looking at the gdml files one first needs to actually fix the detector description. I have put 2 files into github: https://github.com/hanswenzel/CaTS/tree/main/gdml

one is DRSingleCrystal.gdml the other is sarah.gdml. The first one is the corrected gdml file that Sarah sent me the second is the dump of the geometry from Geant4 into a gdml file. Both are valid and can be used by Geant4 as input.

saraheno commented 1 year ago

@hanswenzel I understood part of the problem was that the surface in the gdml was only attached to one volume. to fix it, you must have attached it to a second volume? to which volume did you attach it?

hanswenzel commented 1 year ago

Hi there were actually several problems with the geometry: 1) constants are used when matrices are required, 2) There are 2 types of surfaces a border surface between two physical volumes (e.g. between the crystal and the photodetector) and a skin surface which defines the surface (all sides) of a "logical" volume. The border surface requires 2 physical volumes as input the skin surface one logical volume. See the gdml files that I put into the CaTS github. You will need both types (crystal to surrounding air (skin) and the end of the crystal surface needs to be interfaced with a photo detector surface(border) . It looks to me that the dd4hep implementation needs some work. The produced gdml file is not a valid gdml file. The 2 gdml files in github show the input gdml and the direct dump from Geant4 memory. I had problems in the past were the root gdml parser was slightly different from the Geant4 one creating all kinds of problems.

saraheno commented 1 year ago

@hanswenzel in your previous post, could you defined better "input gdml" and "direct dump from GEANT4 memory"? I'm not sure exactly what they are.

hanswenzel commented 1 year ago

Hi Sarah, the file https://github.com/hanswenzel/CaTS/blob/main/gdml/DRSingleCrystal.gdml is basically the file you provided and that I modified to actually make it a valid gdml file. The file: https://github.com/hanswenzel/CaTS/blob/main/gdml/sarah.gdml is the dump of the Geant4 memory to gdml which is equivalent to the input. I just wanted to make sure that nothing got corrupted between reading and writing and the Geant4 parser and writer are fine. On the other hand the original file produced by the geomconverter didn't produce a valid gdml file and the output is not a representation of a valid Geant4 geometry

Hans .

saraheno commented 1 year ago

Hi Hans, I think the dd4hep team is still waiting for you to reply to https://github.com/AIDASoft/DD4hep/issues/1027#issuecomment-1339810988

On Tue, Dec 6, 2022 at 5:50 PM Hans Wenzel @.***> wrote:

Hi Sarah, the file https://github.com/hanswenzel/CaTS/blob/main/gdml/DRSingleCrystal.gdml is basically the file you provided and that I modified to actually make it a valid gdml file. The file: https://github.com/hanswenzel/CaTS/blob/main/gdml/sarah.gdml is the dump of the Geant4 memory to gdml which is equivalent to the input. I just wanted to make sure that nothing got corrupted between reading and writing and the Geant4 parser and writer are fine. On the other hand the original file produced by the geomconverter didn't produce a valid gdml file and the output is not a representation of a valid Geant4 geometry

Hans .

— Reply to this email directly, view it on GitHub https://github.com/AIDASoft/DD4hep/issues/1027#issuecomment-1340116693, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABHZIOFKXCVXRSNBGIBQLCLWL67LNANCNFSM6AAAAAASQAEGLM . You are receiving this because you were mentioned.Message ID: @.***>

MarkusFrankATcernch commented 1 year ago

@saraheno

I need to know if /ddg4/ConstructGeometry/writeGDML <file-name> produces a correct GDML to know if the problem is in the conversion to Geant4 or in the ROOT GDML writer. This would help to fix the problem.

saraheno commented 1 year ago

my apologies @MarkusFrankATcernch : I do not seem to understand properly how to do that command

image

andresailer commented 1 year ago

Hi @saraheno Run ddsim with --runType shell and put the command in the geant4 shell

Edit: fixed capitalisation

saraheno commented 1 year ago

When I did that, I got

image

the command I gave was ddsim --runtype shell

andresailer commented 1 year ago

--runType with capital T, sorry!

andresailer commented 1 year ago

You also need to give all the other arguments to build the detector, e.g., --compactFile

saraheno commented 1 year ago

do I did: ddsim --steeringFile SCEPCALsteering.py --compact ./DRSingleCrystal.xml --runType shell --part.userParticleHandler='' -G --gun.position="0.,30.,0." --gun.direction "0 -1 0" --gun.energy "1*GeV" --gun.particle="mu-" -N 1 -O out.root

and got

image

saraheno commented 1 year ago

nevermind I'll try again

saraheno commented 1 year ago

okay, it worked. The gdml file created that way is at: https://github.com/saraheno/SingleDualCrystal/blob/main/compact/DRSingleCrystal.gdml

MarkusFrankATcernch commented 1 year ago

@saraheno Thank you. I will check the output GDML. If things are correct, it should be the same as the corrected GDML from Hans. Do you have somewhere Hans' GDML in git ?

saraheno commented 1 year ago

@hanswenzel @MarkusFrankATcernch

I think the file Hans made by correcting my old file is at https://github.com/hanswenzel/CaTS/blob/main/gdml/DRSingleCrystal.gdml

or at

https://github.com/hanswenzel/CaTS/blob/main/gdml/sarah.gdml

hanswenzel commented 1 year ago

I can verify that the last gdml file that sarah provided: https://github.com/saraheno/SingleDualCrystal/blob/main/compact/DRSingleCrystal.gdml is a valid gdml file. you can find my versions in: https://github.com/hanswenzel/CaTS/blob/main/gdml/sarah.gdml

and https://github.com/hanswenzel/CaTS/blob/main/gdml/DRSingleCrystal.gdml I think now that the gdml file is valid one has to check that we get all the properties right (e.g. kill photons once they hit the photo detector)

Screenshot at 2022-12-07 07-57-35

saraheno commented 1 year ago

@hanswenzel when I look at this display, it looks like photons are escaping out the sides even though it is wrapped with a perfectly reflective surface? Is this what is happening?

MarkusFrankATcernch commented 1 year ago

I compared the 2 GDMLs:

hanswenzel commented 1 year ago

I think that's true . Just for fun I replaced the optical surface with one off a tyvek wrapped one. You will see the effect is much different. So what I can do is looking up the best surfaces to use. I guess the slices in the end are supposed to represent photo detectors so it doesn't make sense that they also radiate. 

<!--
<opticalsurface finish="0" model="0" name="_world_DRCrystal_mirrorSurface0xcf11350" type="1" value="1">

  <property name="REFLECTIVITY" ref="REFLECTIVITY0xcf376c0"/>
  <property name="EFFICIENCY" ref="EFFICIENCY0xcf377a0"/>

</opticalsurface>

--> Screenshot at 2022-12-07 08-37-29

chekanov commented 1 year ago

I've just made one small check by adding Al foil between 4 crystals (in depth): https://github.com/chekanov/SimpleECAL/blob/main/compact/ECalBarrel_DualCrystal.xml

            <layer repeat="4" vis="ECalLayerVis">
               <slice material = "Al"            thickness = "0.94*mm" vis="InvisibleNoDaughters"/>
               <slice material = "killMedia" thickness = "0.2" sensitive="yes" limits="cal_limits" vis="CrystalEcalSensitiveVis"/>
               <slice material = "E_PbWO4" thickness = "50.0*mm" sensitive="yes" limits="cal_limits" vis="CrystalEcalSensitiveVis" radiator="yes" />
               <slice material = "killMedia" thickness = "0.2" sensitive="yes" limits="cal_limits" vis="CrystalEcalSensitiveVis"/>
          </layer>

I have to kill the program after 1 min since my server becomes unresponsive. But it runs for 1 min without problems.. For 1 GeV muon, this does not have any sense. The problem disappears when I remove


     scint = PhysicsList(kernel, 'Geant4ScintillationPhysics/ScintillationPhys')
     scint.VerboseLevel = 0
     scint.TrackSecondariesFirst = True
     scint.enableUI()
     seq.adopt(scint)

from SCEPCALsteering.py. So, this problem is not related to the Cherenkov photons.

Comment from Sarah: Changing the constant SCINTILLATIONYIELD from "13.9/keV" to the value 0.00001 (arbitrary small) stops the program after 30 sec. This means the program cannot handle too many photons in these crystals.

MarkusFrankATcernch commented 1 year ago

13.9/keV = 13.9/1e-6 = 13900000 in TGeo (Energies in GeV) and becomes when translated to Geant4 (Energies in MeV): 13.9/TGeo::keV * (TGeo::keV / G4::keV) = 13.9/G4::keV = 13.9 / 1e-3 = 13900

This is quite different from 0.00001.

saraheno commented 1 year ago

My apologies Markus: I am missing your point? I have been telling Sergei that I don't think the code is in an infinite loop, that one just cannot simulate that many photons in that many crystals in a reasonable amount of time (without say the CaTs accelerator). I told him we could prove that by setting the number of photons made very low and then the program would finish an event.

so I think there is nothing really wrong with the code per say. It is just not reasonable without accelerators to put in the real number of photons. Either give up on the scintillation simulation of the crystals, or dial it down to something that can finish in a reasonable amount of time. I think giving up on simulating the scintillation is best. instead, use the energy deposits and then use a single crystal simulation to parameterize the number of collected photons as a function of the energy deposit size and position. That is what CMS does for example.

saraheno commented 1 year ago

also @MarkusFrankATcernch that reminds me: we are approaching the stage where "simhits" are not enough. we'll need rechits. Is there a tutorial or something that will help us understand the recommend code templates for this step?

MarkusFrankATcernch commented 1 year ago

@saraheno I am not sure about the program timing. I think there is nothing wrong with the code as well, but there is simply a lack of understanding what the various parameters really mean and what their value should be when they enter Geant4 -- for some I even had to guess the units for the conversion from TGeo to Geant4. If the parameters entering the game are far from what G4 expects, anything could go wrong.

For rechits or user defined hits, see examples/DDG4_MySensDet Have a look how the dictionary for the MyTrackerhits is build. This should be all you need to e.g. save a vector<MyTrackerHits(*)> to a root tree. There is also the possibility to add a data structure to the Geant4Hit structure. This structure of course when it comes to persistency must be embedded in a user-defined structure, but at least then this string dependency is localized.

MarkusFrankATcernch commented 1 year ago

There is more stuff coming up: I discovered that in Geant4 11 apparently FASTCOMPONENT and SLOWCOMPONENT are no longer material properties. In Geant4 10.07 they were still there...

hanswenzel commented 1 year ago

That is correct one has to be aware that names of some optical properties changed between Geant4 10 and Geant4 11. You can still write gdml files that would work with both versions (the keywords for the different versions would just be ignored). Look e.g. at the https://github.com/hanswenzel/CaTS/blob/main/gdml/simpleLArTPC.gdml for a gdml file that would work both with 10 and 11. And I agree it's a pain in the bud.

saraheno commented 1 year ago

@MarkusFrankATcernch @hanswenzel

Hans mentioned that one can set visual attributes to optical surfaces. Is this possible to do in dd4hep? I don't see this in the OpNovice example?