Incorrect initialization of Pythia 8

[gaudu]

I am currently working on an interface between Pythia 8.3.10 and Apprentice. This interface is meant to generate collision events using Pythia 8.3.10 with parameter values for MultipartonInteractions:pT0Ref, StringZ:aLund and StringZ:bLund such as dictated by Apprentice. Doing so will allow the Apprentice framework to determine the best fitted parameter values to finalize the tune of Pythia 8.3.10 against NA61/SHINE experimental dataset.

In a first attempt of this Pythia 8 + Apprentice interface, I created a branch apprentice_tune and the example program main1009.cc alongside few bash scripts enabling the reading of the parameter values from local Apprentice installation for Pythia 8 to properly initialize with the various parameter iterations.

Current status is that everything seems to be correctly implemented and all input values (MultipartonInteractions:pT0Ref, StringZ:aLund, StringZ:bLund and Beams:eB) are read correctly by the Pythia example program; or so I am confident after printing away and checking logfiles:

Beams:eB = 0.931
Check for (i=2; i<5; i++): i = 2, argv[i] = MultipartonInteractions:pT0Ref=2.8184807600134794
Check for (i=2; i<5; i++): i = 3, argv[i] = StringZ:aLund=0.18085760308865206
Check for (i=2; i<5; i++): i = 4, argv[i] = StringZ:bLund=1.3739677616120314

and

*-------  PYTHIA Flag + Mode + Parm + Word + FVec + MVec + PVec + WVec Settings (changes only)  ------------------* 
 |                                                                                                                 | 
 | Name                                          |                      Now |      Default         Min         Max | 
 |                                               |                          |                                      | 
 | Beams:eA                                      |                158.00000 |     7000.000         0.0             | 
 | Beams:eB                                      |                  0.93149 |     7000.000         0.0             | 
 | Beams:frameType                               |                        2 |            1           1           5 | 
 | Beams:idA                                     |                     -211 |         2212                         | 
 | HadronLevel:all                               |                      off |           on                         | 
 | MultipartonInteractions:pT0Ref                |                  2.81848 |      2.28000     0.50000    10.00000 | 
 | Next:numberCount                              |                        0 |         1000           0             | 
 | Next:numberShowEvent                          |                        0 |            1           0             | 
 | Next:numberShowInfo                           |                        0 |            1           0             | 
 | Next:numberShowLHA                            |                        0 |            1           0             | 
 | Next:numberShowProcess                        |                        0 |            1           0             | 
 | SoftQCD:all                                   |                       on |          off                         | 
 | StringZ:aLund                                 |                  0.18086 |      0.68000         0.0     2.00000 | 
 | StringZ:bLund                                 |                  1.37397 |      0.98000     0.20000     2.00000 | 
 | Tune:ee                                       |                        0 |            7          -1           7 | 
 | Tune:pp                                       |                        0 |           14          -1          34 | 
 |                                                                                                                 | 
 *-------  End PYTHIA Flag + Mode + Parm + Word + FVec + MVec + PVec + WVec Settings  -----------------------------* 

I am nonetheless facing an issue as of now, where further in the Pythia 8 "Multiparton Interactions" initialization the program gets stuck in an infinite loop awaiting an acceptable value for sigmaInteraction. This seems to never come since an input value somewhere is -nan (not a number). This -nan appears in the very beginning of the Pythia 8 initialization:

*----------------------  Initializing Angantyr  ----------------------*
 |                    We collide: pi- on 12C                           |
 |                                                                     |
 |                    Below follows initialization                     |
 |                    of sub-collisions.                               |
 |                                                                     |
 |                   //>________________________________               |
 |          [########[]_________________________________>              |
 |                   \\>                                               |
 |                                                                     |
 *------ HeavyIon fitting of SubCollisionModel to cross sections ------* 
 |                                                                     | 
 |   Fitting parameters for     -nan GeV                               | 

and later on during the "Multiparton Interactions" initialization:

 *-------  PYTHIA Multiparton Interactions Initialization  ---------* 
 |                                                                  | 
 |                   sigmaNonDiffractive =    16.06 mb              | 
 |                                                                  | 
 |    pT0 =  -nan gives sigmaInteraction =     -nan mb: rejected    | 
 |    pT0 =  -nan gives sigmaInteraction =     -nan mb: rejected    | 
 |    pT0 =  -nan gives sigmaInteraction =     -nan mb: rejected    | 
 |    pT0 =  -nan gives sigmaInteraction =     -nan mb: rejected    | 
 |    pT0 =  -nan gives sigmaInteraction =     -nan mb: rejected    | 

[gaudu]

I am now trying to figure out what is the variable behind the -nan in Fitting parameters for -nan GeV. This lead me to the HISubCollisionModel.cc file and to the eCM variable in the SubCollisionModel::evolve(int nGenerations, double eCM) function.

From HISubCollisionModel.hh:

// Use a simplified genetic algorithm to fit the parameters.
  virtual bool evolve(int nGenerations, double eCM);

[gaudu]

I tried to play around with the energy settings in main1009.cc since this is linked to the eCM input into a SubCollisionModel function. I therefore only initialized Pythia by setting the Beams:idA and Beams:idB, leaving out the Beams:frameType, Beams:eA and Beams:eB (expecting the energy to be set as the default 14 TeV):

pythia.readString("Beams:idA = -211");
pythia.readString("Beams:idB = 1000060120");

It indeed ran as expected:

*----------------------  Initializing Angantyr  ----------------------*
 |                    We collide: pi- on 12C                           |
 |                                                                     |
 |                    Below follows initialization                     |
 |                    of sub-collisions.                               |
 |                                                                     |
 |                   //>________________________________               |
 |          [########[]_________________________________>              |
 |                   \\>                                               |
 |                                                                     |
 *------ HeavyIon fitting of SubCollisionModel to cross sections ------* 
 |                                                                     | 
 |   Fitting parameters for  14000.0 GeV                               | 

There must be something that Pythia does not appreciate when the previous combination of settings:

pythia.readString("Beams:frameType = 2");
pythia.readString("Beams:eA = 158.");
pythia.readString("Beams:eB = 0.93149");

[gaudu]

Following this track and modifying the values of Beams:eB:	Beams:eB	eCM in SubCollisionModel
0.9	-nan
pythia.particleData.m0(1000060120)/12 ~ 0.931	-nan
pythia.particleData.m0(2212) ~ 0.938	17.2
1	20.6

Now I found a fix, where I can use the proton mass instead of the carbon mass divided by 12. Question remains, why is this happening? Is this a feature with a scientific purpose or a bug?

[gaudu]

(Temporary?) bug fix implemented in https://github.com/gaudu/pythia8310/commit/38c1215ae1a8969cda1ada6dec35b9db3ffffe38.