M_stopped branch not working?

[AndrewEdmonds11]

So I'm trying to look get the number of stopped protons in the thick silicon detectors and am using the condition:

tree->Draw("M_particleName>>hist3", "M_particleName==\"proton\" && (M_volName == \"ESi1\" || M_volName == \"ESi2\") && M_stopped == 1")

and I get 0 protons.

If I invert the stopping condition then I get 2844 protons.

However, if I look at the number of protons that are in the thick silicon and the corresponding veto, I get:

N(ESi1) == 1446
N(VetoPS1) == 610
N(ESi2) == 1398
N(VetoPS2) == 629

which would make me think that some of the protons are stopping. (I think there may be a geometrical effect where protons might pass through the silicon and miss the veto but I naively think that's tiny).

Also, if I look at the particles that are stopping in ESi1, then I get 28 e+ and 10 tritium.

How do we decide whether a particle has stopped?

I'll have a look at the code to see what's going on here.

[AndrewEdmonds11]

So it seems the track is "stopped" if the track status is fStopButAlive (see here).

Doing a 1000 muon run with verbose output, there is one proton that enters the ESi1:

*********************************************************************************************************
* G4Track Information:   Particle = proton,   Track ID = 10,   Parent ID = 1
*********************************************************************************************************

#Step# X         Y         Z         t         pX        pY        pZ        KineE     dEtot     dEIon     dEdel     StepLeng   TrakLeng   Volume     Process   
0      -3.15 cm  -5.33 mm  -3.15 cm  0     ps  -146  MeV -5.42 MeV 22.5  MeV 11.6  MeV 0     eV  0     eV  0     eV  0      fm  0      fm  Target    initStep
***********************************************************
1      -3.16 cm  -5.33 mm  -3.15 cm  2.03  ps  -141  MeV -5.15 MeV 22.2  MeV 10.8  MeV 810   keV 810   keV -810  keV 94.8   um  94.8   um  Target    UserLimit
***********************************************************
2      -3.58 cm  -5.48 mm  -3.09 cm  96.9  ps  -141  MeV -5.15 MeV 22.2  MeV 10.8  MeV 0.000137 eV  0.000137 eV  -0.000137 eV  4.27   mm  4.36   mm  TargetMount   Transportation
***********************************************************
3      -11.8 cm  -8.48 mm  -1.8  cm  1.94  ns  -141  MeV -5.15 MeV 22.2  MeV 10.8  MeV 0.00266 eV  0.00266 eV  -0.00266 eV  8.3    cm  8.73   cm  ChamberIn   Transportation
***********************************************************
4      -11.8 cm  -8.48 mm  -1.79 cm  1.94  ns  -137  MeV -3.8  MeV 24.6  MeV 10.3  MeV 472   keV 472   keV -472  keV 65.9   um  8.74   cm  dESi1    Transportation
***********************************************************
5      -11.9 cm  -8.51 mm  -1.78 cm  1.96  ns  -137  MeV -3.8  MeV 24.6  MeV 10.3  MeV 3.38e-05 eV  3.38e-05 eV  -3.38e-05 eV  1.02   mm  8.84   cm  DetBlock1   Transportation
***********************************************************
6      -11.9 cm  -8.51 mm  -1.77 cm  1.97  ns  -119  MeV -3.03 MeV 17.2  MeV 7.65  MeV 2.63  MeV 2.63  MeV -2.63 MeV 299    um  8.87   cm  ESi1     Transportation
***********************************************************
7      -11.9 cm  -8.52 mm  -1.77 cm  1.98  ns  -97   MeV -4.73 MeV 8.53  MeV 5.05  MeV 2.61  MeV 2.61  MeV -2.61 MeV 232    um  8.9    cm  ESi1     hIoni     
***********************************************************
8      -11.9 cm  -8.53 mm  -1.77 cm  1.98  ns  -62.4 MeV -946  keV 5.27  MeV 2.09  MeV 2.96  MeV 2.96  MeV -2.96 MeV 168    um  8.91   cm  ESi1     hIoni     
***********************************************************
9      -12   cm  -8.53 mm  -1.77 cm  1.99  ns  -0    eV  0     eV  0     eV  0     eV  2.09  MeV 2.09  MeV -2.09 MeV 51.2   um  8.92   cm  ESi1     hIoni     

Uncommenting the print outs in MonitorSD.cc, I found that the track status is actually fStopAndKill.

Since we also keep track of particles that are killed with M_killed, I propose the following correction:

        bool stopped = false;
        double stop_time = 0;
        bool killed = false;
        double kill_time = 0;
-        if (fTrackStatus == fStopButAlive){
+        if (fTrackStatus == fStopButAlive || fTrackStatus == fStopAndKill){
                stopped = true;
                stop_time = pointIn_time;
        }
-        else if (fTrackStatus == fStopAndKill || fTrackStatus == fKillTrackAndSecondaries){
+       if (fTrackStatus == fStopAndKill || fTrackStatus == fKillTrackAndSecondaries){
                killed = true;
                kill_time = pointIn_time;
        }

Would this be OK, Chen? If so, I can test and commit it tomorrow

[AndrewEdmonds11]

So, I've made the above change and it seems to work. I get 1557 protons stopping in either ESi1 or ESi2 from 1M initial target muons.

I've committed this to the branch I'm currently working on and so will be merged into master when I finish.

[thnam]

I also wonder how to get stopping position (muon on SiR2). I tried:

     if (M_volName->at(ihit) == "SiR2" && M_stopped->at(ihit)==1 &&
         M_particleName->at(ihit) == "mu-"){
       double gx = M_x->at(ihit);
       double gy = M_y->at(ihit);
       double gz = M_z->at(ihit);
      ....

       TVector3 pos(gx, gy, gz);
       pos.RotateY(TMath::Pi()/4); // rotate 45 deg about Y-axis
       }

But the resulted z-distribution does't look right: If I used McTruth_x instead, then there is some more fluctuation which is more reasonable:

And still there is some offset of z that I'm not sure if it can be calculated automatically. Or should we store both global and local coordinates of each hit?

[AndrewEdmonds11]

So apparently we can convert from global to local position (see the Geant4 FAQ).

Looking in MonitorSD.cc, we only have calls to preStepPoint->GetPosition(), which is global. It looks fairly easy to add so I will open another issue.

[wuchen1106]

Sorry for my late response! Andy, you made a good point! I was using M_stopped and M_killed separately but your idea is more reasonable. Nam, you should look at M_Ox which records the post step point. M_x is pre step point which is basically where the particle enters the volume.