Changes to remove the odd entries in the response matrix

[AndrewEdmonds11]

I've finally merged in the changes I made to remove the odd entries in the response matrix (see elog:252 for an explanation and the commit is 8562b04).

Here's the diff in MonitorSD.cc where I just wondered why we weren't also considering stopped particles. (Admittedly, I was suffering from a pretty bad cold at the time so I could be complete wrong).

diff --git a/src/MonitorSD.cc b/src/MonitorSD.cc
index 8d7537d..3d71d2c 100644
--- a/src/MonitorSD.cc
+++ b/src/MonitorSD.cc
@@ -712,7 +712,7 @@ G4bool MonitorSD::ProcessHits(G4Step* aStep,G4TouchableHistory* touchableHistory
 //                     std::cout<<"dt too small, will not push"<<std::endl;
                        willPush = false;
                }
-               if ( m_tid[index] == trackID&& prePoint->GetStepStatus() != fGeomBoundary || killed){ // If this particle was in this volume in last s
+               if ( m_tid[index] == trackID && (prePoint->GetStepStatus() != fGeomBoundary || killed || stopped)){ // If this particle was in this vo
 //                     std::cout<<"Was here last step"<<std::endl;
                        willPush = false;
                }

The other change was to set the tres parameter in output_default to -1 so that no hits are merged.

[benkrikler]

[ Sorry in advance, this post turned into a bit of an essay. It's just that this is important for all the MC outputs. ]

This whole block of code is a bit tricky to get my head round but I think I'm getting there.

First we try to find a previous hit by the same track:

        bool willPush = true;
        //look for last hit in the same volume.
        G4int index = -1;
        G4double dt = 1e10*s;
        for ( int i = 0; i < nHits; i++ ){
                if ( m_volName[i] == VolName && m_volID[i] == ReplicaNo ){
                        G4double dtime = pointOut_time - m_t[i]*unit_t;
                        if ( fabs(dtime) < fabs(dt) ){
                                index = i;
                                dt = dtime;
                        }
                }
        }

So by the end of that, index is either -1 meaning this is the first hit in this volume in this event, or it's equal to the array index of the last hit that:

• Happened in this volume (and replica ID in case we have repeated geometries. We don't in alcap so this is not important here).
• Occurred most recently compared to the new hit we want to add.

Does that sound right?

Next we ask if we really do want to make a new hit or combine this new one with the previous one we've just found.

        if ( index != -1 ){// found a previous hit in the same volume
                if ( fabs(dt) < tres ){// dt too small, will not push
                        willPush = false;
                }
                if ( m_tid[index] == trackID && (prePoint->GetStepStatus() != fGeomBoundary || killed || stopped)){ // If this particle was in this volume in last step, don't generate a new hit
                        willPush = false;
                }
        }

So if we did find such a previous hit (labelled by index) we only make a new one if:

1. the time difference between this new hit and the previous one (index) is greater than tres.

2. and

   a. either the track ID of this new hit is different to the one that made the previous one at index
   b. or the step that made the new hit straddled the volume's boundary and we haven't killed the track and the track didn't stop.

If you set tres to -1 in the config then you completely remove the first condition and only rely on the second. If the track ID is different it makes sense not to combine the hits (though in a real detector this would be pile-up I don't think we want that level of detail at this point, we would add it later on in the simulation).

I'm a bit unsure why the 2.b. criteria are right though. I suppose if a step straddles the volume boundary the previous hit couldn't have come from the same source as this new hit (assuming the lists are chronologically ordered which I think is right to assume). Why do we combine this hit with a previous one if the track was killed or stopped though? I'm not sure I see why that's related to this issue?

I think as well that these criteria for whether we combine the hit should partly be merged into the first loop to look for the previous hit (and set the value of index). We could check the Track ID to be the same in the loop to set index since we have all that information already. If the track ID's are different then we should never combine hits. And in this current setup, if we had two previous hits in the same volume, one with the same track ID the other with a smaller time difference, we'll pick the one with the smaller time difference to consider for merging, when it's probably more important to consider the one with the same track ID since we're not trying to simulate pile-up at this stage.

Chen's original criteria for 2. was:

  if ( m_tid[index] == trackID&& prePoint->GetStepStatus() != fGeomBoundary || killed){ // If this particle was in this volume in last step, don't generate a new hit

Which means we make a new hit if:

1. the track ID was different to the previous hit or the step straddles a geometric boundary
2. and the track wasn't killed

Did you intentonally add the brackets around the parts after the first && in that line? I'm guessing you're aware, but the && takes precedence over ||, so Chen's original implementation was more like:

 (m_tid[index] == trackID && prePoint->GetStepStatus() != fGeomBoundary) || killed

[ true && true || false && false is different to true && (true || false) && false ]

Again, sorry for the loooonng response. I hope this all makes sense...

[AndrewEdmonds11]

Thanks, Ben. It's useful to have it all thought through.

So by the end of that, index is either -1 meaning this is the first hit in this volume in this event, or it's equal to the array index of the last hit that:

• Happened in this volume (and replica ID in case we have repeated geometries. We don't in alcap so this is not important here).
• Occurred most recently compared to the new hit we want to add.

Does that sound right?

Yes, that's what I thought.

So if we did find such a previous hit (labelled by index) we only make a new one if:

• the time difference between this new hit and the previous one (index) is greater than tres.
• and
  a. either the track ID of this new hit is different to the one that made the previous one at index
  b. or the step that made the new hit straddled the volume's boundary and we haven't killed the track and the track didn't stop.

I think you need to swap and/or in b, right?

Why do we combine this hit with a previous one if the track was killed or stopped though?

I'm not sure what the original reason was but isn't this just to ensure that we don't lose the last step. i.e. if you have a track with 10 steps in a volume and the 10th one is when it's killed/stopped, you'll want to add it to the previous steps, right?

We could check the Track ID to be the same in the loop to set index since we have all that information already. If the track ID's are different then we should never combine hits.

Agreed.

Did you intentionally add the brackets around the parts after the first && in that line?

Yes

I'm guessing you're aware, but the && takes precedence over ||

No I wasn't aware. I just put the brackets in to show more easily how I thought the logic was supposed to be (I blame the fact that I was ill 😄)

Evidently, this part of the code will need some refactoring before the second run but can we draw any simple conclusions out of this? What needs changing to get things working for the time being?

[benkrikler]

So if we did find such a previous hit (labelled by index) we only make a new one if:

 the time difference between this new hit and the previous one (index) is greater than tres.
 and a. either the track ID of this new hit is different to the one that made the previous one at index b. or the step that made the new hit straddled the volume's boundary and we haven't killed the track and the track didn't stop.

I think you need to swap and/or in b, right?

I don't think so, since I'm phrasing it under what conditions do we make a new hit whereas the contents of the if block is to check the opposite, when we should combine two hits. !( A && B) = !A || !B

Why do we combine this hit with a previous one if the track was killed or stopped though?

I'm not sure what the original reason was but isn't this just to ensure that we don't lose the last step. i.e. if you have a track with 10 steps in a volume and the 10th one is when it's killed/stopped, you'll want to add it to the previous steps, right?

I'm not sure I understand why Chen did this originally either, but in your new logic if a track is killed or stopped it is always combined to the previous hit (index). I'm not sure why this matters, though I suppose if the other criteria are correct (the relative timing and that index really does point to the last hit) then these additional critieria won't have any impact since the point should be merged anyway.

Did you intentionally add the brackets around the parts after the first && in that line?

Yes

I'm guessing you're aware, but the && takes precedence over ||

No I wasn't aware. I just put the brackets in to show more easily how I thought the logic was supposed to be (I blame the fact that I was ill 😄)

So just to be clear, you were aware that you're changing the logic from what Chen had not just improving legibility? Originally it was:

 (m_tid[index] == trackID && prePoint->GetStepStatus() != fGeomBoundary) || killed

Whereas now we have:

 m_tid[index] == trackID && (prePoint->GetStepStatus() != fGeomBoundary || killed || stopped)

If you ignore the fact we now include stopped in the decision, before the hits were only merged if [the track was killed] or [ the track IDs were the same and the step didn't straddle the edge of the volume [. Now they're merged if they're [stopped or killed or don't straddle a boundary] and they have the same track ID. The more I think about it, I think your way is the right one, though I'm not sure since I'm not sure killed or stopped should be included anyway.

Evidently, this part of the code will need some refactoring before the second run but can we draw any simple conclusions out of this? What needs changing to get things working for the time being?

I would replace the entire block so that we check track IDs in the loop to find the previous hit and then drop killed and stopped from the criteria. So my final implementation for the entire block of code would be:

bool willPush = true;
// Check if there was a previous hit that we should be merging this one against
G4int index = -1;
G4double dt = 1e10*s; 
for ( int i = 0; i < nHits; i++ ){
    if (   m_tid[i] == trackID     // Must have the same track ID
        && m_volName[i] == VolName // Must be in the same logical volume
        && m_volID[i] == ReplicaNo // If we have repeated physical volumes, should have the same replica number
       ){
        // Take the absolute difference between the times of the new hit and the current one we're considering
        G4double dtime = fabs(pointOut_time - m_t[i]*unit_t);  

        // If this is less than the current closest candidate hit, then make this one the new candidate
        if ( dtime < dt ){
            index = i;
            dt = dtime;
            // Set the flag to make a new hit to false if `dt` is less than `tres`
            if( willPush && dt < tres) willPush=false;
        }
    }
}

What do you think?

[AndrewEdmonds11]

I don't think so, since I'm phrasing it under what conditions do we make a new hit whereas the contents of the if block is to check the opposite, when we should combine two hits. !( A && B) = !A || !B

Yes, you're correct - my mistake

though I suppose if the other criteria are correct (the relative timing and that index really does point to the last hit) then these additional critieria won't have any impact since the point should be merged anyway.

Good point.

So just to be clear, you were aware that you're changing the logic from what Chen had not just improving legibility?

No I had no idea what I was doing. I just assumed the things being OR'd together were meant to be together.

What do you think?

I like it. This removes the need of the tres parameter then?

[benkrikler]

This removes the need of the tres parameter then?

No it's still there, in the last line of actual code. We only set willPush to false if dt is less than tres. In the end, I've actually made it so that tres is more important since setting it to -1 will stop the merging altogether right now. If you want to keep that functionality we could change that line to:

if( willPush && (dt < tres || tres == -1 ) ) willPush=false;

or more simply, set tres to be huge in the ReadOutputCard method if a user gives a value equal to -1. That would be slightly more efficient, since we'd only do that check once.

[AndrewEdmonds11]

Ok. Hang on, I've got myself a bit confused here. What's the tres parameter for now?

Originally, I thought it was so that steps from different tracks (i.e. different particles) would be combined together so that we simulate pile-up in the actual detector. But now, since we are inside an if-statement that checks the hits are from the same tracks, what's it for?

[benkrikler]

tres is the maximum time between two hits before you add them together. If one track were to pass through a volume, deposit energy then come back some time later (like it spiralling back in a magnetic field) then you should probably make two separate hits. Or if you have a very slow moving particle that deposits it's energy very slowly, it may be realistic to pair up the energy deposits into separate hits. The tres parameter allows you to control how far apart (in time) the energy deposits can occur for you to combine them together.

[AndrewEdmonds11]

Ok, thanks for clearing that up. I guess for AlCap this isn't so important, right?

Let me know when you've finished implementing this so I can pull in the changes.