Primary particle and a new input mode

[benkrikler]

This is to document a couple of recent additions I've put in:

• Added support for a 'histo' mode for position where you provide a 3D histogram (in global coordinates) ad we generate randomly according to this.
• Store the primary particle's momentum, position and name (proton, mu- etc). These can be controlled through the EventHeader_section in the output macro. The additional flags are 'position', and 'pid' (I've added these to the output_default macro). There was already support for the primary momentum but it seemed to only be meaningful if we read particles in from a root tree. This should now work for all methods of primary particle input.

Hopefully I've not stuck in any bugs but if anyone does use these modes let me know...

Out of curiosity, @thnam and others ( @AndrewEdmonds11 and I have spoken a bit on email) how were you obtaining the primary particle distributions and reading in positions from histograms before? I thought we'd been doing this but the code didn't seem to support it so I guess I mis-understood what was going on.

[AndrewEdmonds11]

I think there may be a problem with the "histo" mode in that the initial proton positions in local coordinates are quite strange.

In the histogram file I use to sample initial positions I have 3D histograms of the muon stopping positions in both local and global coordinates.

The Z-projection of the local coordinate plot looks as expected: .

Projecting the 3D global stopping positions histogram, gives me this:

However, when I look at the initial positions in my proton run I get (in global coordinates): . This is worse binning but looks fairly consistent with the previous plot.

However, when I rotate this to the local Z coordinate I get: (NB I don't know how to rotate the coordinates in the 3D histogram) which looks like some of my protons may be starting outside of the target.

My thinking is that it's either a problem with the global muon stopping position coordinates I'm recording or it's to do with how ROOT gets a random number from a TH3. Any ideas?

[benkrikler]

Part of that shape could be from Root linearly interpolating between bins when it makes the random numbers which gives the triangular shape. Do you know the size of the bins in the 3D plot?

Also, you put the units in micrometres, but they're in cm in the projection of the 3D input histogram and you multiply by 1e3 in the formula, so are the units correct or should it be mm?

In the same way, I notice the mean of your final plot is about a tenth the size of the mean of the first plot, though that could be complete coincidence.

My guess, is that we have a problem with the unit of the output file from the monitors being different to the units used for the input file. Quite why we get this exact shape I don't get though...

[AndrewEdmonds11]

The bin width is 0.1 cm in the 3D plots so I should probably decrease that.

So I get a bit confused with the units. The branches i_x etc. want to be read in in cm, right? I write them out to the 3D distribution as cm. I then wanted the final plot in um, so I should have multiplied by 10^4? This then makes the base of the triangle go between +-700 um, which is even worse!

I'll try decreasing the bin width and see what happens....

[AndrewEdmonds11]

Decreasing the bin width by a factor of 10 obviously increases the number of bins by a factor of 1000 and causes ROOT to crash. I think this was a problem I was having before.

I've tried using THnSparse but it seems to take a long time to run (it's a lot better when I compile the ROOT macro). But this would mean we would have to change the code to use THnSparse->GetRandom(Double_t* rand) rather than TH3->GetRandom3(Double_t& x, Double_t& y, Double_t& z).

The other option (which I think I ended up implementing somewhere) is to just check that the primary vertex is in the Target volume but I seem to remember this slowed the simulation down a bit.

I'll keep plugging away at the first one for the time being and try not to break the code for anyone else...

[benkrikler]

The other thing is to swap to using the root mode which uses a root tree as an input. Then you can load in the positions directly and not worry about this binning.

It's tough to explain here, but I'm thinking that the issue is because the bins are rotated with respect to the target and fairly coarse. If you take the 2nd 2D plot and imagine projecting it back to a plane at 135degrees (so the plane where x=-z), because the bins are rotated 45 degrees compared to this, I think you'll get this funny triangular shape. On the y=x line, all the bins are very full. But once I come away from this line so that now y=x+a where a is less than bin_width*sqrt(2), I start overlapping the full bins with the emptier bins. Since the way the ratio of the contributions from full bins to empty bins changes is linear with a I would get a triangular shape.

If that explanation is right, then this isn't a bug so to speak, but it might suggest this technique is not very good. An improvement would be to use a local histogram and move and rotate it to be back at the target. I'm not sure how easy that is to implement though and I'm not sure if we really need it. So if you are able to switch to the root tree approach then perhaps that's worth doing to avoid this issue complicating the uncertainties at this time.

[AndrewEdmonds11]

So I've got it working with THnSparse :smile:.

Here's the randomly generating global proton positions: which looks very similar to what comes out of the initial MC.

Here are the diffs for PrimaryGeneratorAction.cc and .hh:

diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc
index 9801c53..0e7b253 100644
--- a/src/PrimaryGeneratorAction.cc
+++ b/src/PrimaryGeneratorAction.cc
@@ -188,8 +188,8 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
        }

        if ( EnergyMode == "histo"){
-               G4double mom = EM_hist->GetRandom() * MeV;
-               G4double ekin = sqrt(mom*mom+mass*mass)-mass;
+               G4double ekin = EM_hist->GetRandom() * 1e3 * MeV;
+               //              G4double mom = sqrt(ekin*ekin + 2*ekin*mass);
                particleGun->SetParticleEnergy(ekin);
        }
        else if ( EnergyMode == "root" ){
@@ -531,10 +531,18 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
          SetRandomPosition();  
        }
        else if ( PositionMode == "histo") {
+         if (PM_hist) { // if we have a TH3F
                double x=0,y=0,z=0;
                PM_hist->GetRandom3(x,y,z);
                G4ThreeVector pos_3Vec(x, y, z);
                particleGun->SetParticlePosition(pos_3Vec);
+         }
+         else if (PM_hist_sparse) {
+           double random_pos[3];
+           PM_hist_sparse->GetRandom(random_pos);
+           G4ThreeVector pos_3Vec(random_pos[0], random_pos[1], random_pos[2]);
+           particleGun->SetParticlePosition(pos_3Vec); // should probably try and avoid having same line twice
+         }
        }
        else if ( PositionMode == "turtle" || PositionMode == "muPC" || PositionMode == "collimator") {
          // Already handled in the DirectionMode if block
@@ -863,13 +871,20 @@ void PrimaryGeneratorAction::BuildHistoFromFile(){
                                        "InvalidInput", FatalException,
                                        "Can not find file");
                }
-               if(!obj->InheritsFrom("TH3")){
+               if(obj->InheritsFrom("TH3")){
+                 PM_hist = (TH3*) obj;
+               
+               }
+               else if (obj->InheritsFrom("THnSparse")) {
+                 PM_hist_sparse = (THnSparseF*) obj;
+                 // should also check that it's 3dimenstional here
+               }
+               else {
                        std::cout<<"ERROR: Cannot generate positions from non-3D histogram: "<<PM_hist_histname<<"!!!"<<std::endl;
                        G4Exception("PrimaryGeneratorAction::BuildHistoFromFile()",
                                        "InvalidInput", FatalException,
                                        "Wrong hist type");
                }
-               PM_hist = (TH3*) obj;
        }

 }

and

diff --git a/include/PrimaryGeneratorAction.hh b/include/PrimaryGeneratorAction.hh
index 537bbb9..e832301 100644
--- a/include/PrimaryGeneratorAction.hh
+++ b/include/PrimaryGeneratorAction.hh
@@ -18,6 +18,7 @@
 #include "TF1.h"
 #include "TH2F.h"
 #include "TH3F.h"
+#include "THnSparse.h"

 class G4ParticleGun;
 class G4Event;
@@ -134,6 +135,7 @@ class PrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction,
                G4String                   DM_hist_filename;
                G4String                   DM_hist_histname;
                TH3*                       PM_hist;
+                THnSparseF*                PM_hist_sparse;
                G4String                   PM_hist_filename;
                G4String                   PM_hist_histname;

I think I've coded it in such a way that it doesn't break anything so should I just commit this?