MultipoleMagnet rotation issue

[aciarma]

When rotating the multipole, the result field map is wrong (unphysical). In the following example I build a cylindrical quadrupole along the Z-axis, and rotate it on the Y-axis of multiples of pi/2, and then plot the resulting field map obtained using dumpBfield.

As it can be seen from the plot there are unphysical fields along the length of the quadrupole (1/2pi and 3/2pi rotations), and the field of the pi-rotated quadrupole also shows unrealistic "circular" field lines.

It seems like after the volume rotation, also each "vector" gets its Y-component added of the rotation angle. Could this please be checked? It is very important to have the possibility to correctly implement magnetic elements in the models.

I am working on lxplus, and I get this behaviour both with the key4hep latest and nightly releases.

<fields>
  <field name="QC1L1_field_0" type="MultipoleMagnet" Z="0.0*tesla">
      <position y="0*cm" x="5*cm" z="0*cm"/>
      <rotation x="0" y="0" z="0.0"/>
      <coefficient coefficient="0*tesla"/>
      <coefficient coefficient="(-1)*(45.6)*(-0.273)/0.3*tesla/m"/>
      <shape type="Tube" rmin="0.*cm" rmax="QC1_rmin" dz="2*cm" />
  </field>

  <field name="QC1L1_field_1" type="MultipoleMagnet" Z="0.0*tesla">
      <position y="5*cm" x="0*cm" z="0*cm"/>
      <rotation x="0" y="pi/2." z="0.0"/>
      <coefficient coefficient="0*tesla"/>
      <coefficient coefficient="(-1)*(45.6)*(-0.273)/0.3*tesla/m"/>
      <shape type="Tube" rmin="0.*cm" rmax="QC1_rmin" dz="2*cm" />
  </field>

  <field name="QC1L1_field_2" type="MultipoleMagnet" Z="0.0*tesla">
      <position y="0*cm" x="-5*cm" z="0*cm"/>
      <rotation x="0" y="pi" z="0.0"/>
      <coefficient coefficient="0*tesla"/>
      <coefficient coefficient="(-1)*(45.6)*(-0.273)/0.3*tesla/m"/>
      <shape type="Tube" rmin="0.*cm" rmax="QC1_rmin" dz="2*cm" />
  </field>

  <field name="QC1L1_field_3" type="MultipoleMagnet" Z="0.0*tesla">
      <position y="-5*cm" x="0*cm" z="0*cm"/>
      <rotation x="0" y="pi*3./2." z="0.0"/>
      <coefficient coefficient="0*tesla"/>
      <coefficient coefficient="(-1)*(45.6)*(-0.273)/0.3*tesla/m"/>
      <shape type="Tube" rmin="0.*cm" rmax="QC1_rmin" dz="2*cm" />
  </field>

</fields>

[andresailer]

Hi Andrea,

Can you explain what you expect to see?

[aciarma]

Hi Andre,

the shape of the quadrupole field lines should be "hyperbolic", like what you see in the 'no rot' case. But after the rotation of pi of the element they look like a circle, which is not correct.

I mean, if I take a quadrupole and go around it, the shape of the field lines cannot change! only the direction does of course.

[andresailer]

Can you make a sketch with just a few arrows/lines showing where things should be pointing?

[aciarma]

should look something like this

[andresailer]

Thanks, now I understand what you are pointing out!

[MarkusFrankATcernch]

@aciarma Could you please post the ROOT script to draw the field lines ? Thanks a lot!

[aciarma]

Hi Markus,

I use gnuplot to draw the field lines from dumpBField output (delete or comment non-data lines):

set term postscript enhanced color 20 solid
set output 'field.eps'
set xtics font ",25"
set ytics font ",25"

set palette rgb 33,13,10;
set xlabel "x [cm]" font ",25" offset -1
set ylabel "y [cm]" font ",25" offset -1

z_start = 0

#norm factor for arrows length
an = 10

set xrange [-8:8]
set yrange [-8:8]

#set label "e+ downstream" at -5,3
#set label "e- upstream"   at +2,3
set label "pi rot" at -6,3
set label "no rot"   at +4,3
set label "1/2 pi rot" at +3,+7
set label "3/2 pi rot" at +3,-7
set arrow from -8,-0 to 8,0 nohead lw 1 lt "."
set arrow from -0,-8 to 0,8 nohead lw 1 lt "."
#set title "B_T[T], at z = +2.3 m (QC1L1)"
set title "B_T[T], at z = 0m"
plot 'dumpBField.dat' using 1:(($4!=0&&$5!=0&&$3==z_start)?$2:1/0):($4/an):($5/an):(sqrt($4*$4+$5*$5)) notitle with vectors head filled lc palette lw 3
unset label
unset arrow

!ps2pdf field.eps field.pdf
!rm -f field.eps

[MarkusFrankATcernch]

@aciarma Is this the plot you are looking for?

[aciarma]

Hi @MarkusFrankATcernch, it is quite hard to see the arrows but it seems like a yes! Is there a PR I can try just to be sure?

[MarkusFrankATcernch]

I have created a branch containing the changes: https://github.com/MarkusFrankATcernch/DD4hep/tree/BFieldTrafos

If you could verify against it would be of great help! Thanks a lot!

[aciarma]

Hi Markus, the output of dumpBfield now looks something like this:

#######################################################################################################
      x[cm]            y[cm]            z[cm]          Bx[Tesla]        By[Tesla]        Bz[Tesla]     
Pos: -8.00000000e+00  -8.00000000e+00  -1.00000000e+01  Inverse: -1.30000000e+01  -8.00000000e+00  -1.00000000e+01
Pos: -8.00000000e+00  -8.00000000e+00  -1.00000000e+01  Inverse: +1.00000000e+01  -1.30000000e+01  -8.00000000e+00
Pos: -8.00000000e+00  -8.00000000e+00  -1.00000000e+01  Inverse: +3.00000000e+00  -8.00000000e+00  +1.00000000e+01
Pos: -8.00000000e+00  -8.00000000e+00  -1.00000000e+01  Inverse: -1.00000000e+01  -3.00000000e+00  +8.00000000e+00
 -8.00000000e+00  -8.00000000e+00  -1.00000000e+01  +0.00000000e+00  +0.00000000e+00  +0.00000000e+00  
Pos: -8.00000000e+00  -8.00000000e+00  +0.00000000e+00  Inverse: -1.30000000e+01  -8.00000000e+00  +0.00000000e+00
Pos: -8.00000000e+00  -8.00000000e+00  +0.00000000e+00  Inverse: -4.89858720e-16  -1.30000000e+01  -8.00000000e+00
Pos: -8.00000000e+00  -8.00000000e+00  +0.00000000e+00  Inverse: +3.00000000e+00  -8.00000000e+00  +9.79717439e-16
Pos: -8.00000000e+00  -8.00000000e+00  +0.00000000e+00  Inverse: -4.89858720e-16  -3.00000000e+00  +8.00000000e+00
 -8.00000000e+00  -8.00000000e+00  +0.00000000e+00  +0.00000000e+00  +0.00000000e+00  +0.00000000e+00  
...

what do Pos: and Inverse: mean? Anyway, if I ignore those lines I get a plot similar to yours, which is what I expect.

In order to further test this I'd like to check the effect of this field during particle tracking. I am trying to track through the field few test particles using k4run, but I can't figure out how to save the path of the particles (e.g. every 5mm), SimG4SaveTrajectory only gives me the starting points apparently... I do something like this:

geantservice = SimG4Svc("SimG4Svc", detector='SimG4DD4hepDetector', physicslist="SimG4FtfpBert")
geantservice.g4PostInitCommands += ["/tracking/storeTrajectory 1"]

from Configurables import SimG4Alg, SimG4PrimariesFromEdmTool, SimG4SaveTrajectory

savetrajectorytool = SimG4SaveTrajectory("saveTrajectory")
savetrajectorytool.TrajectoryPoints.Path = "trajectoryPoints"

particle_converter = SimG4PrimariesFromEdmTool("EdmConverter")
particle_converter.GenParticles.Path = "allGenParticles"

geantsim = SimG4Alg("SimG4Alg",
                    outputs= [
        "SimG4SaveTrajectory/saveTrajectory",
        ],
                    eventProvider=particle_converter)

Is there another way to save this information?

[MarkusFrankATcernch]

@aciarma Sorry for the noise in the output....these are leftovers from debugging. I will remove them.

[MarkusFrankATcernch]

@aciarma ....in principle you can save anything to the root file....root is patient. ;-)

But to be more concrete: The Geant4Particle has a field:

      /// User data extension if required
      dd4hep_ptr<ParticleExtension> extension;   //! not persistent. ROOT cannot handle

If ROOT understands now std::unique_ptr<T> (I do not know. It did not at a time), you could remove the "!" from the comment and attach a user object with the proper information. Of course you will have to provide dictionaries for these classes in order to make ROOT happy. This should allow you to save supplementary info, but it is not a out-of-the-box solution.

[aciarma]

That's the point, I am surprised that there is not a "ready to use" way to have truth particle information in ddsim/k4run environments? Or maybe I misunderstood?

[MarkusFrankATcernch]

@aciarma 1) There is truth handling present. It is done in the Geant4ParticleHandler. However, there is currently not a possibility to add user defined information to particles. I should probably re-test if ROOT knows now to handle the persistency of std::unique_ptr 2) Could you please check the head ? I should have removed the unnecessary printout.

[aciarma]

@MarkusFrankATcernch the output is nice and clean now, thank you!

I'll try to look better into the truth particles handling

[andresailer]

Hi @aciarma

You can make whatever volume the particle fly through sensitive, set the max step length to 5*mm in the region, which then gives you the position and momentum at these hits.

Hot to make the volume sensitive depends on what geometry you are using.

[aciarma]

Hi @andresailer I tried using the Box_simpleTrackerSD example, and set in the steering file as following, but it did not work. Is this what you meant by setting the step length?

...
regiontool = SimG4UserLimitRegion("limits", volumeNames=["BoxTracker"],maxStep = 5*mm)
geantservice = SimG4Svc("SimG4Svc", detector='SimG4DD4hepDetector', physicslist="SimG4FtfpBert",
                        #regions=[regiontool],
                        )

from Configurables import SimG4Alg, SimG4SaveTrackerHits, SimG4SaveCalHits, SimG4PrimariesFromEdmTool, SimG4SaveParticleHistory, SimG4SaveTrajectory

savetrackertool = SimG4SaveTrackerHits("saveTrackerHits_box", readoutNames = ["TrackerSiliconHits"])
savetrackertool.SimTrackHits.Path = "TrackerSiliconHits"

particle_converter = SimG4PrimariesFromEdmTool("EdmConverter")
particle_converter.GenParticles.Path = "allGenParticles"

geantsim = SimG4Alg("SimG4Alg",
                    outputs= ["SimG4SaveTrackerHits/saveTrackerHits_box"],
                    eventProvider=particle_converter)
...

[andresailer]

Sorry, but the k4SimGeant4 configurations are not my area of expertise.

I was thinking of something like this (and how it is used in the rest of the XML)

https://github.com/AIDASoft/DD4hep/blob/04ce7bb99fc74e235c6e9d3966124d7b68c91124/examples/ClientTests/compact/SiliconBlock.xml#L54-L57

Only step_length_max, do not set track_length_max

[MarkusFrankATcernch]

I assume this to be fixed.