StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
HeavyIons.cc
1 // HeavyIons.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // Function definitions (not found in the HeavyIons.h header) for the
7 // heavy ion classes classes, and some related global functions.
8 
9 #include "Pythia8/HeavyIons.h"
10 #include <cassert>
11 
12 namespace Pythia8 {
13 
14 //==========================================================================
15 
16 // The abstract HeavyIons class
17 
18 //--------------------------------------------------------------------------
19 
20 // Before doing anything Pythia should add special heavy ion versions
21 // for some groups of settings.
22 
23 void HeavyIons::addSpecialSettings(Settings & settings) {
24  setupSpecials(settings, "Diffraction:");
25  setupSpecials(settings, "MultipartonInteractions:");
26  setupSpecials(settings, "PDF:");
27  setupSpecials(settings, "SigmaDiffractive:");
28  setupSpecials(settings, "BeamRemnants:");
29 }
30 
31 void HeavyIons::setupSpecials(Settings & settings, string match) {
32  map<string,Flag> flags = settings.getFlagMap(match);
33  for ( map<string,Flag>::iterator it = flags.begin();
34  it != flags.end(); ++it )
35  settings.addFlag("HI" + it->second.name, it->second.valDefault);
36  map<string,Mode> modes = settings.getModeMap(match);
37  for ( map<string,Mode>::iterator it = modes.begin();
38  it != modes.end(); ++it )
39  settings.addMode("HI" + it->second.name, it->second.valDefault,
40  it->second.hasMin, it->second.hasMax,
41  it->second.valMin, it->second.valMax, it->second.optOnly);
42  map<string,Parm> parms = settings.getParmMap(match);
43  for ( map<string,Parm>::iterator it = parms.begin();
44  it != parms.end(); ++it )
45  settings.addParm("HI" + it->second.name, it->second.valDefault,
46  it->second.hasMin, it->second.hasMax,
47  it->second.valMin, it->second.valMax);
48  map<string,Word> words = settings.getWordMap(match);
49  for ( map<string,Word>::iterator it = words.begin();
50  it != words.end(); ++it )
51  settings.addWord("HI" + it->second.name, it->second.valDefault);
52  map<string,FVec> fvecs = settings.getFVecMap(match);
53  for ( map<string, FVec>::iterator it = fvecs.begin();
54  it != fvecs.end(); ++it )
55  settings.addFVec("HI" + it->second.name, it->second.valDefault);
56  map<string,MVec> mvecs = settings.getMVecMap(match);
57  for ( map<string,MVec>::iterator it = mvecs.begin();
58  it != mvecs.end(); ++it )
59  settings.addMVec("HI" + it->second.name, it->second.valDefault,
60  it->second.hasMin, it->second.hasMax,
61  it->second.valMin, it->second.valMax);
62  map<string,PVec> pvecs = settings.getPVecMap(match);
63  for ( map<string,PVec>::iterator it = pvecs.begin();
64  it != pvecs.end(); ++it )
65  settings.addPVec("HI" + it->second.name, it->second.valDefault,
66  it->second.hasMin, it->second.hasMax,
67  it->second.valMin, it->second.valMax);
68  map<string,WVec> wvecs = settings.getWVecMap(match);
69  for ( map<string,WVec>::iterator it = wvecs.begin();
70  it != wvecs.end(); ++it )
71  settings.addWVec("HI" + it->second.name, it->second.valDefault);
72 }
73 
74 void HeavyIons::setupSpecials(Pythia & p, string match) {
75  Settings & opts = p.settings;
76  map<string, Flag> flags = opts.getFlagMap(match);
77  for ( map<string, Flag>::iterator it = flags.begin();
78  it != flags.end(); ++it )
79  opts.flag(it->second.name.substr(2), it->second.valNow, true);
80  map<string, Mode> modes = opts.getModeMap(match);
81  for ( map<string, Mode>::iterator it = modes.begin();
82  it != modes.end(); ++it )
83  opts.mode(it->second.name.substr(2), it->second.valNow, true);
84  map<string, Parm> parms = opts.getParmMap(match);
85  for ( map<string, Parm>::iterator it = parms.begin();
86  it != parms.end(); ++it )
87  opts.parm(it->second.name.substr(2), it->second.valNow, true);
88  map<string, Word> words = opts.getWordMap(match);
89  for ( map<string, Word>::iterator it = words.begin();
90  it != words.end(); ++it )
91  opts.word(it->second.name.substr(2), it->second.valNow, true);
92  map<string, FVec> fvecs = opts.getFVecMap(match);
93  for ( map<string, FVec>::iterator it = fvecs.begin();
94  it != fvecs.end(); ++it )
95  opts.fvec(it->second.name.substr(2), it->second.valNow, true);
96  map<string, MVec> mvecs = opts.getMVecMap(match);
97  for ( map<string, MVec>::iterator it = mvecs.begin();
98  it != mvecs.end(); ++it )
99  opts.mvec(it->second.name.substr(2), it->second.valNow, true);
100  map<string, PVec> pvecs = opts.getPVecMap(match);
101  for ( map<string, PVec>::iterator it = pvecs.begin();
102  it != pvecs.end(); ++it )
103  opts.pvec(it->second.name.substr(2), it->second.valNow, true);
104  map<string, WVec> wvecs = opts.getWVecMap(match);
105  for ( map<string, WVec>::iterator it = wvecs.begin();
106  it != wvecs.end(); ++it )
107  opts.wvec(it->second.name.substr(2), it->second.valNow, true);
108 }
109 
110 //--------------------------------------------------------------------------
111 
112 // Sum up info from all used Pythia objects.
113 
114 void HeavyIons::sumUpMessages(Info & in, string tag, const Info * ihi) {
115  if ( !ihi ) return;
116  for ( auto it : ihi->messages ) in.messages[tag + it.first] += it.second;
117 }
118 
119 //--------------------------------------------------------------------------
120 
121 // Reset all process level settings in the given Pythia object. NOTE
122 // must be expanded if new process groups are included in Pythia.
123 
124 void HeavyIons::clearProcessLevel(Pythia & pyt) {
125  string path = pyt.settings.word("xmlPath");
126  pyt.settings.mode("Tune:ee", 0);
127  pyt.settings.mode("Tune:pp", 0);
128  pyt.settings.init(path + "QCDSoftProcesses.xml", true);
129  pyt.settings.init(path + "QCDHardProcesses.xml", true);
130  pyt.settings.init(path + "ElectroweakProcesses.xml", true);
131  pyt.settings.init(path + "OniaProcesses.xml", true);
132  pyt.settings.init(path + "TopProcesses.xml", true);
133  pyt.settings.init(path + "FourthGenerationProcesses.xml", true);
134  pyt.settings.init(path + "HiggsProcesses.xml", true);
135  pyt.settings.init(path + "SUSYProcesses.xml", true);
136  pyt.settings.init(path + "NewGaugeBosonProcesses.xml", true);
137  pyt.settings.init(path + "LeftRightSymmetryProcesses.xml", true);
138  pyt.settings.init(path + "LeptoquarkProcesses.xml", true);
139  pyt.settings.init(path + "CompositenessProcesses.xml", true);
140  pyt.settings.init(path + "HiddenValleyProcesses.xml", true);
141  pyt.settings.init(path + "ExtraDimensionalProcesses.xml", true);
142  pyt.settings.init(path + "DarkMatterProcesses.xml", true);
143  pyt.settings.init(path + "SecondHardProcess.xml", true);
144  pyt.settings.init(path + "PhaseSpaceCuts.xml", true);
145  // NOTE! if new processes are added in separate xml files these have
146  // to be added here.
147 }
148 
149 //--------------------------------------------------------------------------
150 
151 // Update the Info object in the main Pythia object.
152 
153 void HeavyIons::updateInfo() {
154  map<string, int> saveMess = infoPtr->messages;
155  *infoPtr = hiInfo.primInfo;
156  infoPtr->hiInfo = &hiInfo;
157  infoPtr->messages = saveMess;
158  infoPtr->weightContainerPtr->setWeightNominal(hiInfo.weight());
159  infoPtr->sigmaReset();
160  double norm = 1.0/double(hiInfo.NSave);
161  int Nall = 0;
162  double wall = 0.0;
163  double w2all = 0.0;
164  for ( map<int,int>::iterator ip = hiInfo.NPrim.begin();
165  ip != hiInfo.NPrim.end(); ++ip ) {
166  int N = ip->second;
167  if ( !N ) continue;
168  int pc = ip->first;
169  double w = hiInfo.sumPrimW[pc]/millibarn;
170  double w2 = hiInfo.sumPrimW2[pc]/pow2(millibarn);
171  infoPtr->setSigma(pc, hiInfo.NamePrim[pc], N, N, N,
172  w*norm, sqrt(w2*norm)/N, w);
173  Nall += N;
174  wall += w;
175  w2all += w2;
176  }
177  infoPtr->setSigma(0, "sum", hiInfo.NSave, Nall, Nall,
178  wall*norm, sqrt(w2all*norm)/Nall, wall);
179 }
180 
181 //--------------------------------------------------------------------------
182 
183 // Print out statistics from a HeavyIons run.
184 
185 void HeavyIons::stat() {
186  bool showPrL = flag("Stat:showProcessLevel");
187  // bool showPaL = settings.flag("Stat:showPartonLevel");
188  bool showErr = flag("Stat:showErrors");
189  bool reset = flag("Stat:reset");
190  Info & in = *infoPtr;
191  // Header.
192  if ( showPrL ) {
193  cout << "\n *----- HeavyIon Event and Cross Section Statistics ------"
194  << "-------------------------------------------------------*\n"
195  << " | "
196  << " |\n"
197  << " | Subprocess Code | "
198  << " Number of events | sigma +- delta |\n"
199  << " | | "
200  << "Tried Selected Accepted | (estimated) (mb) |\n"
201  << " | | "
202  << " | |\n"
203  << " |------------------------------------------------------------"
204  << "-----------------------------------------------------|\n"
205  << " | | "
206  << " | |\n";
207 
208  vector<int> pc = in.codesHard();
209  for ( int i = 0, N = pc.size(); i < N; ++i ) {
210  cout << " | " << left << setw(45) << in.nameProc(pc[i])
211  << right << setw(5) << pc[i] << " | "
212  << setw(11) << in.nTried(pc[i]) << " "
213  << setw(10) << in.nSelected(pc[i]) << " "
214  << setw(10) << in.nAccepted(pc[i]) << " | "
215  << scientific << setprecision(3)
216  << setw(11) << in.sigmaGen(pc[i])
217  << setw(11) << in.sigmaErr(pc[i]) << " |\n";
218  }
219  if ( pc.empty() ) in.setSigma(0, "sum", hiInfo.NSave, 0, 0, 0.0, 0.0, 0.0);
220 
221  cout << " | | "
222  << " | |\n"
223  << " | " << left << setw(50) << "sum" << right << " | " << setw(11)
224  << in.nTried(0) << " " << setw(10) << in.nSelected(0) << " "
225  << setw(10) << in.nAccepted(0) << " | " << scientific
226  << setprecision(3) << setw(11)
227  << in.sigmaGen(0) << setw(11) << in.sigmaErr(0) << " |\n";
228  cout << " | " << left << setw(50) << "(Estimated total cross section)"
229  << right << " | " << setw(11)
230  << hiInfo.nAttempts() << " " << setw(10) << 0 << " " << setw(10)
231  << 0 << " | " << scientific << setprecision(3) << setw(11)
232  << hiInfo.sigmaTot() << setw(11) << hiInfo.sigmaTotErr() << " |\n";
233  cout << " | " << left << setw(50)
234  << "(Estimated non-diffractive cross section)"
235  << right << " | " << setw(11)
236  << hiInfo.nAttempts() << " " << setw(10) << 0 << " " << setw(10)
237  << 0 << " | " << scientific << setprecision(3) << setw(11)
238  << hiInfo.sigmaND() << setw(11) << hiInfo.sigmaNDErr() << " |\n";
239  // Listing finished.
240  cout << " | "
241  << " |\n"
242  << " *----- End HeavyIon Event and Cross Section Statistics -----"
243  << "-----------------------------------------------------*" << endl;
244  }
245  if ( reset ) hiInfo = HIInfo();
246  if ( showErr ) {
247  for ( int i = 1, np = pythia.size(); i < np; ++i )
248  sumUpMessages(in, "(" + pythiaNames[i] + ")", info[i]);
249  in.errorStatistics();
250  }
251  if ( reset ) in.errorReset();
252 
253 }
254 
255 //--------------------------------------------------------------------------
256 
257 // Check the settings and return false of there are no heavy ion beams.
258 
259 bool HeavyIons::isHeavyIon(Settings & settings) {
260  int idProj = settings.mode("Beams:idA");
261  int idTarg = settings.mode("Beams:idB");
262  return ( abs(idProj/100000000) == 10 ||abs(idTarg/100000000) == 10 );
263 }
264 
265 //==========================================================================
266 
267 // Angantyr is the main HeavyIons model in Pythia.
268 
269 //--------------------------------------------------------------------------
270 
271 // Constructor.
272 
273 Angantyr::Angantyr(Pythia & mainPythiaIn)
274  : HeavyIons(mainPythiaIn), hasSignal(true),
275  bGenPtr(0), projPtr(0), targPtr(0), collPtr(0), recoilerMode(1), bMode(0) {
276  selectMB = make_shared<ProcessSelectorHook>();
277  selectSASD = make_shared<ProcessSelectorHook>();
278  pythia.resize(ALL);
279  info.resize(ALL);
280  pythiaNames.resize(ALL);
281  pythiaNames[HADRON] = "HADRON";
282  pythiaNames[MBIAS] = "MBIAS";
283  pythiaNames[SASD] = "SASD";
284  pythiaNames[SIGPP] = "SIGPP";
285  pythiaNames[SIGPN] = "SIGPN";
286  pythiaNames[SIGNP] = "SIGNP";
287  pythiaNames[SIGNN] = "SIGNN";
288 
289 }
290 
291 //--------------------------------------------------------------------------
292 
293 // Destructor deleting model objects that are not provided from the
294 // outside (via HIUserHooks).
295 
296 Angantyr::~Angantyr() {
297  for ( int i = MBIAS; i < ALL; ++i ) if ( pythia[i] ) delete pythia[i];
298  if ( !( HIHooksPtr && HIHooksPtr->hasImpactParameterGenerator() ) )
299  delete bGenPtr;
300  if ( !( HIHooksPtr && HIHooksPtr->hasProjectileModel() ) )
301  delete projPtr;
302  if ( !( HIHooksPtr && HIHooksPtr->hasTargetModel() ) )
303  delete targPtr;
304  if ( !( HIHooksPtr && HIHooksPtr->hasSubCollisionModel() ) )
305  delete collPtr;
306 }
307 
308 //--------------------------------------------------------------------------
309 
310 // Add a HIUserHooks object to customise the Angantyr model.
311 
312 bool Angantyr::setUserHooksPtr(PythiaObject sel, shared_ptr<UserHooks> uhook) {
313  for ( int i = HADRON; i < ALL; ++i )
314  if ( ( i == sel || ALL == sel ) && !pythia[i]->setUserHooksPtr(uhook) )
315  return false;
316  return true;
317 }
318 
319 //--------------------------------------------------------------------------
320 
321 // Create an EventInfo object connected to a SubCollision from the
322 // last event generated by the given PythiaObject.
323 
324 EventInfo Angantyr::mkEventInfo(Pythia & pyt, Info & infoIn,
325  const SubCollision * coll) {
326  EventInfo ei;
327  ei.coll = coll;
328  ei.event = pyt.event;
329  ei.info = infoIn;
330  ei.code = pyt.info.code();
331  ei.ordering = ( ( HIHooksPtr && HIHooksPtr->hasEventOrdering() )?
332  HIHooksPtr->eventOrdering(ei.event, infoIn):
333  pyt.info.bMPI() );
334  if ( coll ) {
335  ei.projs[coll->proj] = make_pair(1, ei.event.size());
336  ei.targs[coll->targ] = make_pair(2, ei.event.size());
337  }
338 
339  ei.ok = true;
340  return ei;
341  }
342 
343 //--------------------------------------------------------------------------
344 
345 // Initialise Angantyr. Called from within Pythia::init().
346 
347 bool Angantyr::init() {
348 
349  bool print = flag("HeavyIon:showInit");
350 
351  int idProj = mode("Beams:idA");
352  int idTarg = mode("Beams:idB");
353  int idProjP = idProj;
354  int idProjN = 0;
355  int idTargP = idTarg;
356  int idTargN = 0;
357  bool isHIProj = ( abs(idProj/100000000) == 10 );
358  bool isHITarg = ( abs(idTarg/100000000) == 10 );
359  bool isHI = isHIProj || isHITarg || mode("HeavyIon:mode") > 1;
360  if ( isHIProj ) {
361  idProjP = idProj > 0? 2212: -2212;
362  idProjN = idProj > 0? 2112: -2112;
363  }
364  if ( abs(idTarg/100000000) == 10 ) {
365  idTargP = idTarg > 0? 2212: -2212;
366  idTargN = idTarg > 0? 2112: -2112;
367  }
368  if ( mode("HeavyIon:mode") == 1 && !isHI ) {
369  infoPtr->errorMsg("Angantyr Info: No heavy ions requested"
370  " - reverting to normal Pythia behavior.");
371  settingsPtr->mode("HeavyIon:mode", 0);
372  return false;
373  }
374 
375  recoilerMode = mode("Angantyr:SDRecoil");
376  bMode = mode("Angantyr:impactMode");
377 
378  int frame = mode("Beams:frameType");
379  bool dohad = flag("HadronLevel:all");
380  if ( frame > 2 )
381  infoPtr->errorMsg("Angantyr warning: Currently only Beams:frameType "
382  "= 1 or 2 is supported. Assuming 2.");
383  double eAbm = parm("Beams:eA");
384  double eBbm = parm("Beams:eB");
385  if ( frame == 1 ) eAbm = eBbm = parm("Beams:eCM")/2.0;
386  settingsPtr->parm("Beams:eA", eAbm);
387  settingsPtr->parm("Beams:eB", eBbm);
388  settingsPtr->mode("Beams:frameType", 2);
389 
390  settingsPtr->mode("Next:numberCount", 0);
391  settingsPtr->mode("Next:numberShowLHA", 0);
392  settingsPtr->mode("Next:numberShowInfo", 0);
393  settingsPtr->mode("Next:numberShowProcess", 0);
394  settingsPtr->mode("Next:numberShowEvent", 0);
395  settingsPtr->flag("HadronLevel:all", false);
396  settingsPtr->flag("SoftQCD:all", false);
397  settingsPtr->flag("SoftQCD:elastic", false);
398  settingsPtr->flag("SoftQCD:nonDiffractive", false);
399  settingsPtr->flag("SoftQCD:singleDiffractive", false);
400  settingsPtr->flag("SoftQCD:doubleDiffractive", false);
401  settingsPtr->flag("SoftQCD:centralDiffractive", false);
402 
403  for ( int i = MBIAS; i < ALL; ++i ) {
404  pythia[i] = new Pythia(*settingsPtr, *particleDataPtr, false);
405  pythia[i]->settings.mode("HeavyIon:mode", 1);
406  }
407 
408  sigtot.init();
409  sigtot.calc(2212, 2212, sqrt(4.0*eAbm*eBbm));
410 
411  clearProcessLevel(*pythia[MBIAS]);
412  pythia[MBIAS]->settings.flag("SoftQCD:all", true);
413  pythia[MBIAS]->settings.mode("Beams:idA", idProjP);
414  pythia[MBIAS]->settings.mode("Beams:idB", idTargP);
415 
416  clearProcessLevel(*pythia[SASD]);
417  Settings & sdabsopts = pythia[SASD]->settings;
418  sdabsopts.flag("SoftQCD:singleDiffractive", true);
419 
420  setupSpecials(*pythia[SASD], "HIDiffraction:");
421  setupSpecials(*pythia[SASD], "HIMultipartonInteractions:");
422  setupSpecials(*pythia[SASD], "HIPDF:");
423  setupSpecials(*pythia[SASD], "HISigmaDiffractive:");
424  setupSpecials(*pythia[SASD], "HIBeamRemnants:");
425  if ( sdabsopts.mode("Angantyr:SASDmode") > 0 ) {
426  double pT0Ref = sdabsopts.parm("MultipartonInteractions:pT0Ref");
427  double ecmRef = sdabsopts.parm("MultipartonInteractions:ecmRef");
428  double ecmPow = sdabsopts.parm("MultipartonInteractions:ecmPow");
429  double ecm = sqrt(4.0*eBbm*eAbm);
430  sdabsopts.parm("Beams:eCM", ecm);
431  double pT0 = pT0Ref * pow(ecm / ecmRef, ecmPow);
432  sdabsopts.parm("MultipartonInteractions:pT0Ref", pT0);
433  sdabsopts.parm("MultipartonInteractions:ecmRef", ecm);
434  sdabsopts.parm("MultipartonInteractions:ecmPow", 0.0);
435  sdabsopts.word("PDF:PomSet", "11");
436  if ( sdabsopts.mode("Angantyr:SASDmode") == 2 ) {
437  sdabsopts.parm("Diffraction:mRefPomP", ecm);
438  double sigND = sigtot.sigmaND();
439  double mmin = sdabsopts.parm("Diffraction:mMinPert");
440  double powp = sdabsopts.parm("HIDiffraction:mPowPomP");
441  sdabsopts.parm("Diffraction:mPowPomP", powp, true);
442  if ( powp > 0.0 ) sigND /= ((1.0 - pow(mmin/ecm, powp))/powp);
443  else sigND /= log(ecm/mmin);
444  sdabsopts.parm("Diffraction:sigmaRefPomP", sigND, true);
445  }
446  if ( sdabsopts.mode("Angantyr:SASDmode") >= 3 ) {
447  sdabsopts.parm("Diffraction:mRefPomP", ecm);
448  double sigND = sigtot.sigmaND();
449  sdabsopts.parm("Diffraction:sigmaRefPomP", sigND, true);
450  sdabsopts.parm("Diffraction:mPowPomP", 0.0);
451  }
452  }
453  sdabsopts.mode("Beams:idA", idProjP);
454  sdabsopts.mode("Beams:idB", idTargP);
455 
456  clearProcessLevel(*pythia[HADRON]);
457  pythia[HADRON]->settings.flag("ProcessLevel:all", false);
458  pythia[HADRON]->settings.flag("PartonLevel:all", false);
459  pythia[HADRON]->settings.flag("HadronLevel:all", dohad);
460  pythia[HADRON]->settings.mode("Beams:idA", idProj);
461  pythia[HADRON]->settings.mode("Beams:idB", idTarg);
462 
463  pythia[SIGPP]->settings.mode("Beams:idA", idProjP);
464  pythia[SIGPP]->settings.mode("Beams:idB", idTargP);
465  if ( idTargN ) {
466  pythia[SIGPN]->settings.mode("Beams:idA", idProjP);
467  pythia[SIGPN]->settings.mode("Beams:idB", idTargN);
468  }
469  if ( idProjN ) {
470  pythia[SIGNP]->settings.mode("Beams:idA", idProjN);
471  pythia[SIGNP]->settings.mode("Beams:idB", idTargP);
472  }
473  if ( idProjN && idTargN ) {
474  pythia[SIGNN]->settings.mode("Beams:idA", idProjN);
475  pythia[SIGNN]->settings.mode("Beams:idB", idTargN);
476  }
477 
478  if ( HIHooksPtr ) HIHooksPtr->init(idProj, idTarg);
479 
480  if ( HIHooksPtr && HIHooksPtr->hasProjectileModel() )
481  projPtr = HIHooksPtr->projectileModel();
482  else
483  projPtr = new GLISSANDOModel();
484  projPtr->initPtr(idProj, *settingsPtr, *particleDataPtr, *rndmPtr);
485 
486  if ( HIHooksPtr && HIHooksPtr->hasTargetModel() )
487  targPtr = HIHooksPtr->targetModel();
488  else
489  targPtr = new GLISSANDOModel();
490  targPtr->initPtr(idTarg, *settingsPtr, *particleDataPtr, *rndmPtr);
491 
492  if ( HIHooksPtr && HIHooksPtr->hasSubCollisionModel() )
493  collPtr = HIHooksPtr->subCollisionModel();
494  else if ( mode("Angantyr:CollisionModel") == 1 )
495  collPtr = new DoubleStrikman();
496  else if ( mode("Angantyr:CollisionModel") == 2 )
497  collPtr = new DoubleStrikman(1);
498  else if ( mode("Angantyr:CollisionModel") == 3 )
499  collPtr = new BlackSubCollisionModel();
500  else
501  collPtr = new NaiveSubCollisionModel();
502 
503  collPtr->initPtr(*projPtr, *targPtr, sigtot, *settingsPtr,
504  *infoPtr, *rndmPtr);
505  if ( !collPtr->init() ) return false;
506 
507  if ( HIHooksPtr && HIHooksPtr->hasImpactParameterGenerator() )
508  bGenPtr = HIHooksPtr->impactParameterGenerator();
509  else
510  bGenPtr = new ImpactParameterGenerator();
511  bGenPtr->initPtr(*collPtr, *projPtr, *targPtr, *settingsPtr, *rndmPtr);
512 
513  if ( !projPtr->init() ) return false;
514  if ( !targPtr->init() ) return false;
515  if ( !bGenPtr->init() ) return false;
516 
517  string output;
518  if ( hasSignal ) {
519  ostringstream oss;
520  Redirect red(cout, oss);
521  hasSignal = init(SIGPP, "signal process (pp)", 10);
522  output = oss.str();
523  }
524  if ( !hasSignal ) {
525  if ( print ) cout << " Angantyr Info: No signal process specified. "
526  << "Assuming minimum bias." << endl;
527  } else {
528  if ( print )
529  cout << output;
530  if ( idTargN ) init(SIGPN, "signal process (pn)", 10);
531  if ( idProjN ) init(SIGNP, "signal process (np)", 10);
532  if ( idProjN && idTargN ) init(SIGNN, "signal process (nn)", 10);
533 
534  }
535 
536  pythia[MBIAS]->addUserHooksPtr(selectMB);
537  init(MBIAS, "minimum bias processes");
538 
539  pythia[SASD]->addUserHooksPtr(selectSASD);
540  init(SASD, "secondary absorptive processes as single diffraction.");
541 
542  if ( pythia[HADRON]->flag("HadronLevel:all") ) {
543  if ( print )
544  cout << " Angantyr Info: Initializing hadronisation processes." << endl;
545  }
546  settingsPtr->flag("ProcessLevel:all", false);
547 
548  return true;
549 
550 }
551 
552 //--------------------------------------------------------------------------
553 
554 // Initiaize a specific Pythia object and optionally run a number
555 // of events to get a handle of the cross section.
556 
557 bool Angantyr::init(PythiaObject sel, string name, int n) {
558  bool print = flag("HeavyIon:showInit");
559  shared_ptr<InfoGrabber> ihg = make_shared<InfoGrabber>();
560  pythia[sel]->addUserHooksPtr(ihg);
561  if ( print ) cout << " Angantyr Info: Initializing " << name << "." << endl;
562  if ( !pythia[sel]->init() ) return false;
563  info[sel] = ihg->getInfo();
564  if ( n <= 0 ) return true;
565  if ( print ) cout << "Generating a few signal events for " << name
566  << " to build up statistics" << endl;
567  for ( int i = 0; i < 10; ++i ) pythia[sel]->next();
568  return true;
569 }
570 
571 
572 //--------------------------------------------------------------------------
573 
574 // Generate events and return EventInfo objects for different process
575 // types.
576 
577 EventInfo Angantyr::getSignal(const SubCollision & coll) {
578  if ( !hasSignal ) return EventInfo();
579  int pytsel = SIGPP + coll.nucleons();
580  int itry = MAXTRY;
581  while ( itry-- ) {
582  if ( pythia[pytsel]->next() )
583  return mkEventInfo(*pythia[pytsel], *info[pytsel], &coll);
584  }
585  infoPtr->errorMsg("Warning from PyHIa::next: "
586  "Could not setup signal sub collision.");
587  return EventInfo();
588 }
589 
590 EventInfo Angantyr::getMBIAS(const SubCollision * coll, int procid) {
591  int itry = MAXTRY;
592  double bp = -1.0;
593  if ( bMode > 0 && procid == 101 ) bp = coll->bp;
594  HoldProcess hold(selectMB, procid, bp);
595  while ( --itry ) {
596  if ( !pythia[MBIAS]->next() ) continue;
597  assert( pythia[MBIAS]->info.code() == procid );
598  return mkEventInfo(*pythia[MBIAS], *info[MBIAS], coll);
599  }
600  return EventInfo();
601 }
602 
603 EventInfo Angantyr::getSASD(const SubCollision * coll, int procid) {
604  int itry = MAXTRY;
605  double bp = -1.0;
606  if ( bMode > 1 ) bp = coll->bp;
607  HoldProcess hold(selectSASD, procid, bp);
608  while ( --itry ) {
609  if ( !pythia[SASD]->next() ) continue;
610  assert( pythia[SASD]->info.code() == procid );
611  return mkEventInfo(*pythia[SASD], *info[SASD], coll);
612  }
613  return EventInfo();
614 }
615 
616 //--------------------------------------------------------------------------
617 
618 // Generate primary absorptive (non-diffractive) nucleon-nucleon
619 // sub-collisions.
620 
621 bool Angantyr::genAbs(const multiset<SubCollision> & coll,
622  list<EventInfo> & subevents) {
623  // The fully absorptive
624  vector<multiset<SubCollision>::const_iterator> abscoll;
625  // The partly absorptive
626  vector<multiset<SubCollision>::const_iterator> abspart;
627  // The non-diffractive and signal events
628  multiset<EventInfo> ndeve, sigeve;
629 
630  // Select the primary absorptive sub collisions.
631  for ( multiset<SubCollision>::iterator cit = coll.begin();
632  cit != coll.end(); ++cit ) {
633  if ( cit->type != SubCollision::ABS ) continue;
634  if (!cit->proj->done() && !cit->targ->done() ) {
635  abscoll.push_back(cit);
636  if ( bMode > 0 ) {
637  EventInfo ie = getND(*cit);
638  assert( ie.code == 101 );
639  ndeve.insert(ie);
640  }
641  cit->proj->select();
642  cit->targ->select();
643  } else
644  abspart.push_back(cit);
645  }
646 
647  if ( abscoll.empty() ) return true;
648 
649  int Nabs = abscoll.size();
650  int Nadd = abspart.size();
651 
652  if ( bMode == 0 ) {
653  for ( int i = 0; i < Nabs + Nadd; ++i ) {
654  EventInfo ie = getND();
655  assert( ie.code == 101 );
656  ndeve.insert(ie);
657  }
658  }
659  vector<int> Nii(4, 0);
660  vector<double> w(4, 0.0);
661  double wsum = 0.0;
662  double P1 = 1.0;
663  if ( hasSignal ) {
664 
665  // Count how many potential absorpitve collisions there are for
666  // each iso-spin combination.
667  for ( int i = 0, N = abscoll.size(); i < N; ++i )
668  ++Nii[abscoll[i]->nucleons()];
669  for ( int i = 0, N = abspart.size(); i < N; ++i )
670  ++Nii[abspart[i]->nucleons()];
671 
672  if ( Nii[0] )
673  w[0] = pythia[SIGPP]->info.sigmaGen()*millibarn/collPtr->sigND();
674  if ( Nii[1] )
675  w[1] = pythia[SIGPN]->info.sigmaGen()*millibarn/collPtr->sigND();
676  if ( Nii[2] )
677  w[2] = pythia[SIGNP]->info.sigmaGen()*millibarn/collPtr->sigND();
678  if ( Nii[3] )
679  w[3] = pythia[SIGNN]->info.sigmaGen()*millibarn/collPtr->sigND();
680 
681  wsum = Nii[0]*w[0] + Nii[1]*w[1] + Nii[2]*w[2] + Nii[3]*w[3];
682  P1 = 1.0 - pow(1.0 - w[0], Nii[0])*pow(1.0 - w[1], Nii[1])*
683  pow(1.0 - w[2], Nii[2])*pow(1.0 - w[3], Nii[3]);
684 
685  }
686 
687  bool noSignal = hasSignal;
688 
689  // *** THINK *** Is it ok to always pair the hardest events with the
690  // *** most central sub-collisions, or will this introduce a strange
691  // *** bias?
692  multiset<EventInfo>::iterator it = ndeve.begin();
693  EventInfo ei;
694  for ( int i = 0, N = abscoll.size(); i < N; ++i ) {
695  int b = abscoll[i]->nucleons();
696  if ( Nii[b]
697  && ( noSignal || w[b]*(wsum/P1 - 1.0)/(wsum - w[b]) > rndmPtr->flat())
698  && (ei = getSignal(*abscoll[i])).ok ) {
699  noSignal = false;
700  }
701  else
702  ei =*it++;
703  subevents.push_back(ei);
704  if ( !setupFullCollision(subevents.back(), *abscoll[i],
705  Nucleon::ABS, Nucleon::ABS) )
706  return false;
707  }
708 
709  if ( noSignal ) return false;
710 
711  hiInfo.reweight(P1);
712 
713  return true;
714 
715 }
716 
717 //--------------------------------------------------------------------------
718 
719 // Add secondary absorptive sub-collisions to the primary ones.
720 
721 void Angantyr::addSASD(const multiset<SubCollision> & coll) {
722  // Collect absorptively wounded nucleons in secondary
723  // sub-collisions.
724  int ntry = mode("Angantyr:SDTries");
725  if ( settingsPtr->isMode("HI:SDTries") )
726  ntry = mode("HI:SDTries");
727  for ( multiset<SubCollision>::iterator cit = coll.begin();
728  cit != coll.end(); ++cit )
729  if ( cit->type == SubCollision::ABS ) {
730  if ( cit->targ->done() && !cit->proj->done() ) {
731  EventInfo * evp = cit->targ->event();
732  for ( int itry = 0; itry < ntry; ++itry ) {
733  EventInfo add = getSDabsP(*cit);
734  if ( addNucleonExcitation(*evp, add, true) ) {
735  cit->proj->select(*evp, Nucleon::ABS);
736  break;
737  }
738  if ( itry == ntry - 1 ) hiInfo.failedExcitation();
739  }
740  } else if ( cit->proj->done() && !cit->targ->done() ) {
741  EventInfo * evp = cit->proj->event();
742  for ( int itry = 0; itry < ntry; ++itry ) {
743  EventInfo add = getSDabsT(*cit);
744  if ( addNucleonExcitation(*evp, add, true) ) {
745  cit->targ->select(*evp, Nucleon::ABS);
746  break;
747  }
748  if ( itry == ntry - 1 ) hiInfo.failedExcitation();
749  }
750  }
751  }
752 }
753 
754 //--------------------------------------------------------------------------
755 
756 // Add primary double diffraction sub-collisions.
757 
758 bool Angantyr::addDD(const multiset<SubCollision> & coll,
759  list<EventInfo> & subevents) {
760  // Collect full double diffraction collisions.
761  for ( multiset<SubCollision>::iterator cit = coll.begin();
762  cit != coll.end(); ++cit )
763  if ( cit->type == SubCollision::DDE &&
764  !cit->proj->done() && !cit->targ->done() ) {
765  subevents.push_back(getDD(*cit));
766  if ( !setupFullCollision(subevents.back(), *cit,
767  Nucleon::DIFF, Nucleon::DIFF) )
768  return false;
769  }
770  return true;
771 }
772 
773 //--------------------------------------------------------------------------
774 
775 // Add primary single diffraction sub-collisions.
776 
777 bool Angantyr::addSD(const multiset<SubCollision> & coll,
778  list<EventInfo> & subevents) {
779  // Collect full single diffraction collisions.
780  for ( multiset<SubCollision>::iterator cit = coll.begin();
781  cit != coll.end(); ++cit )
782  if ( !cit->proj->done() && !cit->targ->done() ) {
783  if ( cit->type == SubCollision::SDEP ) {
784  subevents.push_back(getSDP(*cit));
785  if ( !setupFullCollision(subevents.back(), *cit,
786  Nucleon::DIFF, Nucleon::ELASTIC) )
787  return false;
788  }
789  if ( cit->type == SubCollision::SDET ) {
790  subevents.push_back(getSDT(*cit));
791  if ( !setupFullCollision(subevents.back(), *cit,
792  Nucleon::ELASTIC, Nucleon::DIFF) )
793  return false;
794  }
795  }
796  return true;
797 }
798 
799 //--------------------------------------------------------------------------
800 
801 // Add all secondary single diffractive sub-collisions to primary
802 // ones.
803 
804 void Angantyr::addSDsecond(const multiset<SubCollision> & coll) {
805  // Collect secondary single diffractive sub-collisions.
806  int ntry = mode("Angantyr:SDTries");
807  if ( settingsPtr->isMode("HI:SDTries") ) ntry = mode("HI:SDTries");
808  for ( multiset<SubCollision>::iterator cit = coll.begin();
809  cit != coll.end(); ++cit ) {
810  if ( !cit->proj->done() &&
811  ( cit->type == SubCollision::SDEP ||
812  cit->type == SubCollision::DDE ) ) {
813  EventInfo * evp = cit->targ->event();
814  for ( int itry = 0; itry < ntry; ++itry ) {
815  EventInfo add = getSDP(*cit);
816  if ( addNucleonExcitation(*evp, add, false) ) {
817  cit->proj->select(*evp, Nucleon::DIFF);
818  break;
819  }
820  if ( itry == ntry - 1 ) hiInfo.failedExcitation();
821  }
822  }
823  if ( !cit->targ->done() &&
824  ( cit->type == SubCollision::SDET ||
825  cit->type == SubCollision::DDE ) ) {
826  EventInfo * evp = cit->proj->event();
827  for ( int itry = 0; itry < ntry; ++itry ) {
828  EventInfo add = getSDT(*cit);
829  if ( addNucleonExcitation(*evp, add, false) ) {
830  cit->targ->select(*evp, Nucleon::DIFF);
831  break;
832  }
833  if ( itry == ntry - 1 ) hiInfo.failedExcitation();
834  }
835  }
836  }
837 }
838 
839 //--------------------------------------------------------------------------
840 
841 // Add all primary central diffraction sub-colliions
842 
843 bool Angantyr::addCD(const multiset<SubCollision> & coll,
844  list<EventInfo> & subevents) {
845  // Collect full central diffraction collisions.
846  for ( multiset<SubCollision>::iterator cit = coll.begin();
847  cit != coll.end(); ++cit )
848  if ( cit->type == SubCollision::CDE &&
849  !cit->proj->done() && !cit->targ->done() ) {
850  subevents.push_back(getCD(*cit));
851  if ( !setupFullCollision(subevents.back(), *cit,
852  Nucleon::ELASTIC, Nucleon::ELASTIC) )
853  return false;
854  }
855  return true;
856 }
857 
858 //--------------------------------------------------------------------------
859 
860 // Add all secondary central diffraction sub-colliions to primary
861 // ones.
862 
863 void Angantyr::addCDsecond(const multiset<SubCollision> & coll) {
864  // Collect secondary central diffractive sub-collisions.
865  for ( multiset<SubCollision>::iterator cit = coll.begin();
866  cit != coll.end(); ++cit ) {
867  if ( !cit->proj->done() && cit->type == SubCollision::CDE ) {
868  EventInfo * evp = cit->targ->event();
869  EventInfo add = getCD(*cit);
870  if ( addNucleonExcitation(*evp, add, false) ) {
871  cit->proj->select(*evp, Nucleon::ELASTIC);
872  }
873  }
874  if ( !cit->targ->done() && cit->type == SubCollision::CDE ) {
875  EventInfo * evp = cit->proj->event();
876  EventInfo add = getCD(*cit);
877  if ( addNucleonExcitation(*evp, add, false) ) {
878  cit->targ->select(*evp, Nucleon::ELASTIC);
879  }
880  }
881  }
882 }
883 
884 //--------------------------------------------------------------------------
885 
886 // Add all primary elastic sub-colliions
887 
888 bool Angantyr::addEL(const multiset<SubCollision> & coll,
889  list<EventInfo> & subevents) {
890  // Collect full elastic collisions.
891  for ( multiset<SubCollision>::iterator cit = coll.begin();
892  cit != coll.end(); ++cit )
893  if ( cit->type == SubCollision::ELASTIC &&
894  !cit->proj->done() && !cit->targ->done() ) {
895  subevents.push_back(getEl(*cit));
896  if ( !setupFullCollision(subevents.back(), *cit,
897  Nucleon::ELASTIC, Nucleon::ELASTIC) )
898  return false;
899  }
900  return true;
901 }
902 
903 //--------------------------------------------------------------------------
904 
905 // Add all secondary elastic sub-colliions to primary ones.
906 
907 void Angantyr::addELsecond(const multiset<SubCollision> & coll) {
908  // Collect secondary elastic sub-collisions.
909  for ( multiset<SubCollision>::iterator cit = coll.begin();
910  cit != coll.end(); ++cit ) {
911  if ( !cit->proj->done() && cit->type == SubCollision::ELASTIC ) {
912  EventInfo * evp = cit->targ->event();
913  EventInfo add = getEl(*cit);
914  if ( addNucleonExcitation(*evp, add, false) ) {
915  cit->proj->select(*evp, Nucleon::ELASTIC);
916  }
917  }
918  if ( !cit->targ->done() && cit->type == SubCollision::ELASTIC ) {
919  EventInfo * evp = cit->proj->event();
920  EventInfo add = getEl(*cit);
921  if ( addNucleonExcitation(*evp, add, false) ) {
922  cit->targ->select(*evp, Nucleon::ELASTIC);
923  }
924  }
925  }
926 }
927 
928 //--------------------------------------------------------------------------
929 
930 // Shift an event in impact parameter from the nucleon-nucleon
931 // sub-collision to the overall nucleus-nucleus frame. It is assumed
932 // that all partonic vertices are given in units of femtometers.
933 
934 EventInfo & Angantyr::shiftEvent(EventInfo & ei) {
935  if ( HIHooksPtr && HIHooksPtr->canShiftEvent() )
936  return HIHooksPtr->shiftEvent(ei);
937 
938  double ymax = ei.event[1].y();
939  Vec4 bmax = ei.coll->proj->bPos();
940  double ymin = ei.event[2].y();
941  Vec4 bmin = ei.coll->targ->bPos();
942  for ( int i = 0, N = ei.event.size(); i < N; ++i ) {
943  Vec4 shift = bmin + (bmax - bmin)*(ei.event[i].y() - ymin)/(ymax - ymin);
944  ei.event[i].vProdAdd( shift * FM2MM);
945  }
946  return ei;
947 }
948 
949 //--------------------------------------------------------------------------
950 
951 // Prepare a primary sub-collision.
952 
953 bool Angantyr::
954 setupFullCollision(EventInfo & ei, const SubCollision & coll,
955  Nucleon::Status ptype, Nucleon::Status ttype) {
956  if ( !ei.ok ) return false;
957  coll.proj->select(ei, ptype);
958  coll.targ->select(ei, ttype);
959  ei.coll = &coll;
960  ei.projs.clear();
961  ei.projs[coll.proj] = make_pair(1, ei.event.size());
962  ei.targs.clear();
963  ei.targs[coll.targ] = make_pair(2, ei.event.size());
964  shiftEvent(ei);
965  ei.event[1].status(-203);
966  ei.event[1].mother1(1);
967  ei.event[1].mother2(0);
968  ei.event[2].status(-203);
969  ei.event[2].mother1(2);
970  ei.event[2].mother2(0);
971  return fixIsoSpin(ei);
972 }
973 
974 //--------------------------------------------------------------------------
975 
976 // Trace a particle back to one of the beams in an event.
977 
978 int Angantyr::getBeam(Event & ev, int i) {
979  if ( int mom = ev[i].mother1() ) {
980  if ( ev[mom].status() != -203 && ev[mom].mother1() < mom )
981  return getBeam(ev, mom);
982  else
983  return mom;
984  }
985  else
986  return i;
987 }
988 
989 //--------------------------------------------------------------------------
990 
991 // Minimum-bias sub-collisions are always generated as p-p events, and
992 // it is assumed to be safe to be assumed that they are iso-spin
993 // invariant so we can just modify the quark contet in the remnants to
994 // get p-n, n-p, and n-n collisions.
995 
996 bool Angantyr::fixIsoSpin(EventInfo & ei) {
997  if ( HIHooksPtr && HIHooksPtr->canFixIsoSpin() )
998  return HIHooksPtr->fixIsoSpin(ei);
999 
1000  // Check if isospin needs fixing.
1001  int pshift = 0, tshift = 0;
1002  if ( ei.event[1].id() == 2212 && ei.coll->proj->id() == 2112 )
1003  pshift = 1;
1004  if ( ei.event[1].id() == -2212 && ei.coll->proj->id() == -2112 )
1005  pshift = -1;
1006  if ( pshift )
1007  ei.event[1].id(pshift*2112);
1008  if ( ei.event[2].id() == 2212 && ei.coll->targ->id() == 2112 )
1009  tshift = 1;
1010  if ( ei.event[2].id() == -2212 && ei.coll->targ->id() == -2112 )
1011  tshift = -1;
1012  if ( tshift )
1013  ei.event[2].id(tshift*2112);
1014 
1015  if ( !pshift && !tshift ) return true;
1016 
1017  // Try to find corresponding remnants that change flavour
1018  for ( int i = ei.event.size() - 1; i > 2 && ( pshift || tshift ); --i ) {
1019  if ( pshift && ( isRemnant(ei, i) || ei.event[i].status() == 14 )
1020  && getBeam(ei.event, i) == 1 ) {
1021  int newid = 0;
1022  if ( ei.event[i].id() == 2*pshift ) newid = 1*pshift;
1023  if ( ei.event[i].id() == 2101*pshift ) newid = 1103*pshift;
1024  if ( ei.event[i].id() == 2103*pshift ) newid = 1103*pshift;
1025  if ( ei.event[i].id() == 2203*pshift ) newid = 2103*pshift;
1026  if ( ei.event[i].id() == 2212*pshift ) newid = 2112*pshift;
1027  if ( newid ) {
1028  ei.event[i].id(newid);
1029  pshift = 0;
1030  continue;
1031  }
1032  }
1033  if ( tshift && ( isRemnant(ei, i) || ei.event[i].status() == 14 )
1034  && getBeam(ei.event, i) == 2 ) {
1035  int newid = 0;
1036  if ( ei.event[i].id() == 2*tshift ) newid = 1*tshift;
1037  if ( ei.event[i].id() == 2101*tshift ) newid = 1103*tshift;
1038  if ( ei.event[i].id() == 2103*tshift ) newid = 1103*tshift;
1039  if ( ei.event[i].id() == 2203*tshift ) newid = 2103*tshift;
1040  if ( ei.event[i].id() == 2212*tshift ) newid = 2112*tshift;
1041  if ( newid ) {
1042  ei.event[i].id(newid);
1043  tshift = 0;
1044  continue;
1045  }
1046  }
1047  }
1048 
1049  if ( !pshift && !tshift ) return true;
1050 
1051  // Try to find any final state quark that we modify, preferably far
1052  // in the beam direction.
1053  int qselp = 0;
1054  int qselt = 0;
1055  double yselp = 0.0;
1056  double yselt = 0.0;
1057  for ( int i = ei.event.size() - 1; i > 2 && ( pshift || tshift ); --i ) {
1058  if ( pshift && ei.event[i].isFinal() && ei.event[i].id() == 2*pshift) {
1059  if ( ei.event[i].y() > yselp ) {
1060  qselp = i;
1061  yselp = ei.event[i].y();
1062  }
1063  }
1064  if ( tshift && ei.event[i].isFinal() && ei.event[i].id() == 2*tshift) {
1065  if ( ei.event[i].y() < yselt ) {
1066  qselt = i;
1067  yselt = ei.event[i].y();
1068  }
1069  }
1070  }
1071  if ( qselp ) {
1072  ei.event[qselp].id(1*pshift);
1073  pshift = 0;
1074  }
1075  if ( qselt ) {
1076  ei.event[qselt].id(1*tshift);
1077  tshift = 0;
1078  }
1079 
1080  return !pshift && !tshift;
1081 
1082 }
1083 
1084 //--------------------------------------------------------------------------
1085 
1086 // Find recoilers in a primary sub-collisions to conserve energy and
1087 // momentum when adding a secondary one. Not actually used yet.
1088 
1089 vector<int> Angantyr::
1090 findRecoilers(const Event & e, bool tside, int beam, int end,
1091  const Vec4 & pdiff, const Vec4 & pbeam) {
1092  vector<int> ret;
1093  multimap<double,int> ordered;
1094  double mtd2 = pdiff.m2Calc() + pdiff.pT2();
1095  int dir = tside? -1: 1;
1096  double ymax = -log(pdiff.pNeg());
1097  if ( tside ) ymax = -log(pdiff.pPos());
1098  for ( int i = beam, N = end; i < N; ++i )
1099  if ( e[i].status() > 0 )
1100  ordered.insert(make_pair(e[i].y()*dir, i));
1101  Vec4 prec;
1102  double pzfree2 = 0.0;
1103  multimap<double,int>::iterator it = ordered.begin();
1104  // cout << "--- find recoilers ----" << endl;
1105  // cout << setw(10) << setprecision(4)
1106  // << log(pdiff.pPos())
1107  // << setw(10) << pdiff.rap()
1108  // << " diffractive system" << pdiff;
1109  while ( it != ordered.end() ) {
1110  if ( it->first > ymax ) break;
1111  int i = (*it++).second;
1112  Vec4 test = prec + e[i].p();
1113  double mtr2 = test.m2Calc() + test.pT2();
1114  double S = (pbeam + test).m2Calc();
1115  double pz2 = 0.25*(pow2(S - mtr2 - mtd2) - 4.0*mtr2*mtd2)/S;
1116  if ( pz2 < pzfree2 ) {
1117  // cout << setw(10) << setprecision(4) << it->first*dir << " failed "
1118  // << (pz2 < 0.0? -sqrt(-pz2): sqrt(pz2)) << endl;
1119  break;
1120  }
1121  // cout << setw(10) << setprecision(4) << it->first*dir << " accept "
1122  // << sqrt(pz2) << " (" << sqrt(S) << ")" << test;
1123  prec = test;
1124  pzfree2 = pz2;
1125  ret.push_back(i);
1126  }
1127  // cout << "--- found recoilers ---" << endl;
1128 
1129  // *** THINK! *** Is this the best way?
1130  return ret;
1131 
1132 }
1133 
1134 //--------------------------------------------------------------------------
1135 
1136 // Add a secondary sub-collision to a primary one.
1137 
1138 bool Angantyr::addNucleonExcitation(EventInfo & ei, EventInfo & sub,
1139  bool colConnect) {
1140  fixIsoSpin(sub);
1141  shiftEvent(sub);
1142  if ( HIHooksPtr && HIHooksPtr->canAddNucleonExcitation() )
1143  return HIHooksPtr->addNucleonExcitation(ei, sub, colConnect);
1144 
1145  typedef map<Nucleon *, pair<int, int> >::iterator NucPos;
1146  bool tside = false;
1147  NucPos recnuc = ei.projs.find(sub.coll->proj);
1148  if ( recnuc != ei.projs.end() ) tside = true;
1149  NucPos rectarg = ei.targs.find(sub.coll->targ);
1150  if ( rectarg != ei.targs.end() ) {
1151  if ( tside ) infoPtr->errorMsg("Warning from Angantyr::"
1152  "addNucleonExcitation: Nucleon already added.");
1153  tside = false;
1154  recnuc = rectarg;
1155  }
1156 
1157  // First get the projectile system to take recoil and their momentum.
1158  int olddiff = tside? 4: 3;
1159  int beam = tside? 2: 1;
1160  int recbeam = recnuc->second.first;
1161  int recend = recnuc->second.second;
1162  Vec4 pbeam = sub.event[beam].p();
1163  Vec4 pdiff = sub.event[olddiff].p();
1164  if ( sub.code == 106 ) pdiff += sub.event[5].p();
1165  vector<int> rec;
1166  Vec4 prec;
1167  if ( HIHooksPtr && HIHooksPtr->canFindRecoilers() )
1168  rec = HIHooksPtr->findRecoilers(ei.event, tside, recbeam, recend,
1169  pdiff, pbeam);
1170  else if ( recoilerMode == 2 )
1171  rec = findRecoilers(ei.event, tside, recbeam, recend, pdiff, pbeam);
1172  else {
1173  if ( tside && ei.code == 104 && ei.event[4].status() > 0 )
1174  rec.push_back(4);
1175  else if ( !tside && ei.code == 103 && ei.event[3].status() > 0 )
1176  rec.push_back(3);
1177  else if ( tside && ei.event[3].status() > 0 &&
1178  ( ei.code == 102 || ei.code == 106 ) )
1179  rec.push_back(3);
1180  else if ( !tside && ei.event[4].status() > 0 &&
1181  ( ei.code == 102 || ei.code == 106 ) )
1182  rec.push_back(4);
1183  else
1184  for ( int i = recbeam, N = recend; i < N; ++i )
1185  if ( isRemnant(ei, i) && getBeam(ei.event, i) == recbeam )
1186  rec.push_back(i);
1187  }
1188  if ( rec.empty() ) return false;
1189  for ( int i = 0, N = rec.size(); i < N; ++i ) prec += ei.event[rec[i]].p();
1190 
1191  // Find the ransform to the recoilers and the diffractive combined cms
1192  pair<RotBstMatrix,RotBstMatrix> R12;
1193  if ( !getTransforms(prec, pdiff, pbeam, R12, ei.code, sub.code) )
1194  return false;
1195 
1196  // Transform the recoilers.
1197  for ( int i = 0, N = rec.size(); i < N; ++i )
1198  ei.event[rec[i]].rotbst(R12.first);
1199 
1200  // Copy the event and transform and offset the particles appropriately.
1201 
1202  int newbeam = ei.event.size();
1203  ei.event.append(sub.event[beam]);
1204  ei.event.back().status(-203);
1205  ei.event.back().mother1(beam);
1206  ei.event.back().mother2(0);
1207  ei.event.back().daughter1(ei.event.size());
1208  int newdiff = ei.event.size();
1209  int nextpos = 5;
1210  ei.event.append(sub.event[olddiff]);
1211  ei.event.back().rotbst(R12.second);
1212  ei.event.back().mother1(newbeam);
1213  ei.event.back().mother2(0);
1214  if ( sub.code == 102 ) {
1215  if ( tside )
1216  ei.targs[sub.coll->targ] = make_pair(newbeam, ei.event.size());
1217  else
1218  ei.projs[sub.coll->proj] = make_pair(newbeam, ei.event.size());
1219  return true;
1220  }
1221 
1222  int idoff = tside? newdiff - olddiff: newdiff - olddiff - 1;
1223  if ( sub.code == 106 ) {
1224  // Special handling of central diffraction.
1225  ++newdiff;
1226  ++nextpos;
1227  idoff = newdiff - 5;
1228  ei.event.append(sub.event[5]);
1229  ei.event.back().rotbst(R12.second);
1230  ei.event.back().mother1(newbeam);
1231  ei.event.back().mother2(0);
1232  }
1233  ei.event.back().daughter1(sub.event[olddiff].daughter1() + idoff);
1234  ei.event.back().daughter2(sub.event[olddiff].daughter2() + idoff);
1235  int coloff = ei.event.lastColTag();
1236  // Add energy to zeroth line and calculate new invariant mass.
1237  ei.event[0].p( ei.event[0].p() + pbeam );
1238  ei.event[0].m( ei.event[0].mCalc() );
1239  for (int i = nextpos; i < sub.event.size(); ++i) {
1240  Particle temp = sub.event[i];
1241 
1242  // Add offset to nonzero mother, daughter and colour indices.
1243  if ( temp.mother1() == olddiff ) temp.mother1(newdiff);
1244  else if ( temp.mother1() > 0 ) temp.mother1(temp.mother1() + idoff );
1245  if ( temp.mother2() == olddiff ) temp.mother2(newdiff);
1246  else if ( temp.mother2() > 0 ) temp.mother2(temp.mother2() + idoff );
1247  if ( temp.daughter1() > 0 ) temp.daughter1( temp.daughter1() + idoff );
1248  if ( temp.daughter2() > 0 ) temp.daughter2( temp.daughter2() + idoff );
1249  if ( temp.col() > 0 ) temp.col( temp.col() + coloff );
1250  if ( temp.acol() > 0 ) temp.acol( temp.acol() + coloff );
1251  temp.rotbst(R12.second);
1252  // Append particle to summed event.
1253  ei.event.append( temp );
1254  }
1255 
1256  addJunctions(ei.event, sub.event, coloff);
1257 
1258  if ( tside )
1259  ei.targs[sub.coll->targ] = make_pair(newbeam, ei.event.size());
1260  else
1261  ei.projs[sub.coll->proj] = make_pair(newbeam, ei.event.size());
1262  return true;
1263 
1264 }
1265 
1266 //--------------------------------------------------------------------------
1267 
1268 // Calculate boosts to shuffle momenta when adding secondary
1269 // sub-collisions.
1270 
1271 bool
1272 Angantyr::getTransforms(Vec4 prec, Vec4 pdiff, const Vec4 & pbeam,
1273  pair<RotBstMatrix,RotBstMatrix> & R12,
1274  int code1, int code2) {
1275  code1 += code2;
1276  // static ofstream os("testhi.out");
1277  // os << "=== " << code1 << "+" << code2 << " ===" << endl;
1278  // os << prec << pbeam << pdiff;
1279  RotBstMatrix Ri;
1280  Ri.toCMframe(pbeam, prec);
1281  Vec4 pr1 = prec;
1282  Vec4 pb1 = pbeam;
1283  Vec4 pd1 = pdiff;
1284  pr1.rotbst(Ri);
1285  pb1.rotbst(Ri);
1286  pd1.rotbst(Ri);
1287  // os << "=>" << endl << pr1 << pb1 << pd1;
1288  Vec4 pr2 = pr1;
1289  if ( pd1.pT() >= abs(pr2.pz()) ) {
1290  // os << "*** failed to rotate ***" << endl;
1291  return false;
1292  }
1293  double the = asin(pd1.pT()/abs(pr2.pz()));
1294  RotBstMatrix R1;
1295  R1.rot(the, pd1.phi());
1296  pr2.rotbst(R1);
1297  // os << "=>" << endl << pr2 << pd1;
1298 
1299  double S = (prec + pbeam).m2Calc();
1300  double mtr2 = pr2.pT2() + pr2.m2Calc();
1301  double mtd2 = pd1.pT2() + pd1.m2Calc();
1302  if ( sqrt(S) <= sqrt(mtr2) + sqrt(mtd2) ) {
1303  // os << "*** failed to boost ***" << endl;
1304  return false;
1305  }
1306  double z2 = 0.25*(mtr2*mtr2 + (mtd2 - S)*(mtd2 - S) - 2.0*mtr2*(mtd2 + S))/S;
1307  if ( z2 <= 0.0 ) {
1308  // os << "*** failed to boost ***" << endl;
1309  return false;
1310  }
1311  double z = sqrt(z2);
1312  double ppo2 = pow2(pr2.pNeg());
1313  double ppn2 = pow2(z + sqrt(z2 + mtr2));
1314  R1.bst(0.0, 0.0, -(ppn2 - ppo2)/(ppn2 + ppo2));
1315 
1316  ppo2 = pow2(pd1.pPos());
1317  ppn2 = pow2(z + sqrt(z2 + mtd2));
1318  RotBstMatrix R2;
1319  R2.bst(0.0, 0.0, (ppn2 - ppo2)/(ppn2 + ppo2));
1320  Vec4 pr3 = pr1;
1321  pr3.rotbst(R1);
1322  Vec4 pd3 = pd1;
1323  pd3.rotbst(R2);
1324  // os << "=>" << endl << pr3 << pd3;
1325  RotBstMatrix Rf = Ri;
1326  Rf.invert();
1327  Vec4 pr4 = pr3;
1328  pr4.rotbst(Rf);
1329  Vec4 pd4 = pd3;
1330  pd4.rotbst(Rf);
1331  // os << "=>" << endl << pr4 << pd4 << pr4 + pd4 << prec + pbeam;
1332 
1333  R12.first = R12.second = Ri;
1334  R12.first.rotbst(R1);
1335  R12.second.rotbst(R2);
1336  R12.first.rotbst(Rf);
1337  R12.second.rotbst(Rf);
1338  prec.rotbst(R12.first);
1339  pdiff.rotbst(R12.second);
1340  // os << prec << pdiff;
1341 
1342  return true;
1343 
1344 }
1345 
1346 //--------------------------------------------------------------------------
1347 
1348 // Add sub-events together taking special care with the status of the
1349 // incoming nucleons, and also handle the junctions correctly.
1350 
1351 void Angantyr::addSubEvent(Event & evnt, Event & sub) {
1352 
1353  int idoff = evnt.size() - 1;
1354  int coloff = evnt.lastColTag();
1355 
1356  for (int i = 1; i < sub.size(); ++i) {
1357  Particle temp = sub[i];
1358 
1359  // Add offset to nonzero mother, daughter and colour indices.
1360  if ( temp.status() == -203 )
1361  temp.status(-13);
1362  else {
1363  if ( temp.mother1() > 0 ) temp.mother1(temp.mother1() + idoff );
1364  if ( temp.mother2() > 0 ) temp.mother2( temp.mother2() + idoff );
1365  }
1366  if ( temp.daughter1() > 0 ) temp.daughter1( temp.daughter1() + idoff );
1367  if ( temp.daughter2() > 0 ) temp.daughter2( temp.daughter2() + idoff );
1368  if ( temp.col() > 0 ) temp.col( temp.col() + coloff );
1369  if ( temp.acol() > 0 ) temp.acol( temp.acol() + coloff );
1370  // Append particle to summed event.
1371  evnt.append( temp );
1372  }
1373 
1374  addJunctions(evnt, sub, coloff);
1375 
1376 }
1377 
1378 void Angantyr::addJunctions(Event & ev, Event & addev, int coloff) {
1379 
1380  // Read out junctions one by one.
1381  Junction tempJ;
1382  int begCol, endCol;
1383  for (int i = 0; i < addev.sizeJunction(); ++i) {
1384  tempJ = addev.getJunction(i);
1385 
1386  // Add colour offsets to all three legs.
1387  for (int j = 0; j < 3; ++j) {
1388  begCol = tempJ.col(j);
1389  endCol = tempJ.endCol(j);
1390  if (begCol > 0) begCol += coloff;
1391  if (endCol > 0) endCol += coloff;
1392  tempJ.cols( j, begCol, endCol);
1393  }
1394  // Append junction to summed event.
1395  ev.appendJunction( tempJ );
1396  }
1397 }
1398 
1399 //--------------------------------------------------------------------------
1400 
1401 // Special function to generatee secondary absorptive events as single
1402 // diffraction. Called from Angantyr::next() and used for debugging
1403 // and tuning purposes.
1404 
1405 bool Angantyr::nextSASD(int procid) {
1406  Nucleon dummy;
1407  double bp = pythia[SASD]->parm("Angantyr:SDTestB");
1408  SubCollision coll(dummy, dummy, bp*collPtr->avNDB(), bp, SubCollision::ABS);
1409  EventInfo ei = getSASD(&coll, procid);
1410  if ( !ei.ok ) return false;
1411  pythia[HADRON]->event = ei.event;
1412  updateInfo();
1413  if ( pythia[HADRON]->flag("HadronLevel:all") ) {
1414  if ( HIHooksPtr && HIHooksPtr->canForceHadronLevel() ) {
1415  if ( !HIHooksPtr->forceHadronLevel(*pythia[HADRON]) ) return false;
1416  } else {
1417  if ( !pythia[HADRON]->forceHadronLevel(false) ) return false;
1418  }
1419  }
1420  return true;
1421 }
1422 
1423 //--------------------------------------------------------------------------
1424 
1425 // Take all sub-events and merge them together.
1426 
1427 bool Angantyr::buildEvent(list<EventInfo> & subevents,
1428  const vector<Nucleon> & proj,
1429  const vector<Nucleon> & targ) {
1430  Event & etmp = pythia[HADRON]->event;
1431  etmp.reset();
1432  etmp.append(projPtr->produceIon(false));
1433  etmp.append(targPtr->produceIon(true));
1434  etmp[0].p(etmp[1].p() + etmp[2].p());
1435  etmp[0].m(etmp[0].mCalc());
1436 
1437  // Start with the signal event(s)
1438  if ( hasSignal ) {
1439  bool found = false;
1440  for ( list<EventInfo>::iterator sit = subevents.begin();
1441  sit != subevents.end(); ++sit ) {
1442  if ( sit->code >= 101 && sit->code <= 106 ) continue;
1443  addSubEvent(etmp, sit->event);
1444  hiInfo.select(sit->info);
1445  hiInfo.addSubCollision(*sit->coll);
1446  subevents.erase(sit);
1447  found = true;
1448  break;
1449  }
1450  if ( !found ) {
1451  infoPtr->errorMsg("Warning from Angantyr::next:"
1452  " Failed to generate signal event.");
1453  return false;
1454  }
1455  } else
1456  hiInfo.select(subevents.begin()->info);
1457 
1458  // Then all the others
1459  for ( list<EventInfo>::iterator sit = subevents.begin();
1460  sit != subevents.end(); ++sit ) {
1461  addSubEvent(etmp, sit->event);
1462  hiInfo.addSubCollision(*sit->coll);
1463  }
1464 
1465  // Finally add all nucleon remnants.
1466  return addNucleusRemnants(proj, targ);
1467 
1468 }
1469 
1470 //--------------------------------------------------------------------------
1471 
1472 // Construct nucleus remnants fron all non-interacting nucleons and
1473 // add them to the main event.
1474 
1475 bool Angantyr::
1476 addNucleusRemnants(const vector<Nucleon> & proj,
1477  const vector<Nucleon> & targ) {
1478  Event & etmp = pythia[HADRON]->event;
1479  int npp = 0;
1480  int nnp = 0;
1481  Vec4 ppsum;
1482  for ( int i = 0, N = proj.size(); i< N; ++i ) {
1483  if ( proj[i].event() ) hiInfo.addProjectileNucleon(proj[i]);
1484  else {
1485  double e = pythia[HADRON]->parm("Beams:eA");
1486  double m = pythia[HADRON]->particleData.m0(proj[i].id());
1487  double pz = sqrt(max(e*e - m*m, 0.0));
1488  if ( proj[i].id() == 2212 ) {
1489  ++npp;
1490  ppsum += Vec4(0.0, 0.0, pz, e);
1491  } else if ( proj[i].id() == 2112 ) {
1492  ++nnp;
1493  ppsum += Vec4(0.0, 0.0, pz, e);
1494  } else
1495  etmp.append(proj[i].id(), 14, 1, 0, 0, 0, 0, 0, 0.0, 0.0, pz, e, m);
1496  }
1497  }
1498  int npt = 0;
1499  int nnt = 0;
1500  Vec4 tpsum;
1501  for ( int i = 0, N = targ.size(); i< N; ++i ) {
1502  if ( targ[i].event() ) hiInfo.addTargetNucleon(targ[i]);
1503  else {
1504  double e = pythia[HADRON]->parm("Beams:eB");
1505  double m = pythia[HADRON]->particleData.m0(targ[i].id());
1506  double pz = -sqrt(max(e*e - m*m, 0.0));
1507  if ( targ[i].id() == 2212 ) {
1508  ++npt;
1509  tpsum += Vec4(0.0, 0.0, pz, e);
1510  } else if ( targ[i].id() == 2112 ) {
1511  ++nnt;
1512  tpsum += Vec4(0.0, 0.0, pz, e);
1513  } else
1514  etmp.append(targ[i].id(), 14, 2, 0, 0, 0, 0, 0, 0.0, 0.0, pz, e, m);
1515  }
1516  }
1517 
1518  Vec4 ptot = etmp[0].p();
1519  for ( int i = 0, N = etmp.size(); i < N; ++i )
1520  if ( etmp[i].status() > 0 ) ptot -= etmp[i].p();
1521 
1522  if ( npp + nnp +npt + nnt == 0 ) return true;
1523  ParticleData & pdt = pythia[HADRON]->particleData;
1524  int idp = 0;
1525  if ( npp + nnp > 1 ) {
1526  idp = 1000000009 + 10000*npp + 10*(nnp + npp);
1527  pdt.addParticle(idp, "NucRem", 0, 3*npp, 0, ppsum.mCalc());
1528  pdt.particleDataEntryPtr(idp)->setHasChanged(false);
1529  }
1530  else if ( npp == 1 ) idp = 2212;
1531  else if ( nnp == 1 ) idp = 2112;
1532  int idt = 0;
1533  if ( npt + nnt > 1 ) {
1534  idt = 1000000009 + 10000*npt + 10*(nnt + npt);
1535  pdt.addParticle(idt, "NucRem", 0, 3*npt, 0, tpsum.mCalc());
1536  pdt.particleDataEntryPtr(idt)->setHasChanged(false);
1537  }
1538  else if ( npt == 1 ) idt = 2212;
1539  else if ( nnt == 1 ) idt = 2112;
1540 
1541  if ( npp + nnp > npt + nnt ) {
1542  if ( npt + nnt > 0 ) {
1543  etmp.append(idt, 14, 2, 0, 0, 0, 0, 0, tpsum, tpsum.mCalc());
1544  ptot -= tpsum;
1545  }
1546  etmp.append(idp, 14, 1, 0, 0, 0, 0, 0, ptot, ptot.mCalc());
1547  } else {
1548  if ( npp + nnp > 0 ) {
1549  etmp.append(idp, 14, 1, 0, 0, 0, 0, 0, ppsum, ppsum.mCalc());
1550  ptot -= ppsum;
1551  }
1552  etmp.append(idt, 14, 2, 0, 0, 0, 0, 0, ptot, ptot.mCalc());
1553  }
1554  return true;
1555 }
1556 
1557 //--------------------------------------------------------------------------
1558 
1559 // The main method called from Pythia::next().
1560 
1561 bool Angantyr::next() {
1562 
1563  if ( flag("Angantyr:SDTest") ) return nextSASD(104);
1564 
1565  int itry = MAXTRY;
1566 
1567  while ( itry-- ) {
1568 
1569  // Generate nuclei, impact paramter and nucleon sub-collisions.
1570  projectile = projPtr->generate();
1571  target = targPtr->generate();
1572 
1573  double bweight = 0.0;
1574  Vec4 bvec = bGenPtr->generate(bweight);
1575  double T = 0.0;
1576  subColls = collPtr->getCollisions(projectile, target, bvec, T);
1577  hiInfo.addAttempt(T, bvec.pT(), bweight);
1578  hiInfo.subCollisionsPtr(&subColls);
1579  if ( flag("Angantyr:GlauberOnly") ) return true;
1580  if ( subColls.empty() ) continue;
1581 
1582 
1583  list<EventInfo> subevents;
1584 
1585  if ( !genAbs(subColls, subevents) ) {
1586  infoPtr->errorMsg("Warning from PyHIia::next: "
1587  "Could not setup signal or ND collisions.");
1588  continue;
1589  }
1590  if ( hasSignal && subevents.empty() ) continue;
1591 
1592  // Collect absorptively wounded nucleons in secondary
1593  // sub-collisions.
1594  addSASD(subColls);
1595 
1596  // Collect full double diffraction collisions.
1597  if ( !addDD(subColls, subevents) ) {
1598  infoPtr->errorMsg("Warning from PyHIia::next:"
1599  " Could not setup DD sub collision.");
1600  continue;
1601  }
1602 
1603  // Collect full single diffraction collisions.
1604  if ( !addSD(subColls, subevents) ) {
1605  infoPtr->errorMsg("Warning from PyHIia::next:"
1606  " Could not setup SD sub collision.");
1607  continue;
1608  }
1609 
1610  // Collect secondary single diffractive sub-collisions.
1611  addSDsecond(subColls);
1612 
1613  // Collect full central diffraction collisions.
1614  if ( !addCD(subColls, subevents) ) {
1615  infoPtr->errorMsg("Warning from PyHIia::next:"
1616  " Could not setup CD sub collisions.");
1617  continue;
1618  }
1619 
1620  // Collect secondary central diffractive sub-collisions.
1621  addCDsecond(subColls);
1622 
1623  // Collect full elastic collisions.
1624  if ( !addEL(subColls, subevents) ) {
1625  infoPtr->errorMsg("Warning from PyHIia::next:"
1626  " Could not setup elastic sub collisions.");
1627  continue;
1628  }
1629 
1630  // Collect secondary elastic sub-collisions.
1631  addELsecond(subColls);
1632 
1633  // Finally bunch all events together.
1634  if ( subevents.empty() ) continue;
1635 
1636  if ( !buildEvent(subevents, projectile, target) ) continue;
1637 
1638  // Finally we hadronise everything, if requested.
1639  if ( pythia[HADRON]->flag("HadronLevel:all") ) {
1640  if ( HIHooksPtr && HIHooksPtr->canForceHadronLevel() ) {
1641  if ( !HIHooksPtr->forceHadronLevel(*pythia[HADRON]) ) continue;
1642  } else {
1643  if ( !pythia[HADRON]->forceHadronLevel(false) ) continue;
1644  }
1645  }
1646 
1647  hiInfo.accept();
1648 
1649  updateInfo();
1650 
1651 
1652  return true;
1653 
1654  }
1655 
1656  infoPtr->errorMsg("Abort from Angantyr::next: Too many "
1657  "attempts to generate a working impact parameter point. "
1658  "Consider reducing HeavyIon:bWidth.");
1659  hiInfo.reject();
1660  return false;
1661 
1662 }
1663 
1664 //==========================================================================
1665 
1666 } // end namespace Pythia8
Pythia8::Angantyr::getTransforms
static bool getTransforms(Vec4 p1, Vec4 p2, const Vec4 &p1p, pair< RotBstMatrix, RotBstMatrix > &R12, int, int)
Definition: HeavyIons.cc:1293
Pythia8::Nucleon::ABS
The nucleon is absorptively wounded.
Definition: HIUserHooks.h:73
Pythia8::SubCollision::ABS
This is an absorptive (non-diffractive) collision.
Definition: HIUserHooks.h:202
Pythia8::SubCollision::CDE
Both excited but with central diffraction.
Definition: HIUserHooks.h:201
Pythia8::ImpactParameterGenerator::init
virtual bool init()
Virtual init method.
Definition: HIUserHooks.cc:238
Pythia8::HeavyIons::HIHooksPtr
HIUserHooks * HIHooksPtr
Definition: HeavyIons.h:108
Pythia8::HeavyIons::sigtot
SigmaTotal sigtot
Definition: HeavyIons.h:104
Pythia8::SubCollision::DDE
Both projectile and target are diffractively excited.
Definition: HIUserHooks.h:200
Pythia8::Angantyr::SIGNN
Optional object for signal processes (nn).
Definition: HeavyIons.h:145
Pythia8::Angantyr::next
virtual bool next()
Produce a collision involving heavy ions.
Definition: HeavyIons.cc:1582
Pythia8::HIUserHooks::hasProjectileModel
virtual bool hasProjectileModel() const
A suser-supplied NucleusModel for the projectile and target.
Definition: HIUserHooks.h:1056
Pythia8::Angantyr::addNucleonExcitation
bool addNucleonExcitation(EventInfo &orig, EventInfo &add, bool colConnect=false)
Definition: HeavyIons.cc:1158
Pythia8::EventInfo::code
int code
The code for the subprocess.
Definition: HIUserHooks.h:830
Pythia8::NucleusModel::initPtr
void initPtr(int idIn, Settings &settingsIn, ParticleData &particleDataIn, Rndm &rndIn)
Init method.
Definition: HIUserHooks.cc:23
Pythia8::SubCollisionModel::avNDB
double avNDB() const
Return the average non-diffractive impact parameter.
Definition: HIUserHooks.h:548
Pythia8::HIInfo::sigmaNDErr
double sigmaNDErr() const
The estimated statistical error on sigmaND().
Definition: HIUserHooks.h:861
Pythia8::HeavyIons::addSpecialSettings
static void addSpecialSettings(Settings &settings)
Definition: HeavyIons.cc:24
Pythia8::HeavyIons::info
vector< Info * > info
Definition: HeavyIons.h:119
Pythia8::Angantyr::getSignal
EventInfo getSignal(const SubCollision &coll)
Generate events from the internal Pythia oblects;.
Definition: HeavyIons.cc:594
Pythia8::EventInfo::ordering
double ordering
The ordering variable of this event.
Definition: HIUserHooks.h:802
Pythia8::HIUserHooks::canFixIsoSpin
virtual bool canFixIsoSpin() const
Definition: HIUserHooks.h:1072
Pythia8::Angantyr::nextSASD
bool nextSASD(int proc)
Generate a single diffractive.
Definition: HeavyIons.cc:1426
Pythia8::Angantyr::SIGPP
Optional object for signal processes (pp).
Definition: HeavyIons.h:142
Pythia8::Angantyr::PythiaObject
PythiaObject
Enumerate the different internal Pythia objects.
Definition: HeavyIons.h:138
Pythia8::EventInfo::projs
map< Nucleon *, pair< int, int > > projs
Definition: HIUserHooks.h:815
Pythia8::HIUserHooks::canFindRecoilers
virtual bool canFindRecoilers() const
Definition: HIUserHooks.h:1092
Pythia8::HIUserHooks::hasImpactParameterGenerator
virtual bool hasImpactParameterGenerator() const
A user-supplied impact parameter generator.
Definition: HIUserHooks.h:1051
Pythia8::HIUserHooks::canAddNucleonExcitation
bool canAddNucleonExcitation() const
Definition: HIUserHooks.h:1081
Pythia8::EventInfo::event
Event event
The Event and info objects.
Definition: HIUserHooks.h:798
Pythia8::Angantyr::MBIAS
Minimum Bias processed.
Definition: HeavyIons.h:140
Pythia8::Angantyr::mkEventInfo
EventInfo mkEventInfo(Pythia &pyt, const SubCollision *coll=0)
Setup an EventInfo object from a Pythia instance.
Definition: HeavyIons.cc:319
Pythia8::SubCollisionModel::sigND
double sigND() const
The non-diffractive (absorptive) cross section.
Definition: HIUserHooks.h:537
Pythia8::HIInfo::failedExcitation
void failedExcitation()
Register a failed nucleon excitation.
Definition: HIUserHooks.h:948
Pythia8::HIUserHooks::canShiftEvent
virtual bool canShiftEvent() const
A user-supplied method for shifting the event in impact parameter space.
Definition: HIUserHooks.h:1076
Pythia8::SubCollisionModel::getCollisions
virtual multiset< SubCollision > getCollisions(vector< Nucleon > &proj, vector< Nucleon > &targ, const Vec4 &bvec, double &T)=0
Definition: HIUserHooks.cc:514
Pythia8::Angantyr::ALL
Indicates all objects.
Definition: HeavyIons.h:146
Pythia8::Angantyr::HADRON
For hadronization only.
Definition: HeavyIons.h:139
Pythia8::SubCollision::targ
Nucleon * targ
The target nucleon.
Definition: HIUserHooks.h:226
Pythia8::HeavyIons::pythia
vector< Pythia * > pythia
Definition: HeavyIons.h:112
Pythia8::HeavyIons::clearProcessLevel
void clearProcessLevel(Pythia &pyt)
Definition: HeavyIons.cc:126
Pythia8::SubCollision::ELASTIC
This is an elastic scattering.
Definition: HIUserHooks.h:197
Pythia8::SubCollision::SDEP
The projectile is diffractively excited.
Definition: HIUserHooks.h:198
Pythia8::NucleusModel::generate
virtual vector< Nucleon > generate() const =0
Pythia8::Nucleon::ELASTIC
The nucleon is elastically scattered.
Definition: HIUserHooks.h:71
Pythia8::Angantyr::SASD
Single diffractive as one side of non-diffractive.
Definition: HeavyIons.h:141
Pythia8::HIInfo::sigmaND
double sigmaND() const
Definition: HIUserHooks.h:856
Pythia8::HeavyIons::sumUpMessages
void sumUpMessages(Info &in, string tag, const Info &other)
Definition: HeavyIons.cc:115
Pythia8::HIUserHooks::hasEventOrdering
virtual bool hasEventOrdering() const
A user-supplied ordering of events in (inverse) hardness.
Definition: HIUserHooks.h:1067
Pythia8::HIUserHooks::hasSubCollisionModel
virtual bool hasSubCollisionModel()
Definition: HIUserHooks.h:1063
Pythia8::HeavyIons::setupSpecials
static void setupSpecials(Settings &settings, string match)
Duplicate setting on the form match: to settings on the form HImatch:
Definition: HeavyIons.cc:32
Pythia8::Angantyr::SIGPN
Optional object for signal processes (pn).
Definition: HeavyIons.h:143
Pythia8::Angantyr::init
virtual bool init()
Initialize Angantyr.
Definition: HeavyIons.cc:340
Pythia8::HIUserHooks::init
virtual void init(int idProjIn, int idTargIn)
Initialize this user hook.
Definition: HIUserHooks.h:1045
Pythia8::EventInfo::ok
bool ok
Is the event properly generated?
Definition: HIUserHooks.h:811
Pythia8::Angantyr::SIGNP
Optional object for signal processes (np).
Definition: HeavyIons.h:144
Pythia8::SubCollision::SDET
The target is diffractively excited.
Definition: HIUserHooks.h:199
Pythia8::EventInfo::coll
const SubCollision * coll
The associated SubCollision object.
Definition: HIUserHooks.h:808
Pythia8::HIInfo::sigmaTot
double sigmaTot() const
The Monte Carlo integrated total cross section in the current run.
Definition: HIUserHooks.h:845
Pythia8::Angantyr::Angantyr
Angantyr(Pythia &mainPythiaIn)
Definition: HeavyIons.cc:272
Pythia8::Nucleon::DIFF
The nucleon is diffractively wounded.
Definition: HIUserHooks.h:72
beam
Definition: beam.h:43
Pythia8::HIInfo::nAttempts
long nAttempts() const
The number of attempted impact parameter points.
Definition: HIUserHooks.h:866
Event
Definition: AgUStep.h:26
Pythia8::Angantyr::setUserHooksPtr
bool setUserHooksPtr(PythiaObject sel, UserHooks *userHooksPtrIn)
Set UserHooks for specific (or ALL) internal Pythia objects.
Definition: HeavyIons.cc:307
Pythia8::Angantyr::addNucleusRemnants
bool addNucleusRemnants(const vector< Nucleon > &proj, const vector< Nucleon > &targ)
Definition: HeavyIons.cc:1497
Pythia8::ImpactParameterGenerator::generate
virtual Vec4 generate(double &weight) const
Definition: HIUserHooks.cc:260
Pythia8::HIInfo::sigmaTotErr
double sigmaTotErr() const
The estimated statistical error on sigmaTot().
Definition: HIUserHooks.h:850
Pythia8::Nucleon::Status
Status
Enum for specifying the status of a nucleon.
Definition: HIUserHooks.h:69
Pythia8::HeavyIons::pythiaNames
vector< string > pythiaNames
The names associated with the secondary pythia objects.
Definition: HeavyIons.h:115
Pythia8::Angantyr::addSubEvent
void addSubEvent(Event &evnt, Event &sub)
Add a sub-event to the final event record.
Definition: HeavyIons.cc:1372
Nucleon
Definition: Nucleon.h:9
Pythia8::HIUserHooks::canForceHadronLevel
virtual bool canForceHadronLevel() const
A user supplied wrapper around the Pythia::forceHadronLevel()
Definition: HIUserHooks.h:1086
Pythia8::HIInfo::weight
double weight() const
The weight for this collision.
Definition: HIUserHooks.h:937
Pythia8::HeavyIons::isHeavyIon
static bool isHeavyIon(Settings &settings)
Definition: HeavyIons.cc:258
Pythia8::SubCollisionModel::init
virtual bool init()
Virtual init method.
Definition: HIUserHooks.cc:278
Generated by   doxygen 1.8.5