StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
DireSplitInfo.cc
1 // DireSplitInfo.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Stefan Prestel, 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 header) for the Dire classes
7 // use as storage containers for splitting information.
8 
9 #include "Pythia8/DireSplitInfo.h"
10 
11 namespace Pythia8 {
12 
13 //==========================================================================
14 
15 // Definition of DireSplitKinematics members.
16 
17 void DireSplitKinematics::store ( const DireSplitKinematics& k) {
18  m2Dip = k.m2Dip;
19  pT2 = k.pT2;
20  pT2Old = k.pT2Old;
21  z = k.z;
22  phi = k.phi;
23  sai = k.sai;
24  xa = k.xa;
25  phi2 = k.phi2;
26  m2RadBef = k.m2RadBef;
27  m2Rec = k.m2Rec;
28  m2RadAft = k.m2RadAft;
29  m2EmtAft = k.m2EmtAft;
30  m2EmtAft2 = k.m2EmtAft2;
31  xBef = k.xBef;
32  xAft = k.xAft;
33 }
34 
35 void DireSplitKinematics::list() {
36  cout << "List DireSplitKinematics:"
37  << scientific << setprecision(3) << "\n"
38  << " m2Dip = " << m2Dip << "\n"
39  << " pT2 = " << pT2 << "\t"
40  << " z = " << z << "\t"
41  << " phi = " << phi << "\n"
42  << " sai = " << sai << "\t"
43  << " xa = " << xa << "\t"
44  << " phi2 = " << phi2 << "\n"
45  << " m2RadBef = " << m2RadBef << " "
46  << " m2Rec = " << m2Rec << " "
47  << " m2RadAft = " << m2RadAft << " "
48  << " m2EmtAft = " << m2EmtAft << " "
49  << " m2EmtAft2t = " << m2EmtAft2 << "\n";
50 }
51 
52 //==========================================================================
53 
54 // Definition of DireSplitInfo members.
55 
56 void DireSplitInfo::init(const Event& state) {
57  if (iRadBef>0) particleSave.push_back(DireSplitParticle(state[iRadBef]));
58  else particleSave.push_back(DireSplitParticle());
59  if (iRecBef>0) particleSave.push_back(DireSplitParticle(state[iRecBef]));
60  else particleSave.push_back(DireSplitParticle());
61  if (iRadAft>0) particleSave.push_back(DireSplitParticle(state[iRadAft]));
62  else particleSave.push_back(DireSplitParticle());
63  if (iRecAft>0) particleSave.push_back(DireSplitParticle(state[iRecAft]));
64  else particleSave.push_back(DireSplitParticle());
65  if (iEmtAft>0) particleSave.push_back(DireSplitParticle(state[iEmtAft]));
66  else particleSave.push_back(DireSplitParticle());
67  if (iEmtAft2>0)particleSave.push_back(DireSplitParticle(state[iEmtAft2]));
68  else particleSave.push_back(DireSplitParticle());
69 }
70 
71 void DireSplitInfo::store (const DireSplitInfo& s) {
72  clear();
73  kinSave.clear();
74  particleSave.resize(0);
75  extras.clear();
76  iRadBef = s.iRadBef;
77  iRecBef = s.iRecBef;
78  iRadAft = s.iRadAft;
79  iRecAft = s.iRecAft;
80  iEmtAft = s.iEmtAft;
81  iEmtAft2 = s.iEmtAft2;
82  for (int i=0; i < int(s.particleSave.size()); ++i)
83  particleSave.push_back(s.particleSave[i]);
84  kinSave.store(s.kinSave);
85  side = s.side;
86  type = s.type;
87  system = s.system;
88  systemRec = s.systemRec;
89  splittingSelName = s.splittingSelName;
90  for ( unordered_map<string,double>::const_iterator it = s.extras.begin();
91  it != s.extras.end(); ++it )
92  extras.insert(make_pair(it->first,it->second));
93  useForBranching = s.useForBranching;
94  terminateEvolution = s.terminateEvolution;
95  iSiblings = s.iSiblings;
96 }
97 
98 void DireSplitInfo::save () {
99  kinSaveStore = kinSave;
100  particleSaveStore = particleSave;
101  extrasStore = extras;
102  iRadBefStore = iRadBef;
103  iRecBefStore = iRecBef;
104  iRadAftStore = iRadAft;
105  iRecAftStore = iRecAft;
106  iEmtAftStore = iEmtAft;
107  iEmtAft2Store = iEmtAft2;
108  sideStore = side;
109  typeStore = type;
110  systemStore = system;
111  systemRecStore = systemRec;
112  splittingSelNameStore = splittingSelName;
113  useForBranchingStore = useForBranching;
114  terminateEvolutionStore = terminateEvolution;
115  iSiblingsStore = iSiblings;
116 }
117 
118 void DireSplitInfo::restore () {
119  kinSave = kinSaveStore;
120  particleSave = particleSaveStore;
121  extras = extrasStore;
122  iRadBef = iRadBefStore;
123  iRecBef = iRecBefStore;
124  iRadAft = iRadAftStore;
125  iRecAft = iRecAftStore;
126  iEmtAft = iEmtAftStore;
127  iEmtAft2 = iEmtAft2Store;
128  side = sideStore;
129  type = typeStore;
130  system = systemStore;
131  systemRec = systemRecStore;
132  splittingSelName = splittingSelNameStore;
133  useForBranching = useForBranchingStore;
134  terminateEvolution = terminateEvolutionStore;
135  iSiblings = iSiblingsStore;
136 }
137 
138 void DireSplitInfo::storeInfo(string name, int typeIn, int systemIn,
139  int systemRecIn,
140  int sideIn, int iPosRadBef, int iPosRecBef,
141  const Event& state, int idEmtAft,
142  int idRadAft, int nEmissions, double m2Dip, double pT2, double pT2Old,
143  double z, double phi, double m2Bef, double m2s, double m2r, double m2i,
144  double sa1, double xa, double phia1, double m2j, double xBef, double xAft) {
145  clear();
146  storeName(name);
147  storeType(typeIn);
148  storeSystem(systemIn);
149  storeSystemRec(systemRecIn);
150  storeSide(sideIn);
151  storeRadRecBefPos(iPosRadBef, iPosRecBef);
152  storeRadBef(state[iPosRadBef]);
153  storeRecBef(state[iPosRecBef]);
154  setEmtAft(idEmtAft);
155  setRadAft(idRadAft);
156  if (nEmissions == 2) set2to4kin( m2Dip, pT2, z, phi, sa1, xa, phia1,
157  m2Bef, m2s, m2r, m2i, m2j);
158  else set2to3kin( m2Dip, pT2, z, phi, m2Bef, m2s, m2r, m2i);
159  storeExtras(
160  unordered_map<string,double>(create_unordered_map<string,double>
161  ("iRadBef",iPosRadBef)("iRecBef",iPosRecBef)("idRadAft",idRadAft)) );
162  set_pT2Old(pT2Old);
163  set_xBef(xBef);
164  set_xAft(xAft);
165 }
166 
167 void DireSplitInfo::storePosAfter( int iRadAftIn, int iRecAftIn, int iEmtAftIn,
168  int iEmtAft2In) {
169  iRadAft = iRadAftIn;
170  iRecAft = iRecAftIn;
171  iEmtAft = iEmtAftIn;
172  iEmtAft2 = iEmtAft2In;
173 }
174 
175 void DireSplitInfo::clear() {
176  iRadBef = iRecBef = iRadAft = iRecAft = iEmtAft = iEmtAft2 = side
177  = type = system = systemRec = 0;
178  splittingSelName = "";
179  useForBranching = terminateEvolution = false;
180  for (int i= 0; i < int(particleSave.size()); ++i) particleSave[i].clear();
181  kinSave.clear();
182  extras.clear();
183 }
184 
185 void DireSplitInfo::list() {
186  cout << "List DireSplitInfo: "
187  << " name = " << splittingSelName << "\n"
188  << " [ id(radBef)= " << radBef()->id
189  << " id(recBef)= " << recBef()->id << " ] --> "
190  << " { id(radAft)= " << radAft()->id
191  << " id(emtAft)= " << emtAft()->id
192  << " id(emtAft2)= " << emtAft2()->id
193  << " id(recAft)= " << recAft()->id
194  << " } \n";
195  kinSave.list();
196  cout << "\n";
197 }
198 
199 //==========================================================================
200 
201 // Definition of color chain members.
202 
203 DireSingleColChain::DireSingleColChain(int iPos, const Event& state,
204  PartonSystems* partonSysPtr) {
205 
206  int colSign = (iPos > 0) ? 1 : -1;
207  iPos = abs(iPos);
208  int type = state[iPos].colType();
209  int iStart = iPos;
210  int iSys = partonSysPtr->getSystemOf(iPos,true);
211  int sizeSystem = partonSysPtr->sizeAll(iSys);
212  if (!state[iPos].isFinal() || colSign < 0) type *= -1;
213  addToChain(iPos, state);
214 
215  do {
216 
217  int icol = colEnd();
218  if (type < 0) icol = acolEnd();
219 
220  bool foundRad = false;
221  for ( int i = 0; i < sizeSystem; ++i) {
222  int j = partonSysPtr->getAll(iSys, i);
223  if ( j == iPos || state[j].colType() == 0) continue;
224  if (!state[j].isFinal()
225  && state[j].mother1() != 1
226  && state[j].mother1() != 2) continue;
227  int jcol = state[j].col();
228  int jacl = state[j].acol();
229  if ( type < 0) swap(jcol,jacl);
230  if ( !state[j].isFinal()) swap(jcol,jacl);
231  if ( jacl == icol ) {
232  iPos = j;
233  foundRad = true;
234  break;
235  }
236  }
237  // Found next color index in evolving system.
238  if (foundRad) addToChain( iPos, state);
239  // Try to use find next color index in system of mother (if
240  // different), and then exit.
241  else {
242  bool foundMotRad = false;
243 
244  // Try to find incoming particle in other systems, i.e. if the current
245  // system arose from a resonance decay.
246  int sizeSys = partonSysPtr->sizeSys();
247  int in1=0;
248  int iParentInOther = 0;
249  int nSys = partonSysPtr->sizeAll(iSys);
250  for (int iInSys = 0; iInSys < nSys; ++iInSys){
251  int iNow = partonSysPtr->getAll(iSys,iInSys);
252  for (int iOtherSys = 0; iOtherSys < sizeSys; ++iOtherSys){
253  if (iOtherSys == iSys) continue;
254  int nOtherSys = partonSysPtr->sizeAll(iOtherSys);
255  for (int iInOtherSys = 0; iInOtherSys < nOtherSys; ++iInOtherSys){
256  int iOtherNow = partonSysPtr->getAll(iOtherSys,iInOtherSys);
257  if (state[iNow].isAncestor(iOtherNow)) {
258  iParentInOther = iOtherNow;
259  }
260  }
261  }
262  }
263  in1 = iParentInOther;
264 
265  int jcol = state[in1].col();
266  int jacl = state[in1].acol();
267  if ( !state[in1].isFinal() || type < 0) swap(jcol,jacl);
268  if ( !state[in1].isFinal() && type < 0) swap(jcol,jacl);
269  if ( jacl == icol ) {
270  iPos = in1;
271  foundMotRad = true;
272  }
273 
274  if (foundMotRad) { addToChain( iPos, state); break;}
275  }
276 
277  } while ( abs(state[iPosEnd()].colType()) != 1 && iPosEnd() != iStart );
278  if (iPosEnd() == iStart) chain.pop_back();
279 
280 }
281 
282 void DireSingleColChain::addToChain(const int iPos, const Event& state){
283  int col = state[iPos].col();
284  int acl = state[iPos].acol();
285  original_chain.push_back( make_pair(iPos, make_pair(col, acl)));
286  if ( !state[iPos].isFinal()) swap(col,acl);
287  chain.push_back( make_pair(iPos, make_pair(col, acl)));
288 }
289 
290 bool DireSingleColChain::isInChain( int iPos) {
291  for (int i=0; i< size(); ++i)
292  if (chain[i].first == iPos) return true;
293  return false;
294 }
295 
296 int DireSingleColChain::posInChain( int iPos) {
297  for (int i=0; i< size(); ++i)
298  if (chain[i].first == iPos) return i;
299  return -1;
300 }
301 
302 bool DireSingleColChain::colInChain( int col) {
303  for (int i=0; i< size(); ++i)
304  if ( chain[i].second.first == col
305  || chain[i].second.second == col) return true;
306  return false;
307 }
308 
309 DireSingleColChain DireSingleColChain::chainFromCol(int iPos, int col,
310  int nSteps, const Event& state) {
311  DireSingleColChain ret;
312  int iSteps = 0;
313  int iPosInChain = posInChain(iPos);
314 
315  // For gluon, just match both colors.
316  if (state[iPos].id() == 21) {
317 
318  if (iPosInChain == 0) {
319 
320  ret.addToChain(chain[iPosInChain].first, state);
321  if ( iPosInChain+1 < size() && chain[iPosInChain+1].first > 0
322  && !ret.isInChain(chain[iPosInChain+1].first))
323  ret.addToChain(chain[iPosInChain+1].first, state);
324  if ( iPosInChain+2 < size() && chain[iPosInChain+2].first > 0
325  && !ret.isInChain(chain[iPosInChain+2].first))
326  ret.addToChain(chain[iPosInChain+2].first, state);
327 
328  } else if (iPosInChain == size()-1) {
329 
330  if ( iPosInChain-2 >= 0
331  && chain[iPosInChain-2].first > 0
332  && !ret.isInChain(chain[iPosInChain-2].first))
333  ret.addToChain(chain[iPosInChain-2].first, state);
334  if ( iPosInChain-1 >= 0 && iPosInChain-1 < size()
335  && chain[iPosInChain-1].first > 0
336  && !ret.isInChain(chain[iPosInChain-1].first))
337  ret.addToChain(chain[iPosInChain-1].first, state);
338  ret.addToChain(chain[iPosInChain].first, state);
339 
340  } else {
341 
342  if ( iPosInChain-1 >= 0 && iPosInChain-1 < size()
343  && chain[iPosInChain-1].first > 0
344  && !ret.isInChain(chain[iPosInChain-1].first))
345  ret.addToChain(chain[iPosInChain-1].first, state);
346  if ( iPosInChain >= 0 && iPosInChain < size()
347  && chain[iPosInChain].first > 0
348  && !ret.isInChain(chain[iPosInChain].first))
349  ret.addToChain(chain[iPosInChain].first, state);
350  if ( iPosInChain+1 < size() && chain[iPosInChain+1].first > 0
351  && !ret.isInChain(chain[iPosInChain+1].first))
352  ret.addToChain(chain[iPosInChain+1].first, state);
353  }
354  return ret;
355  }
356 
357  // Loop through, find color and attach subsequent particles to chain.
358  for (int i=0; i< size(); ++i) {
359  if ( iSteps == 0 && size() - 1 - i > nSteps
360  && chain[i].second.first != col
361  && chain[i].second.second != col) continue;
362  iSteps++;
363  if (chain[i].first > 0 && !ret.isInChain(chain[i].first))
364  ret.addToChain(chain[i].first, state);
365  if (iSteps > nSteps) break;
366  }
367 
368  // Done
369  return ret;
370 }
371 
372 string DireSingleColChain::listPos() const {
373  ostringstream os;
374  for (int i=0; i< size(); ++i) os << " " << chain[i].first;
375  return os.str();
376 }
377 
378 // List functions by N. Fischer.
379 void DireSingleColChain::print() const {
380  int i = 0;
381  int length = size();
382  int max = length;
383 
384  // first line: positions
385  for (i=0; i<length; i++) {
386  cout << setw(i == 0 ? 5 : 10) << chain[i].first;
387  }
388  cout << endl;
389  // second line: color-anticolor connections (horizontal)
390  i = 0;
391  max = (length%2 == 0 ? length : length-1);
392  while (i < max) {
393  if (i == 0) cout << " ";
394  if (i < max-1) cout << (i%2 == 0 ? " _____________" : " ");
395  i++;
396  }
397  cout << endl;
398  // third line: color-anticolor connections (vertical)
399  i = 0;
400  while (i < max) {
401  if (i == 0) cout << " ";
402  cout << "|";
403  if (i < max-1) cout << (i%2 == 0 ? " " : " ");
404  i++;
405  }
406  cout << endl;
407  // fourth line: colors and anticolors
408  for (i=0; i<length; i++) {
409  cout << setw(4) << chain[i].second.first;
410  cout << setw(4) << chain[i].second.second;
411  cout << " ";
412  }
413  cout << endl;
414  // fifth line: color-anticolor connections (vertical & horizontal)
415  i = 0;
416  max = (length%2 == 0 ? length-2 : length-1);
417  while (i < max) {
418  if (i == 0) cout << " ";
419  cout << "|";
420  if (i < max-1) cout << (i%2 == 0 ? "_____________" : " ");
421  i++;
422  }
423  cout << endl;
424  // sixth line: if first gluon is connected to last gluon
425  max = 10*(length-1)-5;
426  if (chain[0].second.second == chain[length-1].second.first
427  && chain[0].second.second != 0) {
428  cout << " |";
429  for (i=0; i<max; i++) cout << "_";
430  cout << "|";
431  }
432  cout << endl;
433 }
434 
435 // List functions by N. Fischer.
436 void DireSingleColChain::list () const {
437  if (size() > 0) cout << " ";
438  for (int i=0; i<size(); i++) {
439  cout << "[" << chain[i].second.second << "]";
440  cout << " " << chain[i].first << " ";
441  cout << "(" << chain[i].second.first << ")";
442  if (i < size()-1) cout << " --- ";
443  }
444  cout << endl;
445 }
446 
447 // List functions by N. Fischer.
448 string DireSingleColChain::list2 () const {
449  ostringstream os;
450  if (size() > 0) os << " ";
451  for (int i=0; i<size(); i++) {
452  os << "[" << chain[i].second.second << "]";
453  os << " " << chain[i].first << " ";
454  os << "(" << chain[i].second.first << ")";
455  if (i < size()-1) os << " --- ";
456  }
457  return os.str();
458 }
459 
460 //==========================================================================
461 
462 // Definition of DireColChain members
463 
464 DireSingleColChain DireColChains::chainOf (int iPos) {
465  for (int i=0; i< size(); ++i)
466  if ( chains[i].isInChain(iPos) ) return chains[i];
467  return DireSingleColChain();
468 }
469 
470 DireSingleColChain DireColChains::chainFromCol (int iPos, int col, int nSteps,
471  const Event& state) {
472  for (int i=0; i< size(); ++i) {
473  if ( chains[i].colInChain(col) ) {
474  return chains[i].chainFromCol(iPos, col, nSteps, state);}
475  }
476  return DireSingleColChain();
477 }
478 
479 int DireColChains::check(int iSys, const Event& state,
480  PartonSystems* partonSysPtr) {
481 
482  int sizeSystem = partonSysPtr->sizeAll(iSys);
483  int nFinal = 0;
484  for ( int i = 0; i < sizeSystem; ++i) {
485  int j = partonSysPtr->getAll(iSys, i);
486  if (!state[j].isFinal()) continue;
487  nFinal++;
488  if ( state[j].colType() == 0) continue;
489  if ( chainOf(j).size() < 2) return j;
490  }
491 
492  for ( int i = 0; i < sizeSystem; ++i) {
493  int j = partonSysPtr->getAll(iSys, i);
494  if ( state[j].colType() == 0) continue;
495  if ( state[j].mother1() != 1 && state[j].mother1() != 2) continue;
496  if ( nFinal > 0 && chainOf(j).size() < 2) return j;
497  }
498 
499  // Done.
500  return -1;
501 
502 }
503 
504 void DireColChains::list() {
505  cout << "\n --------- Begin DIRE Color Chain Listing -----------------"
506  << "--------------------------------------------------------------"
507  << "----------" << endl << endl;
508 
509  for (int i=0; i < size(); ++i){
510  cout << " Chain " << setw(4) << i << "\n" << endl;
511  chains[i].print();
512  if (i < size()-1)
513  cout << " **********************************************************"
514  << "***********************************************************"
515  << "**************" << endl;
516  }
517  // Done.
518  cout << " ---------- End DIRE Color Chain Listing -----------------"
519  << "--------------------------------------------------------------"
520  << "----------" << endl;
521 
522 }
523 
524 //==========================================================================
525 
526 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5