PYTHIA  8.312
VinciaHistory.h
1 // VinciaHistory.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 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 // Authors: Helen Brooks, Christian T Preuss
7 // This file contains the VinciaHistory class.
8 
9 #ifndef VINCIA_History_H
10 #define VINCIA_History_H
11 
12 #include "Pythia8/History.h"
13 #include "Pythia8/Info.h"
14 #include "Pythia8/Vincia.h"
15 #include "Pythia8/VinciaAntennaFunctions.h"
16 #include "Pythia8/VinciaMergingHooks.h"
17 
18 namespace Pythia8 {
19 
20 //==========================================================================
21 
22 // Convenient shorthand for storing ordered list of chains.
23 
24 struct PseudoChain {
25  // Ordered list of concatenated chains.
26  vector<int> chainlist;
27  // Index unique up to chain content (not ordering).
28  int index;
29  // Charge Index.
30  int cindex;
31  // Does any of the chains contain an initial state parton.
32  bool hasInitial;
33  // Flavour at start of first and end of last.
34  int flavStart;
35  int flavEnd;
36  // Charge.
37  int charge;
38 };
39 
40 //==========================================================================
41 
42 // A class for associating disconnected colour chains with either
43 // resonances or the hard QCD scattering. Once fixed, the sector
44 // history is deterministic.
45 
46 class ColourFlow {
47 
48  public:
49 
50  // Constructor.
51  ColourFlow() : nChains(0), nBeamChainsMin(0), nBeamChainsMax(0), nRes(0) {
52  for(int i=0; i<4; i++) {countChainsByChargeIndex[i]=0;
53  countResByChargeIndex[i]=0;}}
54 
55  // Methods.
56 
57  // Add or assign chains.
58  void addChain(int charge, int flavStart, int flavEnd, bool hasInitialIn);
59  void selectResChains(int index,int iorder, int id);
60  void selectBeamChains(int index,int iorder);
61 
62  // Initialise from hard process information.
63  bool initHard(map<int, map<int,int> >& countRes,
64  shared_ptr<VinciaMergingHooks> vinMergingHooksPtr);
65 
66  // Get information about current status of (pseudo)chains.
67  bool checkChains(int cIndex);
68  bool checkChains();
69  int getNChainsLeft();
70  int maxLength();
71  int minLength();
72 
73  // Print information about (pseudo)chains.
74  void print(bool printpsch = false);
75 
76  // Members.
77 
78  // Chains that arise from the decay of a resonance.
79  map<int, vector<PseudoChain> > resChains;
80 
81  // Remaining list of ordered chains after resonances stripped off.
82  vector<PseudoChain> beamChains;
83 
84  // Maps to all viable combinations of chains (of each charge)
85  // from an index which depends on the identities of chains in it
86  // - will be empty when resChains and beamChains set.
87  map< int , vector< PseudoChain > > pseudochains;
88 
89  // Easy way to look up all pseudochains which contain a given chain.
90  // Map from chain number to pseudochain index.
91  map<int,vector<int> > chainToIndices;
92 
93  // Map from chain to the flavour at the start or end of that chain.
94  map<int,int> chainStartToFlav;
95  map<int,int> chainEndToFlav;
96  // Keep track of chains which contain an initial state parton.
97  map<int,bool> hasInitial;
98  // Keep track of charge information of chains.
99  map<int,int> chainToCharge;
100 
101  private:
102 
103  // Add or assign chains.
104  void addChain(int oldIndex, int chainsIndex, int iChain,
105  vector<int> & newIndices);
106  void selectPseudochain(vector<int> & psch);
107  void selectChain(int iChain);
108 
109  // Create empty vectors in resChains, and count each type.
110  void addResonances(vector<int>& idsIn, map<int, map<int,int> >& idCounter,
111  int charge, bool fc);
112 
113  // Convert charge information to an index.
114  int getChargeIndex(int charge, bool fc);
115 
116  // Convert pseudochain to unique ID.
117  int getID(PseudoChain& psc) {
118  int id = 0;
119  int iChains = psc.chainlist.size()-1;
120  for (int iPwr(iChains); iPwr>=0; --iPwr)
121  id += (psc.chainlist.at(iPwr)+1)*pow(10,iChains-iPwr);
122  return id;
123  }
124 
125  // List of found pseudochains to not double count.
126  vector<int> pseudochainIDs;
127 
128  // Counters.
129  int nChains;
130  int nBeamChainsMin, nBeamChainsMax, nRes;
131  map<int,int> countChainsByChargeIndex;
132  map<int,int> countResByChargeIndex;
133 
134 };
135 
136 //==========================================================================
137 
138 // Class for a single step in the history of a process.
139 
140 class HistoryNode {
141 
142  public:
143 
144  // Constructors.
146  HistoryNode(Event& stateIn, vector< vector<int> > chainsIn,
147  double scaleIn) : HistoryNode() {
148  state = stateIn;
149  clusterableChains = chainsIn;
150  qEvolNow = scaleIn;
151  hasRes = false;
152  iRes = 0;
153  idRes = 0;
154  isInitPtr=false;
155  };
156 
157  // Methods.
158  void initPtr(VinciaCommon* vinComPtrIn, Resolution* resPtrIn,
159  AntennaSetFSR* antSetPtrIn) {
160  resPtr = resPtrIn;
161  vinComPtr = vinComPtrIn;
162  antSetFSRptr = antSetPtrIn;
163  isInitPtr=true;
164  nMinQQbar = 0;
165  }
166 
167  // Set clusterList.
168  int getNClusterings(shared_ptr<VinciaMergingHooks> vinMergingHooksPtr,
169  Logger* loggerPtr, int verboseIn);
170  void setClusterList(shared_ptr<VinciaMergingHooks> vinMergingHooksPtr,
171  Logger* loggerPtr, int verboseIn);
172 
173  // Perform the clusterings according to the resolution criterion.
174  bool cluster(HistoryNode& nodeClus,
175  Logger* loggerPtr, int verboseIn);
176 
177  // Get energy fractions (used in PDF ratios).
178  double xA() const {return 2. * state[3].e() / state[0].e();}
179  double xB() const {return 2. * state[4].e() / state[0].e();}
180 
181  // Get flavours of beams (used in PDF ratios).
182  int idA() const {return state[3].id();}
183  int idB() const {return state[4].id();}
184 
185  // Get colour types of beams (used in PDF ratios).
186  int colTypeA() const {return state[3].colType();}
187  int colTypeB() const {return state[4].colType();}
188 
189  // Get evolution scale (used in trial shower).
190  double getEvolNow() const {return qEvolNow;}
191 
192  // Setter methods.
193  void setEvolScale(double scaleIn) {qEvolNow = scaleIn;}
194 
195  // Current state.
197 
198  // Resonance info.
199  bool hasRes;
200  int iRes;
201  int idRes;
202 
203  // Minimal number of qqbar pairs.
205 
206  // List of unclusterable lists of colour-connected partons.
207  vector<vector<int>> clusterableChains;
208 
209  // Information corresponding to the last clustering.
211 
212  private:
213 
214  // Perform a clustering.
215  bool doClustering(VinciaClustering& clus, Event& clusEvent,
216  vector<vector<int>>& clusChains, Logger* loggerPtr, int verboseIn);
217 
218  // Methods to calculate resolution and evolution scales.
219  double calcResolution(VinciaClustering& clusIn) {
220  return resPtr->q2sector(clusIn);}
221  double calcEvolScale(VinciaClustering& clusIn) {
222  return resPtr->q2evol(clusIn);}
223 
224  // Members.
225 
226  // Vincia pointers.
227  Resolution* resPtr;
228  VinciaCommon* vinComPtr;
229  AntennaSetFSR* antSetFSRptr;
230 
231  bool isInitPtr;
232 
233  // The value of the evolution scale.
234  double qEvolNow;
235 
236  // List of next possible clusterings.
237  // Map is from corresponding resolution criterion.
238  map<double, VinciaClustering> clusterList;
239 
240 };
241 typedef map<int, vector<HistoryNode> > HistoryNodes;
242 
243 //==========================================================================
244 
245 // History class for the Vincia shower.
246 
248 
249 public:
250 
251  // Constructor.
252  VinciaHistory(Event &stateIn,
253  BeamParticle* beamAPtrIn, BeamParticle* beamBPtrIn,
254  MergingHooksPtr mergingHooksPtrIn,
255  PartonLevel* trialPartonLevelPtrIn,
256  ParticleData* particleDataPtrIn,
257  Info* infoPtrIn);
258 
259  // Check if a valid history was found.
260  bool isValid() {return foundValidHistory;}
261 
262  // Check if history failed merging scale cut.
263  bool isBelowMS() {return failedMSCut;}
264 
265  // Perform a trial shower and get the ckkwl weight.
266  double getWeightCKKWL();
267 
268  // What was the multiplicity of this event?
269  int getNClusterSteps();
270 
271  // Get the shower starting scale.
272  double getRestartScale();
273 
274  // Should we overwrite the original event?
275  // E.g. if an MPI was generated.
276  bool hasNewProcess() {return hasNewProcessSav;}
277  Event getNewProcess() {return newProcess;}
278  bool doAbort() {return aborted;}
279 
280 private:
281 
282  // Loop over colPerms and set historyBest.
283  void findBestHistory();
284 
285  // Find all viable colour orderings and return number.
286  unsigned int countPerms();
287 
288  // Fetch the colour chains from the hard event.
289  bool getColChains();
290 
291  // Find all selections of chains for colour singlets.
292  bool assignResChains(map<int, map<int,int>>& idCounter,
293  vector<ColourFlow>& flowsSoFar);
294 
295  // Find all selections of chains for colour singlets.
296  bool assignResFromEvent(map<int, map<int,int>>& idCounter,
297  vector<ColourFlow>& flowsSoFar);
298 
299  // Find all selections of chains for everything that is
300  // not a colour singlet resonance.
301  bool assignBeamChains(vector<ColourFlow>& flowsSoFar);
302 
303  // Make a single selection in all possible ways.
304  bool assignNext(vector<ColourFlow>& flowsSoFar, bool isRes = false,
305  int id = 0, int cIndex = 0);
306 
307  // Make a specific selection for a resonance.
308  bool assignThis(vector<ColourFlow>& flowsSoFar, int id, int cIndex,
309  vector<int>& chains);
310 
311  // Is this (completed flow) compatible with the number of
312  // resonances of each type?
313  bool check(ColourFlow& flow);
314 
315  // Construct history for a given colour permutation.
316  tuple<bool, double, HistoryNodes> findHistoryPerm(ColourFlow& flow);
317 
318  // Check if history failed merging scale cut.
319  bool checkMergingCut(HistoryNodes& history);
320 
321  // Initialise history nodes for each system.
322  HistoryNodes initHistoryNodes(ColourFlow& flow );
323 
324  // Translate abstract book-keeping of colourordering into
325  // systems of particles.
326  map<int, vector<vector<int>>> getSystems(ColourFlow& flow,
327  map<int, int>& sysToRes);
328 
329  // Decide if state is the Born topology.
330  bool isBorn(const HistoryNode& nodeIn, bool isRes);
331 
332  // Set up beams for given history node.
333  bool setupBeams(const HistoryNode* node, double scale2);
334 
335  // Calculate criterion for testing whether to keep history.
336  double calcME2guess(vector<HistoryNode>& history,bool isRes);
337  // Calculate ME for Born-level event.
338  double calcME2Born(const HistoryNode& bornNode, bool isRes);
339 
340  // Calculate the antenna function for a given clustering.
341  double calcAntFun(const VinciaClustering& clusNow);
342 
343  // Calculate PDF ratio to multiply CKKW-L weight.
344  double calcPDFRatio(const HistoryNode* nodeNow,
345  double pT2now, double pT2next);
346 
347  // Calculate alphaS ratio to multiply CKKW-L weight.
348  double calcAlphaSRatio(const HistoryNode& node);
349 
350  // Return the kinematic maximum for this event.
351  double getStartScale(Event& event, bool isRes);
352 
353  // Perform a trial branching and return scale.
354  double qNextTrial(double qStart, Event& evtIn);
355 
356  // Print colour chains.
357  void printChains();
358 
359  // Verbosity.
360  int verbose;
361 
362  // Beams -- for PDFs.
363  BeamParticle beamA, beamB;
364 
365  // MergingHooks.
366  shared_ptr<VinciaMergingHooks> vinMergingHooksPtr{};
367 
368  // PartonLevel pointer.
369  PartonLevel* trialPartonLevel{};
370 
371  // Particle data pointer.
372  ParticleData* particleDataPtr{};
373 
374  // Other Pythia pointers.
375  Info* infoPtr{};
376  Logger* loggerPtr{};
377 
378  // Vincia pointers.
379  shared_ptr<VinciaFSR> fsrShowerPtr{};
380  shared_ptr<VinciaISR> isrShowerPtr{};
381  MECs* mecsPtr{};
382  VinciaCommon* vinComPtr{};
383  Resolution* resPtr{};
384  AntennaSetFSR* antSetFSRptr{};
385 
386  // Check we found a valid history at all.
387  bool foundValidHistory;
388 
389  // Check we vetoed due failing the merging scale cut.
390  bool failedMSCut;
391 
392  // This is the best PS history so far.
393  HistoryNodes historyBest;
394 
395  // Criterion to minimise (best so far).
396  // Call calcME2guess(..) to evaluate.
397  double ME2guessBest;
398 
399  // Initial colour permutation, stored as
400  // list of sets of colour connected partons.
401  vector<vector<int>> colChainsSav;
402 
403  // Track if chain contains initial state quarks.
404  map<int, bool> chainHasInitial;
405 
406  // Keep track of resonances in the event record.
407  map<int, vector<int>> resIDToIndices;
408  map<int, vector<int>> resIndexToChains;
409 
410  // All colour flows compatible with Born process.
411  vector<ColourFlow> colPerms;
412 
413  // ME generated event.
414  Event state;
415 
416  // Number of quarks and gluon loops in event.
417  int nQSave, nGluonLoopsSave;
418 
419  // The merging scale and whether it is the evolution variable.
420  double qms;
421  bool msIsEvolVar;
422 
423  // The maximum multiplicity of our me-generator.
424  int nMax, nMaxRes;
425 
426  // Possible new hard process info (if MPI was generated).
427  bool hasNewProcessSav;
428  Event newProcess;
429  double newProcessScale;
430 
431  // Flag to signal if something went wrong.
432  bool aborted;
433 
434 };
435 
436 //==========================================================================
437 
438 } // end namespace Pythia8
439 
440 #endif // Pythia8_VinciaHistory_H
History class for the Vincia shower.
Definition: VinciaHistory.h:247
VinciaClustering lastClustering
Information corresponding to the last clustering.
Definition: VinciaHistory.h:210
void initPtr(VinciaCommon *vinComPtrIn, Resolution *resPtrIn, AntennaSetFSR *antSetPtrIn)
Methods.
Definition: VinciaHistory.h:158
Definition: Info.h:45
Event state
Current state.
Definition: VinciaHistory.h:196
The Event class holds all info on the generated event.
Definition: Event.h:453
Convenient shorthand for storing ordered list of chains.
Definition: VinciaHistory.h:24
Definition: BeamParticle.h:133
Simple struct to store information about a 3 -> 2 clustering.
Definition: VinciaCommon.h:278
bool hasRes
Resonance info.
Definition: VinciaHistory.h:199
vector< PseudoChain > beamChains
Remaining list of ordered chains after resonances stripped off.
Definition: VinciaHistory.h:82
map< int, int > chainStartToFlav
Map from chain to the flavour at the start or end of that chain.
Definition: VinciaHistory.h:94
int charge
Charge.
Definition: VinciaHistory.h:37
double xA() const
Get energy fractions (used in PDF ratios).
Definition: VinciaHistory.h:178
bool hasNewProcess()
Definition: VinciaHistory.h:276
bool isValid()
Check if a valid history was found.
Definition: VinciaHistory.h:260
Definition: Logger.h:23
map< int, vector< PseudoChain > > resChains
Members.
Definition: VinciaHistory.h:79
Definition: VinciaCommon.h:494
map< int, vector< int > > chainToIndices
Definition: VinciaHistory.h:91
Definition: PartonLevel.h:45
double getEvolNow() const
Get evolution scale (used in trial shower).
Definition: VinciaHistory.h:190
void setEvolScale(double scaleIn)
Setter methods.
Definition: VinciaHistory.h:193
int colTypeA() const
Get colour types of beams (used in PDF ratios).
Definition: VinciaHistory.h:186
bool isBelowMS()
Check if history failed merging scale cut.
Definition: VinciaHistory.h:263
ColourFlow()
Constructor.
Definition: VinciaHistory.h:51
Class for a single step in the history of a process.
Definition: VinciaHistory.h:140
Definition: VinciaHistory.h:46
A simple class for containing evolution variable definitions.
Definition: VinciaCommon.h:382
int cindex
Charge Index.
Definition: VinciaHistory.h:30
map< int, vector< PseudoChain > > pseudochains
Definition: VinciaHistory.h:87
int nMinQQbar
Minimal number of qqbar pairs.
Definition: VinciaHistory.h:204
The AntennaSetFSR class. Simple container of FF and RF antenna functions.
Definition: VinciaAntennaFunctions.h:1125
int idA() const
Get flavours of beams (used in PDF ratios).
Definition: VinciaHistory.h:182
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
vector< vector< int > > clusterableChains
List of unclusterable lists of colour-connected partons.
Definition: VinciaHistory.h:207
Definition: VinciaAntennaFunctions.h:1258
map< int, bool > hasInitial
Keep track of chains which contain an initial state parton.
Definition: VinciaHistory.h:97
HistoryNode()
Constructors.
Definition: VinciaHistory.h:145
bool hasInitial
Does any of the chains contain an initial state parton.
Definition: VinciaHistory.h:32
int flavStart
Flavour at start of first and end of last.
Definition: VinciaHistory.h:34
vector< int > chainlist
Ordered list of concatenated chains.
Definition: VinciaHistory.h:26
This class holds a map of all ParticleDataEntries.
Definition: ParticleData.h:422
map< int, int > chainToCharge
Keep track of charge information of chains.
Definition: VinciaHistory.h:99
int index
Index unique up to chain content (not ordering).
Definition: VinciaHistory.h:28