nPrimaryParticles is wrong when running with single pion generator

[WenjieWu-Sci]

Run the program with single pion (pi+) generator (macros/single_particle/LAr_piplus.mac), the output is

number of primary vertices  : 1
Recorded tracks       : 24463
Tracks from FSL       : 24463
Tracks from FSL (2nd) : 0
number of primary particles : 9 , in which number of particles from fsl : 0
       PDG       Angle  TrackLength ShowerLength  ShowerWidthInLAr      EInLAr   EInHadCal   dEdxInLAr ProngType          Pz
       211       0.000     1617.166     6340.677           295.329   16099.895     113.337       2.683         0   25005.181
        22       0.474      275.084     1627.162            37.134     482.561       0.000       0.297         4     429.364
        22       0.605      213.532     1547.823            48.047     329.927       0.000       0.213         4     290.211
        22       0.423      972.331     1105.041            54.831     158.629       0.000       0.144         4     144.650
        22       0.686      237.900      858.349            25.740     138.116       0.000       0.161         4     106.894
        22       0.386      192.785     2047.834            62.943     619.644       0.000       0.303         4     576.882
        22       0.073      226.138     1154.668            28.369     147.386       0.000       0.128         4     146.999
        22       0.061      351.664     2193.675            54.501     842.286       0.000       0.384         4     840.722
        22       0.314      170.231     1773.946            56.978     221.952       0.000       0.125         4     212.921

nPrimaryParticle, tracksFromFSL are incorrectly assigned. Anything related to FSL should only work with neutrino interaction generator, not for single particle generator. Needs to find a way to fix this.

[mvicenzi]

Some investigations...

1. I suppose the second issue comes from: https://github.com/WenjieWu-Sci/FLArE/blob/main/src/AnalysisManager.cc#L340-L360

   // find all the tracks originate from the final state lepton, include FSL itself (TID=1)
   tracksFromFSL.insert(1);
   for (auto x : allTracksPTPair) {
   if (tracksFromFSL.find(x.first) != tracksFromFSL.end()) {
     tracksFromFSL.insert(x.second);
   }
   }

   Here FSL is basically defined as the track with TID = 1. This is true for GENIE interactions, but the vertex created with particle gun sets the track id as 1 as well, causing the pion to be considered a FSL. We should add here a proper PDG check on the TID = 1 track and skip if failed.

2. The issue with nPrimaryParticle comes from where it gets updated (https://github.com/WenjieWu-Sci/FLArE/blob/main/src/AnalysisManager.cc#L328), since just above primary_particle_info->Print(); is called only once for the correct primary.

   // update number of primary particles
   // including decay products from tau and pizero
   nPrimaryParticle = countPrimaryParticle;

   In a test I ran, I see that these wrong primaries have PID=3, so looking at the various ways countPrimaryParticle is set (https://github.com/WenjieWu-Sci/FLArE/blob/main/src/AnalysisManager.cc#L574 and onwards) I believe they are probably coming from being considered FSL secondaries of some kind... so they might be a consequence of the wrong FSL track list.

[WenjieWu-Sci]

I think I've found a workaround, validating with some examples. I'll update as soon as the validation looks reasonable.

[mvicenzi]

While producing plots for #14 I noticed a different error related to nPrimaryParticle happening when the starting position of the muon is in the world volume. In this case nPrimaryParticle=0.

PrimaryParticleInformation: PDG code 13
Particle unique ID : 0
MC momentum : (0,0,5104.564987832) MeV
MC vertex : (0,0,30000) mm

number of primary vertices  : 1
Recorded tracks       : 21
Tracks from FSL       : 19
Tracks from FSL (2nd) : 0
number of primary particles : 0 , in which number of particles from fsl : 0
Can't find the primary particle of the hit, something is wrong 13 edep(0.652446) TID(1) PID(0) creator process (PrimaryParticle) position(0.033951,0.053436,34502.102168)
Can't find the primary particle of the hit, something is wrong 13 edep(1.091510) TID(1) PID(0) creator process (PrimaryParticle) position(0.034705,0.053354,34507.102168)
Can't find the primary particle of the hit, something is wrong 13 edep(2.339208) TID(1) PID(0) creator process (PrimaryParticle) position(0.022532,0.215490,34725.000000)
Can't find the primary particle of the hit, something is wrong 13 edep(1.148989) TID(1) PID(0) creator process (PrimaryParticle) position(-0.030680,0.200767,34943.794435)
Can't find the primary particle of the hit, something is wrong 13 edep(0.713474) TID(1) PID(0) creator process (PrimaryParticle) position(-0.026756,0.206044,34948.794435)
...

I will investigate if the fix also solves this... tomorrow :)

[WenjieWu-Sci]

It won't be solved by this fix.

1. If the start position of the particle is not in a sensitive detector, the first step won't be recorded and it will fail this condition to assign a value to countPrimaryParticle. If all events are generated within the detector it should be fine. But we might need to deal with it eventually because we're going to study backgrounds generated outside the detector. One workaround I can think of for now is changing this line to: if (countPrimaryParticle != 0) nPrimaryParticle = countPrimaryParticle; In this way, if the particle is generated outside a SD, nPrimaryParticle would be assigned directly by count_particles from the generator. I'll give it a try.
2. The error with "Can't find the primary particle of the hit ......" is caused by the gap between sensitive detectors. Since not all the steps are recorded in the hitCollection, some of the tracks can't find their parent tracks in the hitCollection and hence can't trace back to theirs primary tracks. I'm not sure if this is a real problem and I don't know how to resolve it without converting all the components to sensitive detectors. (Actually, maybe it can be done with a UserSteppingAction.)

[WenjieWu-Sci]

One workaround I can think of for now is changing this line to: if (countPrimaryParticle != 0) nPrimaryParticle = countPrimaryParticle; In this way, if the particle is generated outside a SD, nPrimaryParticle would be assigned directly by count_particles from the generator. I'll give it a try.

No, this can't resolve this issue. I think the way to do this is via a UserSteppingAction. Either way, we probably should open another issue for this. It's another problem rather than the issue in this thread.