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