Using the following command, the number of secondaries and their respective species have been written out using mgdraw.f : WRITE(98,*) NP,KPART(1:NP)

In the fort.98 file, I could see this list. It will be helpful if you kindly explain what does the number refer to; for e.g., in the first line, I guess 16 is total number of secondaries, then there will be 16 such no of species i.e. their ID numbers. Is it so ? There are also some negative numbers.

Can you please link the corresponding manual pages where the meaning of each number is listed (i.e. which number refers to what particle?)

As a test case, I kept 10 histories. So for each cycle there will be 10 no of primary projectile proton (1 GeV proton beam is incident on Cu target). In this case I assumed there will be 10 sets of entries in the fort.98 file. However, I got 161 sets in first cycle, 192 sets in the second cycle and so on…
Can you please explain the reason ? Is it because JTRACK = 1 does not account for incident primary proton, rather, it takes into account any proton ??

If the projectile is a proton:
IF (JTRACK.EQ.1) THEN
WRITE(98,*) NP,KPART(1:NP)

So, if I want to store the information for primary proton beam only, then should I write the condition as:

If the projectile is a primary proton:
IF (JTRACK.EQ.1 .AND. (ABS(ETRACK) .EQ. TKI(1))) THEN
WRITE(98,*) NP,KPART(1:NP)

Of course: secondary particles re-interact in turn. One primary history does not necessarily feature one interaction only, as can be read in the output file’s final tables.

If you want to isolate the primary proton interaction, you should rather put the condition:

IF (JTRACK.EQ.1 .AND. LTRACK.EQ.1) THEN

where LTRACK is the generation number of the interacting particle.

What you wrote actually makes no sense, since you are comparing the interacting proton’s TOTAL energy (ETRACK, which is always positive and does not require ABS) with the KINETIC energy of the first product of the reaction (that can be any secondary particle).
Moreover, from a coding point of view, checking if two real numbers are equal is a discouraged practice, due to intrinsic precision limits that may always imply a negative answer.

I do not see here the input file, but I understand from your mgdraw comments that you are biasing nuclear interactions (with LAM-BIAS). If so, the primary proton, when surviving the biased nuclear interaction, may undergo a new one in the same history, hence the data set excess. For an event-by-event analysis as the one you are setting up, I suggest to drop the biasing.