Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DSourceComboer.cc
Go to the documentation of this file.
4 
5 //Abandon all hope, ye who enter here.
6 //Seriously, it will take you at LEAST a month to understand this.
7 //If you really want to see what's going on, run with -PCOMBO:DEBUG_LEVEL=5000 and prepare to embrace the pain
8 
9 namespace DAnalysis
10 {
11 
12 /*
13 FAQ:
14 Q) If an event has the minimum # tracks, how can it fail to create combos for that event?
15 A) It may be that one of the tracks failed cuts for the PIDs that you need, but passed for others. Thus the total #tracks is OK.
16  Then, say you need 1 pi+ & 1 proton, and you detected 2 tracks. The other track may have passed dE/dx cuts for both proton & pi+, so it registers as both.
17  Also, one track could have both a positively & negatively charged hypothesis, and thus would count for both charges.
18 
19 FAQ:
20 Q) How can tracks have t1_detector() == SYS_NULL? I thought the PreSelect cuts were supposed to remove those?
21 A) Not exactly. If ANY of the hypos for a track has at least one hit in any detector, ALL hypos are saved.
22 
23 FAQ:
24 Q) How can I speed this up?
25 A) You can try reducing the #z-bins by increasing their widths. However, much sure you also increase the uncertainty on the timing & invariant-mass cuts for photons as well.
26 */
27 
28 /****************************************************** COMBOING STRATEGY ******************************************************
29 *
30 * Creating all possible combos can be very time- and memory-intensive if not done properly.
31 * For example, consider a 4pi0 analysis and 20 (N) reconstructed showers (it happens).
32 * If you make all possible pairs of photons (for pi0's), you get 19 + 18 + 17 + ... 1 = (N - 1)*N/2 = 190 pi0 combos.
33 * Now, consider that you have 4 pi0s: On the order of 190^4/16: On the order of a 80 million combos (although less once you guard against photon reuse)
34 *
35 * So, we must do everything we can to reduce the # of possible combos in ADVANCE of actually attempting to make them.
36 * And, we have to make sure we don't do anything twice (e.g. two different users have 4pi0s in their channel).
37 * The key to this being efficient (besides splitting the BCAL photons into vertex-z bins and placing timing cuts) is combo re-use.
38 *
39 * For example, suppose a channel needs 3 pi0s.
40 * First this will build all combos for 1 pi0, then all combos for 2 pi0s, then 3. Placing mass cuts along the way.
41 * The results after each of these steps is saved. That way, if someone then requests 2 pi0s, we merely have to return the results from the previous work.
42 * Also, if someone later requests 4pi0s, then we just take the 3pi0 results and expand them by 1 pi0.
43 * Or, if someone requests p3pi, we take the 1 pi0 combos and combine them with a proton, pi+, and pi-. Etc., etc.
44 *
45 * For more details on how this is done, see the comments in the Create_SourceCombos_Unknown function
46 * But ultimately, this results in a clusterfuck of recursive calls.
47 * Also, because of how the combo-info classes are structured (decaying PID NOT a member), you have be extremely careful not to get into an infinite loop.
48 * So, modify this code at your own peril. Just try not to take the rest of the collaboration down with you.
49 *
50 * Now, technically, when we construct combos for a (e.g.) pi0, we are saving 2 different results:
51 * The combos of 2 photons, and which of those combos survive the pi0 mass cut.
52 * That way, if later someone wants to build an eta, all we have to do is take 2-photon combos and place eta mass cuts.
53 *
54 * Combos are created on-demand, used immediately, and once they are cut the memory is recycled for the next combo in that event.
55 *
56 *
57 * The BCAL photons are evaluated in different vertex-z bins for calculating their kinematics (momentum & timing).
58 * This is because their kinematics have a strong dependence on vertex-z, while the FCAL showers do not (see above derivations).
59 * Whereas the FCAL photons have only a small dependence, so their kinematics are regardless of vertex-z.
60 * For more discussion the above, see the derivations in the DSourceComboTimeHandler and DSourceComboP4Handler classes.
61 *
62 *
63 *
64 *
65 *
66 * Note that combos are constructed separately for different beam bunches.
67 * This is because photons only survive their timing cuts for certain beam bunches.
68 * Comboing only within a given beam bunch reduces the #photons we need to combo, and is thus faster.
69 *
70 * When comboing, first all of the FCAL showers alone are used to build the requested combos.
71 * Then, the BCAL showers surviving the timing cuts within the input vertex-z bin are used to build the requested combos.
72 * Finally, combos are created using a mix of these BCAL & FCAL showers.
73 * The results from this comboing is saved for all cases, that way they can be easily retrieved and combined as needed for similar requests.
74 *
75 *
76 *******************************************************************************************************************************/
77 
78 /****************************************************** DESIGN MOTIVATION ******************************************************
79 *
80 *
81 *
82 * 1) Re-use comboing results between DReactions.
83 * If working on each DReaction individually, it is difficult (takes time & memory) to figure out what has already been done, and what to share
84 * So instead, first break down the DReactions to their combo-building components, and share those components.
85 * Then build combos out of the components, and distribute the results for each DReaction.
86 *
87 * 2) Reduce the time spent trying combos that we can know in advance won't work.
88 * We can do this by placing cuts IMMEDIATELY on:
89 * a) Time difference between charged tracks
90 * b) Time difference between photons and possible RF bunches (discussed more below).
91 * c) Invariant mass cuts for various decaying particles (e.g. pi0, eta, omega, phi, lambda, etc.)
92 * Also, when building combos of charged tracks, we could only loop over PIDs of the right type, rather than all hypotheses
93 *
94 * 3) The only way to do both 1) and 2) is to make the loose time & mass cuts reaction-independent.
95 * Users can always specify reaction-dependent tighter cuts later, but they cannot specify looser ones.
96 * However, these cuts should be tweakable on the command line in case someone wants to change them.
97 *
98 *******************************************************************************************************************************/
99 
100 
101 /*****
102  * COMBOING PHOTONS AND RF BUNCHES
103  *
104  * So, this is tricky.
105  * Start out by allowing ALL beam bunches, regardless of what the charged tracks want.
106  * Then, as each photon is chosen, reduce the set of possible photons to choose next: only those that agree on at least one RF bunch
107  * As combos are made, the valid RF bunches are saved along with the combo
108  * That way, as combos are combined with other combos/particles, we make sure that only valid possibilities are chosen.
109  *
110  * We can't start with those only valid for the charged tracks because:
111  * When we generate combos for a given info, we want to generate ALL combos at once.
112  * E.g. some charged tracks may want pi0s with beam bunch = 1, but another group might want pi0s with bunch 1 OR 2.
113  * Dealing with the overlap is a nightmare. This avoids the problem entirely.
114  *
115  * BEWARE: Massive-neutral-particle momentum depends on the RF bunch. So a cut on the invariant mass with a neutron is effectively a cut on the RF bunches
116  * Suppose: Sigma+ -> pi+ n
117  * You first generate combos for -> pi+ n, and save them for the use X -> pi+, n
118  * We then re-use the combos for the use Sigma+ -> pi+ n
119  * But then a cut on the Sigma+ mass reduces the #valid RF bunches. So now we need a new combo!
120  * We could decouple the RF bunches from the combo: e.g. save in map from combo_use -> rf bunches
121  * However, this would result in many duplicate entries: e.g. X -> 2g, pi0 -> 2g, eta -> 2g, etc.
122  * Users choosing final-state neutrons or KLongs is pretty rare compared to everything else: we are better off just creating new combos
123  *
124  * BEWARE: Massive-neutral-particle momentum depends on the RF bunch. So a cut on the invariant mass with a neutron is effectively a cut on the RF bunches.
125  * So we can't actually vote on RF bunches until we choose our massive-neutral particles!!!
126  */
127 
128 
130 {
131 //COMPARE:
132  //DEFINE DEFAULT dE/dx CUTS
133  //CDC Proton
134  ddEdxCuts_TF1FunctionStrings[Proton][SYS_CDC].first = "exp(-1.0*[0]*x + [1]) + [2]"; //low bound
135  ddEdxCuts_TF1Params[Proton][SYS_CDC].first = {4.0, 2.25, 1.0};
136  ddEdxCuts_TF1FunctionStrings[Proton][SYS_CDC].second = "[0]"; //high bound
137  ddEdxCuts_TF1Params[Proton][SYS_CDC].second = {9.9E9};
138 
139  //CDC Pi+
140  ddEdxCuts_TF1FunctionStrings[PiPlus][SYS_CDC].first = "[0]"; //low bound
141  ddEdxCuts_TF1Params[PiPlus][SYS_CDC].first = {-9.9E9};
142  ddEdxCuts_TF1FunctionStrings[PiPlus][SYS_CDC].second = "exp(-1.0*[0]*x + [1]) + [2]"; //high bound
143  ddEdxCuts_TF1Params[PiPlus][SYS_CDC].second = {7.0, 3.0, 6.2};
144 
145  //CDC K+
146  ddEdxCuts_TF1FunctionStrings[KPlus][SYS_CDC].first = "[0]"; //low bound
147  ddEdxCuts_TF1Params[KPlus][SYS_CDC].first = {-9.9E9};
148  ddEdxCuts_TF1FunctionStrings[KPlus][SYS_CDC].second = "exp(-1.0*[0]*x + [1]) + [2]"; //high bound
149  ddEdxCuts_TF1Params[KPlus][SYS_CDC].second = {7.0, 3.0, 6.2};
150 
151  //CDC e-
152  ddEdxCuts_TF1FunctionStrings[Electron][SYS_CDC].first = "[0]"; //low bound
153  ddEdxCuts_TF1Params[Electron][SYS_CDC].first = {-9.9E9};
154  ddEdxCuts_TF1FunctionStrings[Electron][SYS_CDC].second = "[0]"; //high bound
155  ddEdxCuts_TF1Params[Electron][SYS_CDC].second = {5.5};
156 
157  //pbar
160 
161  //Pi-
164 
165  //K-
168 
169  //e+
172 
173  //DEFINE DEFAULT E/p CUTS //vs p, cut away everything above if hadron, everything below if lepton
174  //e- FCAL
177 
178  //e- BCAL
181 
182  //e+
185 
186  /*
187  //mu-
188  dEOverPCuts_TF1FunctionStrings.emplace(MuonMinus, dEOverPCuts_TF1FunctionStrings[Electron]);
189  dEOverPCuts_TF1Params.emplace(MuonMinus, dEOverPCuts_TF1Params[Electron]);
190 
191  //mu+
192  dEOverPCuts_TF1FunctionStrings.emplace(MuonPlus, dEOverPCuts_TF1FunctionStrings[MuonMinus]);
193  dEOverPCuts_TF1Params.emplace(MuonPlus, dEOverPCuts_TF1Params[MuonMinus]);
194  */
195 }
196 
198 {
199  //PARAM EXAMPLES:
200  //COMBO_DEDXCUT:Low_14_1=0.75_0.5_1.0 //Cut protons (14) in the CDC (1) with the following parameters for the low-side cut
201  //COMBO_DEDXCUT:High_9_256_FUNC="[0] + [1]*x" //Cut pi-'s (9) in the SC (256) according to the functional form for the high-side cut //x = track momentum
202 
203  map<string, string> locParameterMap; //parameter key - filter, value
204  gPARMS->GetParameters(locParameterMap, "COMBO_DEDXCUT:"); //gets all parameters with this filter at the beginning of the key
205  for(auto locParamPair : locParameterMap)
206  {
207  if(dDebugLevel)
208  cout << "param pair: " << locParamPair.first << ", " << locParamPair.second << endl;
209 
210  //High or low cut?
211  auto locFirstUnderscoreIndex = locParamPair.first.find('_');
212  auto locSideString = locParamPair.first.substr(0, locFirstUnderscoreIndex);
213  auto locHighSideFlag = (locSideString == "Low") ? false : true;
214 
215  //Figure out which particle was specified
216  auto locSecondUnderscoreIndex = locParamPair.first.find('_', locFirstUnderscoreIndex + 1);
217  auto locParticleString = locParamPair.first.substr(locFirstUnderscoreIndex + 1, locSecondUnderscoreIndex);
218  istringstream locPIDtream(locParticleString);
219  int locPIDInt;
220  locPIDtream >> locPIDInt;
221  if(locPIDtream.fail())
222  continue;
223  Particle_t locPID = (Particle_t)locPIDInt;
224 
225  //Figure out which detector was specified
226  auto locFuncIndex = locParamPair.first.find("_FUNC");
227  auto locDetectorString = locParamPair.first.substr(locSecondUnderscoreIndex + 1, locFuncIndex);
228  istringstream locDetectorStream(locDetectorString);
229  int locSystemInt;
230  locDetectorStream >> locSystemInt;
231  if(locDetectorStream.fail())
232  continue;
233  DetectorSystem_t locSystem = (DetectorSystem_t)locSystemInt;
234 
235  if(dDebugLevel)
236  cout << "dE/dx cut: pid, detector, high-side flag = " << locPID << ", " << locSystem << ", " << locHighSideFlag << endl;
237 
238  //get the parameter, with hack so that don't get warning message about no default
239  string locKeyValue;
240  string locFullParamName = string("COMBO_DEDXCUT:") + locParamPair.first; //have to add back on the filter
241  gPARMS->SetDefaultParameter(locFullParamName, locKeyValue);
242 
243  //If functional form, save it and continue
244  if(locFuncIndex != string::npos)
245  {
246  if(!locHighSideFlag)
247  ddEdxCuts_TF1FunctionStrings[locPID][locSystem].first = locKeyValue;
248  else
249  ddEdxCuts_TF1FunctionStrings[locPID][locSystem].second = locKeyValue;
250  continue;
251  }
252 
253  //is cut parameters: extract and save
254  //CUT PARAMETERS SHOULD BE SEPARATED BY SPACES
255  //get rid of previous cut values
256  if(!locHighSideFlag)
257  ddEdxCuts_TF1Params[locPID][locSystem].first.clear();
258  else
259  ddEdxCuts_TF1Params[locPID][locSystem].second.clear();
260  while(true)
261  {
262  auto locUnderscoreIndex = locKeyValue.find('_');
263  auto locValueString = locKeyValue.substr(0, locUnderscoreIndex);
264 
265  istringstream locValuetream(locValueString);
266  double locParameter;
267  locValuetream >> locParameter;
268  if(locValuetream.fail())
269  continue; //must be for a different use
270  if(dDebugLevel)
271  cout << "param: " << locParameter << endl;
272 
273  //save locParameter and truncate locKeyValue (or break if done)
274  if(!locHighSideFlag)
275  ddEdxCuts_TF1Params[locPID][locSystem].first.push_back(locParameter);
276  else
277  ddEdxCuts_TF1Params[locPID][locSystem].second.push_back(locParameter);
278  if(locUnderscoreIndex == string::npos)
279  break;
280  locKeyValue = locKeyValue.substr(locUnderscoreIndex + 1);
281  }
282  }
283 }
284 
286 {
287  //PARAM EXAMPLES:
288  //COMBO_EOVERP:14_32_FUNC="[0] + [1]*x" //Cut protons (14) in the FCAL (32) according to the functional form //x = track momentum
289  //COMBO_EOVERP:14_32=0.75_0.5 //Cut protons (14) in the FCAL (32) with the following parameters
290 
291  map<string, string> locParameterMap; //parameter key - filter, value
292  gPARMS->GetParameters(locParameterMap, "COMBO_EOVERP:"); //gets all parameters with this filter at the beginning of the key
293  for(auto locParamPair : locParameterMap)
294  {
295  if(dDebugLevel)
296  cout << "param pair: " << locParamPair.first << ", " << locParamPair.second << endl;
297 
298  //Figure out which particle was specified
299  auto locUnderscoreIndex = locParamPair.first.find('_');
300  auto locParticleString = locParamPair.first.substr(0, locUnderscoreIndex);
301  istringstream locPIDtream(locParticleString);
302  int locPIDInt;
303  locPIDtream >> locPIDInt;
304  if(locPIDtream.fail())
305  continue;
306  Particle_t locPID = (Particle_t)locPIDInt;
307 
308  //Figure out which detector was specified
309  auto locFuncIndex = locParamPair.first.find("_FUNC");
310  auto locDetectorString = locParamPair.first.substr(locUnderscoreIndex + 1, locFuncIndex);
311  istringstream locDetectorStream(locDetectorString);
312  int locSystemInt;
313  locDetectorStream >> locSystemInt;
314  if(locDetectorStream.fail())
315  continue;
316  DetectorSystem_t locSystem = (DetectorSystem_t)locSystemInt;
317 
318  if(dDebugLevel)
319  cout << "E/p cut: pid, detector = " << locPID << ", " << locSystem << endl;
320 
321  //get the parameter, with hack so that don't get warning message about no default
322  string locKeyValue;
323  string locFullParamName = string("COMBO_EOVERP:") + locParamPair.first; //have to add back on the filter
324  gPARMS->SetDefaultParameter(locFullParamName, locKeyValue);
325 
326  //If functional form, save it and continue
327  if(locFuncIndex != string::npos)
328  {
329  dEOverPCuts_TF1FunctionStrings[locPID][locSystem] = locKeyValue;
330  continue;
331  }
332 
333  //is cut parameters: extract and save
334  //CUT PARAMETERS SHOULD BE SEPARATED BY SPACES
335  dEOverPCuts_TF1Params[locPID][locSystem].clear(); //get rid of previous cut values
336  while(true)
337  {
338  auto locUnderscoreIndex = locKeyValue.find('_');
339  auto locValueString = locKeyValue.substr(0, locUnderscoreIndex);
340 
341  istringstream locValuetream(locValueString);
342  double locParameter;
343  locValuetream >> locParameter;
344  if(locValuetream.fail())
345  continue; //must be for a different use
346  if(dDebugLevel)
347  cout << "param: " << locParameter << endl;
348 
349  //save locParameter and truncate locKeyValue (or break if done)
350  dEOverPCuts_TF1Params[locPID][locSystem].push_back(locParameter);
351  if(locUnderscoreIndex == string::npos)
352  break;
353  locKeyValue = locKeyValue.substr(locUnderscoreIndex + 1);
354  }
355  }
356 }
357 
359 {
360  //No idea why this lock is necessary, but it crashes without it. Stupid ROOT.
361  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
362 
363  //dE/dx
364  for(auto& locPIDPair : ddEdxCuts_TF1Params)
365  {
366  auto& locSystemMap = locPIDPair.second;
367  for(auto& locSystemPair : locSystemMap)
368  {
369  auto& locParamsPair = locSystemPair.second;
370  auto& locParamVector_Low = locParamsPair.first;
371  auto& locParamVector_High = locParamsPair.second;
372  if(locParamVector_Low.empty() || locParamVector_High.empty())
373  continue;
374 
375  //Get cut strings
376  if(ddEdxCuts_TF1FunctionStrings.find(locPIDPair.first) == ddEdxCuts_TF1FunctionStrings.end())
377  continue;
378  auto locSystemStringMap = ddEdxCuts_TF1FunctionStrings[locPIDPair.first];
379  if(locSystemStringMap.find(locSystemPair.first) == locSystemStringMap.end())
380  continue;
381  auto locCutFuncString_Low = locSystemStringMap[locSystemPair.first].first;
382  auto locCutFuncString_High = locSystemStringMap[locSystemPair.first].second;
383 
384  //Create TF1 low-side, Set cut values
385  auto locFunc_Low = new TF1("df_dEdxCut_Low", locCutFuncString_Low.c_str(), 0.0, 12.0);
386  if(dPrintCutFlag)
387  jout << "dE/dx Cut PID, System, low-side func form, params: " << ParticleType(locPIDPair.first) << ", " << SystemName(locSystemPair.first) << ", " << locCutFuncString_Low;
388  ddEdxCutMap[locPIDPair.first][locSystemPair.first].first = locFunc_Low;
389  for(size_t loc_i = 0; loc_i < locParamVector_Low.size(); ++loc_i)
390  {
391  locFunc_Low->SetParameter(loc_i, locParamVector_Low[loc_i]);
392  if(dPrintCutFlag)
393  jout << ", " << locParamVector_Low[loc_i];
394  }
395  if(dPrintCutFlag)
396  jout << endl;
397 
398  //Create TF1 high-side, Set cut values
399  auto locFunc_High = new TF1("df_dEdxCut_High", locCutFuncString_High.c_str(), 0.0, 12.0);
400  if(dPrintCutFlag)
401  jout << "dE/dx Cut PID, System, High-side func form, params: " << ParticleType(locPIDPair.first) << ", " << SystemName(locSystemPair.first) << ", " << locCutFuncString_High;
402  ddEdxCutMap[locPIDPair.first][locSystemPair.first].second = locFunc_High;
403  for(size_t loc_i = 0; loc_i < locParamVector_High.size(); ++loc_i)
404  {
405  locFunc_High->SetParameter(loc_i, locParamVector_High[loc_i]);
406  if(dPrintCutFlag)
407  jout << ", " << locParamVector_High[loc_i];
408  }
409  if(dPrintCutFlag)
410  jout << endl;
411  }
412  }
413 
414  //E/p
415  for(auto& locPIDPair : dEOverPCuts_TF1Params)
416  {
417  auto& locSystemMap = locPIDPair.second;
418  for(auto& locSystemPair : locSystemMap)
419  {
420  auto& locParamVector = locSystemPair.second;
421 
422  //Get cut strings
423  if(dEOverPCuts_TF1FunctionStrings.find(locPIDPair.first) == dEOverPCuts_TF1FunctionStrings.end())
424  continue;
425  auto locSystemStringMap = dEOverPCuts_TF1FunctionStrings[locPIDPair.first];
426  if(locSystemStringMap.find(locSystemPair.first) == locSystemStringMap.end())
427  continue;
428  auto locCutFuncString = locSystemStringMap[locSystemPair.first];
429 
430  //Create TF1, Set cut values
431  auto locFunc = new TF1("df_EOverPCut", locCutFuncString.c_str(), 0.0, 12.0);
432  if(dPrintCutFlag)
433  jout << "E/p Cut PID, System, func form, params: " << ParticleType(locPIDPair.first) << ", " << SystemName(locSystemPair.first) << ", " << locCutFuncString;
434  dEOverPCutMap[locPIDPair.first][locSystemPair.first] = locFunc;
435  for(size_t loc_i = 0; loc_i < locParamVector.size(); ++loc_i)
436  {
437  locFunc->SetParameter(loc_i, locParamVector[loc_i]);
438  if(dPrintCutFlag)
439  jout << ", " << locParamVector[loc_i];
440  }
441  if(dPrintCutFlag)
442  jout << endl;
443  }
444  }
445 
446  japp->RootUnLock(); //RELEASE ROOT LOCK!!
447 }
448 
449 /********************************************************************* CONSTRUCTOR **********************************************************************/
450 
451 DSourceComboer::DSourceComboer(JEventLoop* locEventLoop)
452 {
453  dResourcePool_SourceCombo.Set_ControlParams(100, 50, 1000, 20000, 0);
455  dCreatedCombos.reserve(100000);
456  dCreatedComboVectors.reserve(1000);
457 
458  //Get preselect tag, debug level
459  gPARMS->SetDefaultParameter("COMBO:SHOWER_SELECT_TAG", dShowerSelectionTag);
460  gPARMS->SetDefaultParameter("COMBO:DEBUG_LEVEL", dDebugLevel);
461  gPARMS->SetDefaultParameter("COMBO:PRINT_CUTS", dPrintCutFlag);
462  gPARMS->SetDefaultParameter("COMBO:MAX_NEUTRALS", dMaxNumNeutrals);
463 
464  //GET THE GEOMETRY
465  DApplication* locApplication = dynamic_cast<DApplication*>(locEventLoop->GetJApplication());
466  DGeometry* locGeometry = locApplication->GetDGeometry(locEventLoop->GetJEvent().GetRunNumber());
467 
468  //TARGET INFORMATION
469  double locTargetCenterZ = 65.0;
470  locGeometry->GetTargetZ(locTargetCenterZ);
471  dTargetCenter.SetXYZ(0.0, 0.0, locTargetCenterZ);
472 
473  //SETUP CUTS
478 
479  //GET THE REACTIONS
480  auto locReactions = DAnalysis::Get_Reactions(locEventLoop);
481 
482  //CREATE DSourceComboINFO'S
483  vector<const DReactionVertexInfo*> locVertexInfos;
484  locEventLoop->Get(locVertexInfos);
485  for(const auto& locVertexInfo : locVertexInfos)
486  Create_SourceComboInfos(locVertexInfo);
487 
488  //TRANSFER INFOS FROM SET TO VECTOR
489  dSourceComboInfos.reserve(dSourceComboInfoSet.size());
490  std::copy(dSourceComboInfoSet.begin(), dSourceComboInfoSet.end(), std::back_inserter(dSourceComboInfos));
491  dSourceComboInfoSet.clear(); //free up the memory
492 
493  //CREATE HANDLERS
501 
502  //save rf bunch cuts
503  if(gPARMS->Exists("COMBO:NUM_PLUSMINUS_RF_BUNCHES"))
504  {
505  size_t locNumPlusMinusRFBunches;
506  gPARMS->GetParameter("COMBO:NUM_PLUSMINUS_RF_BUNCHES", locNumPlusMinusRFBunches);
507  for(const auto& locReaction : locReactions)
508  dRFBunchCutsByReaction.emplace(locReaction, locNumPlusMinusRFBunches);
509  }
510  else //by reaction
511  {
512  for(const auto& locReaction : locReactions)
513  {
514  auto locNumBunches = locReaction->Get_NumPlusMinusRFBunches();
515  pair<bool, double> locMaxPhotonRFDeltaT = locReaction->Get_MaxPhotonRFDeltaT(); //DEPRECATED!!!
516  if(locMaxPhotonRFDeltaT.first)
517  locNumBunches = size_t(locMaxPhotonRFDeltaT.second/dSourceComboTimeHandler->Get_BeamBunchPeriod() - 0.499999999);
518  dRFBunchCutsByReaction.emplace(locReaction, locNumBunches);
519  }
520  }
521 
522  //save max bunch cuts
523  for(const auto& locVertexInfo : locVertexInfos)
524  {
525  dMaxRFBunchCuts.emplace(locVertexInfo, 0);
526  for(const auto& locReaction : locVertexInfo->Get_Reactions())
527  {
528  if(dRFBunchCutsByReaction[locReaction] > dMaxRFBunchCuts[locVertexInfo])
529  dMaxRFBunchCuts[locVertexInfo] = dRFBunchCutsByReaction[locReaction];
530  }
531  }
532 
533  //Make sure this matches DConstructionStage!!!
534  vector<string> locBuildStages_Event = {"Input", "Min # Particles", "Max # Particles", "In Skim", "Charged Combos", "Charged RF Bunch", "Full Combos", "Neutral RF Bunch",
535  "No-Vertex RF Bunch", "Heavy-Neutral IM", "Beam Combos", "MM/Dangling-Vertex Timing", "MM/Dangling-Vertex IM Cuts", "Accurate-Photon IM", "Reaction Beam-RF Cuts", "Missing Mass"};
536  vector<string> locBuildStages_Combo{"In Skim Events"}; //has same bin content as "In Skim"
537  locBuildStages_Combo.insert(locBuildStages_Combo.end(), locBuildStages_Event.begin() + 4, locBuildStages_Event.end());
538 
539  //initialize success tracking
540  for(auto locReaction : locReactions)
541  {
542  for(auto locStage = static_cast<DConstructionStageType>(DConstructionStage::Min_Particles); locStage <= static_cast<DConstructionStageType>(DConstructionStage::Missing_Mass); ++locStage)
543  dNumCombosSurvivedStageTracker[locReaction][static_cast<DConstructionStage>(locStage)] = 0;
544  }
545 
546  //Setup hists
547  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!
548  {
549  vector<DetectorSystem_t> locdEdxSystems {SYS_CDC, SYS_FDC, SYS_START, SYS_TOF};
550  vector<Particle_t> locPIDs {Electron, Positron, MuonPlus, MuonMinus, PiPlus, PiMinus, KPlus, KMinus, Proton, AntiProton};
551  vector<DetectorSystem_t> locEOverPSystems {SYS_BCAL, SYS_FCAL};
552 
553  //get and change to the base (file/global) directory
554  TDirectory* locCurrentDir = gDirectory;
555 
556  string locDirName = "Independent";
557  TDirectoryFile* locDirectoryFile = static_cast<TDirectoryFile*>(gDirectory->GetDirectory(locDirName.c_str()));
558  if(locDirectoryFile == NULL)
559  locDirectoryFile = new TDirectoryFile(locDirName.c_str(), locDirName.c_str());
560  locDirectoryFile->cd();
561 
562  locDirName = "Combo_Construction";
563  locDirectoryFile = static_cast<TDirectoryFile*>(gDirectory->GetDirectory(locDirName.c_str()));
564  if(locDirectoryFile == NULL)
565  locDirectoryFile = new TDirectoryFile(locDirName.c_str(), locDirName.c_str());
566  locDirectoryFile->cd();
567 
568  locDirName = "PID";
569  locDirectoryFile = static_cast<TDirectoryFile*>(gDirectory->GetDirectory(locDirName.c_str()));
570  if(locDirectoryFile == NULL)
571  locDirectoryFile = new TDirectoryFile(locDirName.c_str(), locDirName.c_str());
572  locDirectoryFile->cd();
573 
574  for(auto& locPID : locPIDs)
575  {
576  locDirName = ParticleType(locPID);
577  locDirectoryFile = static_cast<TDirectoryFile*>(gDirectory->GetDirectory(locDirName.c_str()));
578  if(locDirectoryFile == NULL)
579  locDirectoryFile = new TDirectoryFile(locDirName.c_str(), locDirName.c_str());
580  locDirectoryFile->cd();
581 
582  for(auto& locSystem : locdEdxSystems)
583  {
584  string locHistName = string("dEdxVsP_") + string(SystemName(locSystem));
585  auto locHist = gDirectory->Get(locHistName.c_str());
586  if(locHist == nullptr)
587  {
588  string locUnits = ((locSystem == SYS_CDC) || (locSystem == SYS_FDC)) ? "(keV/cm)" : "(MeV/cm)";
589  string locHistTitle = ParticleName_ROOT(locPID) + string(", ") + string(SystemName(locSystem)) + string(";p (GeV/c);dE/dX ") + locUnits;
590  dHistMap_dEdx[locPID][locSystem] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), 400, 0.0, 12.0, 400, 0.0, 25.0);
591  }
592  else
593  dHistMap_dEdx[locPID][locSystem] = static_cast<TH2*>(locHist);
594  }
595 
596  for(auto& locSystem : locEOverPSystems)
597  {
598  string locHistName = string("EOverP_") + string(SystemName(locSystem));
599  auto locHist = gDirectory->Get(locHistName.c_str());
600  if(locHist == nullptr)
601  {
602  string locHistTitle = ParticleName_ROOT(locPID) + string(", ") + string(SystemName(locSystem)) + string(";p (GeV/c);E_{Shower}/p_{Track} (c)");
603  dHistMap_EOverP[locPID][locSystem] = new TH2I(locHistName.c_str(), locHistTitle.c_str(), 400, 0.0, 12.0, 400, 0.0, 4.0);
604  }
605  else
606  dHistMap_EOverP[locPID][locSystem] = static_cast<TH2*>(locHist);
607  }
608  gDirectory->cd("..");
609  }
610  locCurrentDir->cd();
611 
612  //construction stage tracking
613  for(auto locReaction : locReactions)
614  {
615  string locReactionName = locReaction->Get_ReactionName();
616 
617  locDirName = locReactionName;
618  locDirectoryFile = static_cast<TDirectoryFile*>(gDirectory->GetDirectory(locDirName.c_str()));
619  if(locDirectoryFile == NULL)
620  locDirectoryFile = new TDirectoryFile(locDirName.c_str(), locDirName.c_str());
621  locDirectoryFile->cd();
622 
623  string locHistName = "ComboConstruction_NumEventsSurvived";
624  auto locHist = gDirectory->Get(locHistName.c_str());
625  if(locHist == nullptr)
626  {
627  string locHistTitle = locReactionName + string(";;# Events Survived Stage");
628  dNumEventsSurvivedStageMap[locReaction] = new TH1D(locHistName.c_str(), locHistTitle.c_str(), locBuildStages_Event.size(), -0.5, locBuildStages_Event.size() - 0.5);
629  for(size_t loc_i = 0; loc_i < locBuildStages_Event.size(); ++loc_i)
630  dNumEventsSurvivedStageMap[locReaction]->GetXaxis()->SetBinLabel(loc_i + 1, locBuildStages_Event[loc_i].c_str());
631  }
632  else
633  dNumEventsSurvivedStageMap[locReaction] = static_cast<TH1*>(locHist);
634 
635  locHistName = "ComboConstruction_NumCombosSurvived";
636  locHist = gDirectory->Get(locHistName.c_str());
637  if(locHist == nullptr)
638  {
639  string locHistTitle = locReactionName + string(";;# Combos Survived Stage");
640  dNumCombosSurvivedStageMap[locReaction] = new TH1D(locHistName.c_str(), locHistTitle.c_str(), locBuildStages_Combo.size(), -0.5, locBuildStages_Combo.size() - 0.5);
641  for(size_t loc_i = 0; loc_i < locBuildStages_Combo.size(); ++loc_i)
642  dNumCombosSurvivedStageMap[locReaction]->GetXaxis()->SetBinLabel(loc_i + 1, locBuildStages_Combo[loc_i].c_str());
643  }
644  else
645  dNumCombosSurvivedStageMap[locReaction] = static_cast<TH1*>(locHist);
646 
647  locHistName = "ComboConstruction_NumCombosSurvived2D";
648  locHist = gDirectory->Get(locHistName.c_str());
649  if(locHist == nullptr)
650  {
651  string locHistTitle = locReactionName + string(";;# Combos Survived Stage");
652  dNumCombosSurvivedStage2DMap[locReaction] = new TH2D(locHistName.c_str(), locHistTitle.c_str(), locBuildStages_Combo.size(), -0.5, locBuildStages_Combo.size() - 0.5, 1000, 0, 1000);
653  for(size_t loc_i = 0; loc_i < locBuildStages_Combo.size(); ++loc_i)
654  dNumCombosSurvivedStage2DMap[locReaction]->GetXaxis()->SetBinLabel(loc_i + 1, locBuildStages_Combo[loc_i].c_str());
655  }
656  else
657  dNumCombosSurvivedStage2DMap[locReaction] = static_cast<TH2*>(locHist);
658 
659  gDirectory->cd("..");
660  }
661  locCurrentDir->cd();
662  }
663  japp->RootUnLock(); //RELEASE ROOT LOCK!!
664 }
665 
667 {
668  auto locNumPreComboStages = 3; //"In Skim" will be first for #combos
669  japp->WriteLock("DSourceComboer_Survival");
670  {
671  for(auto& locReactionPair : dNumCombosSurvivedStageTracker)
672  {
673  auto& locReaction = locReactionPair.first;
674  for(auto& locStagePair : locReactionPair.second)
675  {
676  auto locNumCombos = locStagePair.second;
677  if(locNumCombos == 0)
678  break;
679 
680  auto locStageIndex = static_cast<std::underlying_type<DConstructionStage>::type>(locStagePair.first);
681  dNumEventsSurvivedStageMap[locReaction]->Fill(locStageIndex);
682  if(locStageIndex < locNumPreComboStages)
683  continue;
684 
685  //fill combo hists
686  auto locBin = locStageIndex - locNumPreComboStages + 1;
687  auto locBinContent = dNumCombosSurvivedStageMap[locReaction]->GetBinContent(locBin) + locNumCombos;
688  dNumCombosSurvivedStageMap[locReaction]->SetBinContent(locBin, locBinContent);
689 
690  if(dNumCombosSurvivedStage2DMap[locReaction]->GetYaxis()->FindBin(locNumCombos) <= dNumCombosSurvivedStage2DMap[locReaction]->GetNbinsY())
691  dNumCombosSurvivedStage2DMap[locReaction]->Fill(locBin, locNumCombos);
692  }
693  }
694  }
695  japp->Unlock("DSourceComboer_Survival");
696 
697  //Reset for next event
698  for(auto& locReactionPair : dNumCombosSurvivedStageTracker)
699  {
700  for(auto& locStagePair : locReactionPair.second)
701  locStagePair.second = 0;
702  }
703 }
704 
705 /******************************************************************* CREATE DSOURCOMBOINFO'S ********************************************************************/
706 
708 {
709  //FULL combo use: Segregate each step into (up to) 3 combos: a fully charged, a fully neutral, and a mixed
710  //that way we will combo each separately before combining them horizontally: maximum re-use, especially of time-intensive neutral comboing
711  //However, an exception: if a any-# of a single neutral PID (e.g. pi0, n, or g), promote it to the level where the charged/neutral/mixed are combined
712  //Charged is similar, but not the same: if a single DECAYING-to-charged particle, promote it as well
713  //Not so for a single detected charged particle though: We want to keep charged separate because that's what defines the vertices: Easier lookup
714 
715  /*
716  * suppose reaction is 0) g, p -> omega, p
717  * 1) omega -> 3pi
718  * 2) pi0 -> 2g
719  *
720  * It will have uses/infos like:
721  * 0: X -> A, 1 (mixed + charged)
722  * A: X -> p (charged)
723  * 1: omega -> B, 2 (mixed)
724  * B: X -> pi+, pi- (charged)
725  * 2: pi0 -> 2g (neutral)
726  */
727 
728  /*
729  * suppose reaction is 0) g, p -> K0, Sigma+
730  * 1) K0 -> 3pi
731  * 2) pi0 -> 2g
732  * 3) Sigma+ -> pi+, n
733  *
734  * It will have uses/infos like:
735  * 0: X -> A, 1, 3 (mixed -> charged, mixed, mixed)
736  * A: X -> p (charged)
737  * 1: K0 -> B, 2 (mixed -> charged, neutral)
738  * B: X -> pi+, pi- (charged)
739  * 2: pi0 -> 2g (neutral)
740  * 3: Sigma+ -> C, n (mixed -> charged, n)
741  * C: X -> pi+ (charged)
742  */
743 
744  if(dDebugLevel > 0)
745  cout << "CREATING DSourceComboInfo OBJECTS FOR DREACTION " << locReactionVertexInfo->Get_Reaction()->Get_ReactionName() << endl;
746 
747  //We will register what steps these combos are created for
748  map<size_t, DSourceComboUse> locStepComboUseMap; //size_t = step index
749 
750  //loop over steps in reverse order
751  auto locReaction = locReactionVertexInfo->Get_Reaction();
752  auto locReactionSteps = locReaction->Get_ReactionSteps();
753  for(auto locStepIterator = locReactionSteps.rbegin(); locStepIterator != locReactionSteps.rend(); ++locStepIterator)
754  {
755  auto locStep = *locStepIterator;
756  auto locStepIndex = locReaction->Get_NumReactionSteps() - std::distance(locReactionSteps.rbegin(), locStepIterator) - 1;
757  if(dDebugLevel >= 5)
758  cout << "Step index " << locStepIndex << endl;
759 
760  //create combo uses for all charged, all neutral, then for any mixed decays
761  map<Particle_t, unsigned char> locChargedParticleMap = Build_ParticleMap(locReaction, locStepIndex, d_Charged);
762  map<Particle_t, unsigned char> locNeutralParticleMap = Build_ParticleMap(locReaction, locStepIndex, d_Neutral);
763 
764  //get combo infos for final-state decaying particles //if not present, ignore parent
765  auto locFinalStateDecayingComboUsesPair = Get_FinalStateDecayingComboUses(locReaction, locStepIndex, locStepComboUseMap);
766  auto locIncludeParentFlag = locFinalStateDecayingComboUsesPair.first;
767  auto& locFurtherDecays = locFinalStateDecayingComboUsesPair.second;
768 
769  //split up further-decays into all-charged, all-neutral, and mixed
770  map<DSourceComboUse, unsigned char> locFurtherDecays_Charged, locFurtherDecays_Neutral, locFurtherDecays_Mixed;
771  for(const auto& locDecayPair : locFurtherDecays)
772  {
773  auto locChargeContent = dComboInfoChargeContent[std::get<2>(locDecayPair.first)];
774  if(locChargeContent == d_Charged)
775  locFurtherDecays_Charged.emplace(locDecayPair);
776  else if(locChargeContent == d_Neutral)
777  locFurtherDecays_Neutral.emplace(locDecayPair);
778  else
779  locFurtherDecays_Mixed.emplace(locDecayPair);
780  }
781 
782  //exclude parent if production
783  if((locStepIndex == 0) && DAnalysis::Get_IsFirstStepBeam(locReaction)) //decay
784  locIncludeParentFlag = false;
785 
786  //create combo uses for each case
787  auto locInitPID = locIncludeParentFlag ? locStep->Get_InitialPID() : Unknown;
788  bool locNoChargedFlag = (locChargedParticleMap.empty() && locFurtherDecays_Charged.empty());
789  bool locNoNeutralFlag = (locNeutralParticleMap.empty() && locFurtherDecays_Neutral.empty());
790 
791  //determine if we need to subtract a target particle when calculating the invariant mass (e.g. rescattering)
792  auto locTargetToInclude = (locStepIndex != 0) ? locStep->Get_TargetPID() : Unknown;
793 
794  //determine if there is a missing decay product, such that we can't do invariant mass cuts
795  bool locMissingDecayProductFlag = false;
796  if((locStepIndex != 0) || !DAnalysis::Get_IsFirstStepBeam(locReaction)) //decay
797  locMissingDecayProductFlag = DAnalysis::Check_IfMissingDecayProduct(locReaction, locStepIndex);
798 
799  if(dDebugLevel >= 5)
800  cout << "locIncludeParentFlag, init pid, missing-product flag, to-include target pid: " << locIncludeParentFlag << ", " << locInitPID << ", " << locMissingDecayProductFlag << ", " << locTargetToInclude << endl;
801 
802  //default to unknown use
803  DSourceComboUse locPrimaryComboUse(Unknown, DSourceComboInfo::Get_VertexZIndex_ZIndependent(), nullptr, false, Unknown);
804  if(locNoChargedFlag && locNoNeutralFlag) //only mixed
805  locPrimaryComboUse = Make_ComboUse(locInitPID, {}, locFurtherDecays_Mixed, locMissingDecayProductFlag, locTargetToInclude);
806  else if(locNoNeutralFlag && locFurtherDecays_Mixed.empty()) //only charged
807  locPrimaryComboUse = Make_ComboUse(locInitPID, locChargedParticleMap, locFurtherDecays_Charged, locMissingDecayProductFlag, locTargetToInclude);
808  else if(locNoChargedFlag && locFurtherDecays_Mixed.empty()) //only neutral
809  locPrimaryComboUse = Make_ComboUse(locInitPID, locNeutralParticleMap, locFurtherDecays_Neutral, locMissingDecayProductFlag, locTargetToInclude);
810  else //some combination
811  {
812  auto locFurtherDecays_All = locFurtherDecays_Mixed;
813  map<Particle_t, unsigned char> locParticleMap_All = {};
814  //create a combo for each charged group, with init pid = unknown
815  if(!locNoChargedFlag)
816  {
817  //if lone Charged decaying particle, promote to be parallel with mixed
818  if(locChargedParticleMap.empty() && (locFurtherDecays_Charged.size() == 1) && (locFurtherDecays_Charged.begin()->second == 1))
819  locFurtherDecays_All.emplace(locFurtherDecays_Charged.begin()->first, 1);
820  else //multiple Charged decaying particles, group together separately (own use)
821  {
822  auto locComboUse_Charged = Make_ComboUse(Unknown, locChargedParticleMap, locFurtherDecays_Charged, false, Unknown);
823  locFurtherDecays_All.emplace(locComboUse_Charged, 1);
824  }
825  }
826  if(!locNoNeutralFlag)
827  {
828  //if lone neutral PID, promote to be parallel with mixed
829  if(locNeutralParticleMap.empty() && (locFurtherDecays_Neutral.size() == 1))
830  locFurtherDecays_All.emplace(locFurtherDecays_Neutral.begin()->first, locFurtherDecays_Neutral.begin()->second); //decaying
831  else if(locFurtherDecays_Neutral.empty() && (locNeutralParticleMap.size() == 1))
832  locParticleMap_All.emplace(locNeutralParticleMap.begin()->first, locNeutralParticleMap.begin()->second); //detected
833  else //multiple neutral particles, group together separately (own use)
834  {
835  auto locComboUse_Neutral = Make_ComboUse(Unknown, locNeutralParticleMap, locFurtherDecays_Neutral, false, Unknown);
836  locFurtherDecays_All.emplace(locComboUse_Neutral, 1);
837  }
838  }
839 
840  locPrimaryComboUse = Make_ComboUse(locInitPID, locParticleMap_All, locFurtherDecays_All, locMissingDecayProductFlag, locTargetToInclude);
841  }
842 
843  locStepComboUseMap.emplace(locStepIndex, locPrimaryComboUse);
844  }
845 
846  //Register the results!!
847  for(const auto& locStepVertexInfo : locReactionVertexInfo->Get_StepVertexInfos())
848  dSourceComboUseReactionMap.emplace(locStepVertexInfo, locStepComboUseMap[locStepVertexInfo->Get_StepIndices().front()]);
849  for(const auto& locUseStepPair : locStepComboUseMap)
850  dSourceComboInfoStepMap.emplace(std::make_pair(locReactionVertexInfo->Get_StepVertexInfo(locUseStepPair.first), locUseStepPair.second), locUseStepPair.first);
851  for(auto locTempReaction : locReactionVertexInfo->Get_Reactions())
852  dSourceComboUseReactionStepMap.emplace(locTempReaction, locStepComboUseMap);
853 
854  if(dDebugLevel > 0)
855  cout << "DSourceComboInfo OBJECTS CREATED" << endl;
856 }
857 
858 pair<bool, map<DSourceComboUse, unsigned char>> DSourceComboer::Get_FinalStateDecayingComboUses(const DReaction* locReaction, size_t locStepIndex, const map<size_t, DSourceComboUse>& locStepComboUseMap) const
859 {
860  //get combo infos for final-state decaying particles //if one is not present, ignore parent
861  auto locIncludeParentFlag = true; //unless changed below
862  map<DSourceComboUse, unsigned char> locFurtherDecays;
863  auto locStep = locReaction->Get_ReactionStep(locStepIndex);
864  for(size_t loc_i = 0; loc_i < locStep->Get_NumFinalPIDs(); ++loc_i)
865  {
866  int locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, locStepIndex, loc_i);
867  if(locDecayStepIndex < 0)
868  continue;
869 
870  auto locUseIterator = locStepComboUseMap.find(size_t(locDecayStepIndex));
871  if(locUseIterator == locStepComboUseMap.end())
872  locIncludeParentFlag = false;
873  else
874  {
875  //save decay
876  auto& locSourceComboUse = locUseIterator->second;
877  auto locDecayIterator = locFurtherDecays.find(locSourceComboUse);
878  if(locDecayIterator == locFurtherDecays.end())
879  locFurtherDecays.emplace(locSourceComboUse, 1);
880  else
881  ++(locDecayIterator->second);
882  }
883  }
884 
885  return std::make_pair(locIncludeParentFlag, locFurtherDecays);
886 }
887 
888 map<Particle_t, unsigned char> DSourceComboer::Build_ParticleMap(const DReaction* locReaction, size_t locStepIndex, Charge_t locCharge) const
889 {
890  //build map of charged particles
891  map<Particle_t, unsigned char> locNumParticles;
892  auto locParticles = locReaction->Get_FinalPIDs(locStepIndex, false, false, locCharge, true); //no missing or decaying, include duplicates
893  for(const auto& locPID : locParticles)
894  {
895  auto locPIDIterator = locNumParticles.find(locPID);
896  if(locPIDIterator != locNumParticles.end())
897  ++(locPIDIterator->second);
898  else
899  locNumParticles.emplace(locPID, 1);
900  }
901 
902  return locNumParticles;
903 }
904 
905 DSourceComboUse DSourceComboer::Make_ComboUse(Particle_t locInitPID, const map<Particle_t, unsigned char>& locNumParticles, const map<DSourceComboUse, unsigned char>& locFurtherDecays, bool locMissingDecayProductFlag, Particle_t locTargetToInclude)
906 {
907  //convert locFurtherDecays map to a vector
908  vector<pair<DSourceComboUse, unsigned char>> locDecayVector;
909  locDecayVector.reserve(locFurtherDecays.size());
910  std::copy(locFurtherDecays.begin(), locFurtherDecays.end(), std::back_inserter(locDecayVector));
911 
912  //convert locNumParticles map to a vector
913  vector<pair<Particle_t, unsigned char>> locParticleVector;
914  locParticleVector.reserve(locNumParticles.size());
915  std::copy(locNumParticles.begin(), locNumParticles.end(), std::back_inserter(locParticleVector));
916 
917  //make or get the combo info
918  auto locComboInfo = MakeOrGet_SourceComboInfo(locParticleVector, locDecayVector, 0);
919  auto locComboUse = DSourceComboUse(locInitPID, DSourceComboInfo::Get_VertexZIndex_ZIndependent(), locComboInfo, locMissingDecayProductFlag, locTargetToInclude);
920  if(dDebugLevel >= 5)
921  {
922  cout << "CREATED COMBO USE:" << endl;
923  DAnalysis::Print_SourceComboUse(locComboUse);
924  }
925  return locComboUse;
926 }
927 
928 const DSourceComboInfo* DSourceComboer::MakeOrGet_SourceComboInfo(const vector<pair<Particle_t, unsigned char>>& locNumParticles, const vector<pair<DSourceComboUse, unsigned char>>& locFurtherDecays, unsigned char locNumTabs)
929 {
930  //to be called (indirectly) by constructor: during the stage when primarily making
931  //create the object on the stack
932  DSourceComboInfo locSearchForInfo(locNumParticles, locFurtherDecays);
933 
934  //then search through the set to retrieve the pointer to the corresponding object if it already exists
935  auto locInfoIterator = dSourceComboInfoSet.find(&locSearchForInfo);
936  if(locInfoIterator != dSourceComboInfoSet.end())
937  return *locInfoIterator; //it exists: return it
938 
939  //doesn't exist, make it and insert it into the sorted vector in the correct spot
940  auto locComboInfo = new DSourceComboInfo(locNumParticles, locFurtherDecays);
941  dSourceComboInfoSet.insert(locComboInfo);
942  dComboInfoChargeContent.emplace(locComboInfo, DAnalysis::Get_ChargeContent(locComboInfo));
943  if(dDebugLevel >= 5)
944  {
945  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
946  cout << "CREATED COMBO INFO:" << endl;
947  DAnalysis::Print_SourceComboInfo(locComboInfo, locNumTabs);
948  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
949  cout << "charge content = " << dComboInfoChargeContent[locComboInfo] << endl;
950  }
951  if(DAnalysis::Get_HasMassiveNeutrals(locComboInfo))
952  dComboInfosWithMassiveNeutrals.insert(locComboInfo);
953  if(DAnalysis::Get_HasPhotons(locComboInfo))
954  dComboInfosWithPhotons.insert(locComboInfo);
955  return locComboInfo;
956 }
957 
958 const DSourceComboInfo* DSourceComboer::GetOrMake_SourceComboInfo(const vector<pair<Particle_t, unsigned char>>& locNumParticles, const vector<pair<DSourceComboUse, unsigned char>>& locFurtherDecays, unsigned char locNumTabs)
959 {
960  //to be called when making combos: during the stage when primarily getting
961  //create the object on the stack
962  DSourceComboInfo locSearchForInfo(locNumParticles, locFurtherDecays);
963 
964  //then search through the vector to retrieve the pointer to the corresponding object if it already exists
965  auto locIteratorPair = std::equal_range(dSourceComboInfos.begin(), dSourceComboInfos.end(), &locSearchForInfo, DCompare_SourceComboInfos());
966  if(locIteratorPair.first != locIteratorPair.second)
967  return *(locIteratorPair.first); //it exists: return it
968 
969  //doesn't exist, make it and insert it into the sorted vector in the correct spot
970  auto locComboInfo = new DSourceComboInfo(locNumParticles, locFurtherDecays);
971  dSourceComboInfos.emplace(locIteratorPair.first, locComboInfo);
972  dComboInfoChargeContent.emplace(locComboInfo, DAnalysis::Get_ChargeContent(locComboInfo));
973  if(dDebugLevel >= 5)
974  {
975  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
976  cout << "CREATED COMBO INFO:" << endl;
977  DAnalysis::Print_SourceComboInfo(locComboInfo, locNumTabs);
978  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
979  cout << "charge content = " << dComboInfoChargeContent[locComboInfo] << endl;
980  }
981  if(DAnalysis::Get_HasMassiveNeutrals(locComboInfo))
982  dComboInfosWithMassiveNeutrals.insert(locComboInfo);
983  if(DAnalysis::Get_HasPhotons(locComboInfo))
984  dComboInfosWithPhotons.insert(locComboInfo);
985  return locComboInfo;
986 }
987 
988 DSourceComboUse DSourceComboer::Create_ZDependentSourceComboUses(const DReactionVertexInfo* locReactionVertexInfo, const DSourceCombo* locReactionChargedCombo)
989 {
990  //this creates new uses, with the specific vertex-z bins needed
991  //note that the use can have a different structure from the charged!! (although not likely)
992  //E.g. if something crazy like 2 KShorts -> 3pi, each at a different vertex-z bin, then they will no longer be grouped together vertically (separate uses: horizontally instead)
993 
994  //see if they've already been created. if so, just return it.
995  auto locVertexZBins = dSourceComboVertexer->Get_VertexZBins(locReactionVertexInfo, locReactionChargedCombo, nullptr, true);
996  auto locCreationPair = std::make_pair(locReactionVertexInfo, locVertexZBins);
997  auto locUseIterator = dSourceComboUseVertexZMap.find(locCreationPair);
998  if(locUseIterator != dSourceComboUseVertexZMap.end())
999  return locUseIterator->second; //already created! we are done
1000 
1001  auto locReaction = locReactionVertexInfo->Get_Reaction();
1002 
1003  //loop over vertex infos in reverse-step order
1004  unordered_map<size_t, DSourceComboUse> locCreatedUseMap; //size_t: step index
1005  auto locStepVertexInfos = DAnalysis::Get_StepVertexInfos_ReverseOrderByStep(locReactionVertexInfo);
1006  for(const auto& locStepVertexInfo : locStepVertexInfos)
1007  {
1008  //for this vertex, get the vertex z bin
1009  auto locVertexZBin = (locReactionChargedCombo != nullptr) ? dSourceComboVertexer->Get_VertexZBin(locStepVertexInfo, locReactionChargedCombo, nullptr, true) : dSourceComboTimeHandler->Get_VertexZBin_TargetCenter();
1010  if(dDebugLevel >= 20)
1011  cout << "step, vertex z-bin = " << locStepVertexInfo->Get_StepIndices().front() << ", " << locVertexZBin << endl;
1012 
1013  //loop over the steps at this vertex z bin, in reverse order
1014  auto locStepIndices = locStepVertexInfo->Get_StepIndices();
1015  for(auto locStepIterator = locStepIndices.rbegin(); locStepIterator != locStepIndices.rend(); ++locStepIterator)
1016  {
1017  auto locStepIndex = *locStepIterator;
1018  auto locStepOrigUse = dSourceComboUseReactionStepMap[locReaction][locStepIndex];
1019 
1020  //build new use for the further decays, setting the vertex z-bins
1021  auto locNewComboUse = Build_NewZDependentUse(locReaction, locStepIndex, locVertexZBin, locStepOrigUse, locCreatedUseMap);
1022  locCreatedUseMap.emplace(locStepIndex, locNewComboUse);
1023  }
1024  }
1025 
1026  dSourceComboUseVertexZMap.emplace(locCreationPair, locCreatedUseMap[0]);
1027  return locCreatedUseMap[0];
1028 }
1029 
1030 DSourceComboUse DSourceComboer::Build_NewZDependentUse(const DReaction* locReaction, size_t locStepIndex, signed char locVertexZBin, const DSourceComboUse& locOrigUse, const unordered_map<size_t, DSourceComboUse>& locCreatedUseMap)
1031 {
1032  //each step can be broken up into combo infos with a depth of 2 (grouping charges separately)
1033  auto locStep = locReaction->Get_ReactionStep(locStepIndex);
1034  auto locOrigInfo = std::get<2>(locOrigUse);
1035  if(dComboInfoChargeContent[locOrigInfo] == d_Charged)
1036  {
1037  dZDependentUseToIndependentMap.emplace(locOrigUse, locOrigUse);
1038  return locOrigUse; //no need to change!: no neutrals anyway
1039  }
1040 
1041  map<DSourceComboUse, unsigned char> locNewFurtherDecays;
1042  auto locOrigFurtherDecays = locOrigInfo->Get_FurtherDecays();
1043  for(const auto& locDecayPair : locOrigFurtherDecays)
1044  {
1045  const auto& locOrigDecayUse = locDecayPair.first;
1046  auto locDecayPID = std::get<0>(locOrigDecayUse);
1047  if(locDecayPID != Unknown)
1048  {
1049  //these decays are represented by other steps, and have already been saved
1050  for(unsigned char locInstance = 1; locInstance <= locDecayPair.second; ++locInstance)
1051  {
1052  auto locParticleIndex = DAnalysis::Get_ParticleIndex(locStep, locDecayPID, locInstance);
1053  auto locDecayStepIndex = DAnalysis::Get_DecayStepIndex(locReaction, locStepIndex, locParticleIndex);
1054  const auto& locSavedDecayUse = locCreatedUseMap.find(locDecayStepIndex)->second; //is same as locOrigDecayUse, except different zbins along chain
1055 
1056  //save the use for this decay
1057  auto locUseIterator = locNewFurtherDecays.find(locSavedDecayUse);
1058  if(locUseIterator != locNewFurtherDecays.end())
1059  ++(locUseIterator->second);
1060  else
1061  locNewFurtherDecays.emplace(locSavedDecayUse, 1);
1062  }
1063  }
1064  else //is unknown (and guaranteed to be size 1 since has unknown parent)
1065  {
1066  //must dig down, but only one level: their decays must terminate at new steps (or end)
1067  auto locNewComboUse = Build_NewZDependentUse(locReaction, locStepIndex, locVertexZBin, locOrigDecayUse, locCreatedUseMap);
1068  //save the use for this decay
1069  auto locUseIterator = locNewFurtherDecays.find(locNewComboUse);
1070  if(locUseIterator != locNewFurtherDecays.end())
1071  ++(locUseIterator->second);
1072  else
1073  locNewFurtherDecays.emplace(locNewComboUse, 1);
1074  }
1075  }
1076 
1077  //build and save new info, use, and return
1078  vector<pair<DSourceComboUse, unsigned char>> locFurtherDecayVector;
1079  locFurtherDecayVector.reserve(locNewFurtherDecays.size());
1080  std::copy(locNewFurtherDecays.begin(), locNewFurtherDecays.end(), std::back_inserter(locFurtherDecayVector));
1081  auto locNewComboInfo = locNewFurtherDecays.empty() ? locOrigInfo : GetOrMake_SourceComboInfo(locOrigInfo->Get_NumParticles(), locFurtherDecayVector, 0);
1082 
1083  DSourceComboUse locNewComboUse(std::get<0>(locOrigUse), locVertexZBin, locNewComboInfo, std::get<3>(locOrigUse), std::get<4>(locOrigUse));
1084  if(dDebugLevel >= 30)
1085  {
1086  cout << "NEW Z-DEPENDENT USE:" << endl;
1087  Print_SourceComboUse(locNewComboUse);
1088  cout << "FROM ORIG USE:" << endl;
1089  Print_SourceComboUse(locOrigUse);
1090  }
1091  dZDependentUseToIndependentMap.emplace(locNewComboUse, locOrigUse);
1092  return locNewComboUse;
1093 }
1094 
1095 /********************************************************************** SETUP FOR NEW EVENT ***********************************************************************/
1096 
1097 void DSourceComboer::Reset_NewEvent(JEventLoop* locEventLoop)
1098 {
1099  //check if it's actually a new event
1100  auto locEventNumber = locEventLoop->GetJEvent().GetEventNumber();
1101  if(locEventNumber == dEventNumber)
1102  return; //nope
1103  dEventNumber = locEventNumber;
1104  if(dDebugLevel >= 5) //for the last event
1105  {
1106  cout << "Total # of Combos Allocated (All threads): " << dResourcePool_SourceCombo.Get_NumObjectsAllThreads() << endl;
1107  cout << "Total # of Combo Vectors Allocated (All threads): " << dResourcePool_SourceComboVector.Get_NumObjectsAllThreads() << endl;
1109  }
1110 
1112 
1113  /************************************************************* RECYCLE AND RESET **************************************************************/
1114 
1115  //RECYCLE COMBO & VECTOR POINTERS
1116  //be careful! don't recycle combos with a use pid != unknown, because they are just copies! not unique pointers!
1117 
1118  //HANDLERS AND VERTEXERS
1123 
1124  //PARTICLES
1125  dNumChargedTracks = 0;
1126  dTracksByPID.clear();
1127  dTracksByCharge.clear();
1128  dShowersByBeamBunchByZBin.clear();
1129 
1130  //RECYCLE THE DSOURCECOMBO OBJECTS
1132  decltype(dCreatedCombos)().swap(dCreatedCombos); //should have been reset anyway, but just in case
1133  Recycle_Vectors();
1134 
1135  //COMBOING RESULTS:
1136  dSourceCombosByUse_Charged.clear(); //BEWARE, CONTAINS VECTORS
1137  dMixedCombosByUseByChargedCombo.clear(); //BEWARE, CONTAINS VECTORS
1139  dVertexPrimaryComboMap.clear();
1140  dValidRFBunches_ByCombo.clear();
1141  dNPhotonsToComboMap.clear();
1142 
1143  //COMBOING RESUME/SEARCH-AFTER TRACKING
1146 
1147  /************************************************************ SETUP FOR NEW EVENT *************************************************************/
1148 
1149  //GET JANA OBJECTS
1150  vector<const DNeutralShower*> locNeutralShowers;
1151  locEventLoop->Get(locNeutralShowers, dShowerSelectionTag.c_str());
1152 
1153  vector<const DChargedTrack*> locChargedTracks;
1154  locEventLoop->Get(locChargedTracks, "Combo");
1155 
1156  vector<const DBeamPhoton*> locBeamPhotons;
1157  locEventLoop->Get(locBeamPhotons);
1158 
1159  const DEventRFBunch* locInitialRFBunch = nullptr;
1160  locEventLoop->GetSingle(locInitialRFBunch);
1161 
1162  const DDetectorMatches* locDetectorMatches = nullptr;
1163  locEventLoop->GetSingle(locDetectorMatches, "Combo");
1164 
1165 //COMPARE:
1166 const DVertex* locVertex = nullptr;
1167 locEventLoop->GetSingle(locVertex);
1168 dSourceComboVertexer->Set_Vertex(locVertex);
1169 
1170  vector<const DESSkimData*> locESSkimDataVector;
1171  locEventLoop->Get(locESSkimDataVector);
1172  dESSkimData = locESSkimDataVector.empty() ? NULL : locESSkimDataVector[0];
1173 
1174  //SETUP NEUTRAL SHOWERS
1175  dSourceComboTimeHandler->Setup(locNeutralShowers, locInitialRFBunch, locDetectorMatches);
1178  for(auto& locZBinPair : dShowersByBeamBunchByZBin)
1179  {
1180  auto& locShowerByBunchMap = locZBinPair.second;
1181  if(dDebugLevel >= 20)
1182  cout << "Register zbin: " << int(locZBinPair.first) << endl;
1183  for(auto& locBunchPair : locShowerByBunchMap)
1184  Build_ParticleIndices(Gamma, locBunchPair.first, locBunchPair.second, locZBinPair.first);
1185  }
1186 
1187  //SETUP BEAM PARTICLES
1188  dSourceComboTimeHandler->Set_BeamParticles(locBeamPhotons);
1189 
1190  //SETUP TRACKS
1191  dNumChargedTracks = locChargedTracks.size();
1192  for(const auto& locChargedTrack : locChargedTracks)
1193  {
1194  for(const auto& locChargedHypo : locChargedTrack->dChargedTrackHypotheses)
1195  {
1196  if(dDebugLevel >= 5)
1197  cout << "track, hypo, pid, t1 system = " << locChargedTrack << ", " << locChargedHypo << ", " << locChargedHypo->PID() << ", " << locChargedHypo->t1_detector() << endl;
1198  if(!Cut_dEdxAndEOverP(locChargedHypo))
1199  continue;
1200  if(dDebugLevel >= 5)
1201  cout << "passed cuts, register" << endl;
1202  dTracksByPID[locChargedHypo->PID()].push_back(locChargedTrack);
1203  dTracksByCharge[ParticleCharge(locChargedHypo->PID()) > 0].push_back(locChargedTrack); //will insert duplicates
1204  }
1205  }
1206 
1207  //sort by pid & create indices
1208  for(auto& locPIDPair : dTracksByPID)
1209  {
1210  std::sort(locPIDPair.second.begin(), locPIDPair.second.end());
1211  Build_ParticleIndices(locPIDPair.first, {}, locPIDPair.second, DSourceComboInfo::Get_VertexZIndex_ZIndependent());
1212  }
1213  //sort & remove duplicates in tracks-by-charge
1214  for(auto& locChargePair : dTracksByCharge)
1215  {
1216  auto& locVector = locChargePair.second;
1217  std::sort(locVector.begin(), locVector.end());
1218  locVector.erase(std::unique(locVector.begin(), locVector.end()), locVector.end()); //remove duplicates
1219  }
1220 
1221  if(dDebugLevel > 0)
1222  {
1223  cout << "TRACKS BY PID:" << endl;
1224  for(const auto& locPIDPair : dTracksByPID)
1225  {
1226  cout << "PID, pointers: " << locPIDPair.first << ", ";
1227  for(const auto& locTrack : locPIDPair.second)
1228  cout << locTrack << ", ";
1229  cout << endl;
1230  }
1231  cout << "TRACKS BY CHARGE:" << endl;
1232  for(const auto& locChargePair : dTracksByCharge)
1233  {
1234  cout << "charge, pointers: " << (locChargePair.first ? 1 : -1) << ", ";
1235  for(const auto& locTrack : locChargePair.second)
1236  cout << locTrack << ", ";
1237  cout << endl;
1238  }
1239  }
1240 
1241  //Fill histograms
1242  Fill_CutHistograms();
1243 }
1244 
1245 bool DSourceComboer::Cut_dEdxAndEOverP(const DChargedTrackHypothesis* locChargedTrackHypothesis)
1246 {
1247  auto locPID = locChargedTrackHypothesis->PID();
1248  auto locTrackTimeBased = locChargedTrackHypothesis->Get_TrackTimeBased();
1249  auto locP = locTrackTimeBased->momentum().Mag();
1250  bool locPassedCutFlag = true;
1251 
1252  //CDC dE/dx
1253 //cout << "PID, p, dedx, #hits = " << locPID << ", " << locP << ", " << locTrackTimeBased->ddEdx_CDC*1.0E6 << ", " << locTrackTimeBased->dNumHitsUsedFordEdx_CDC << endl;
1254  if(locTrackTimeBased->dNumHitsUsedFordEdx_CDC > 0)
1255  {
1256  auto locdEdx = locTrackTimeBased->ddEdx_CDC_amp*1.0E6;
1257  if(!Cut_dEdx(locPID, SYS_CDC, locP, locdEdx))
1258  locPassedCutFlag = false;
1259  }
1260  else if((locPID == KPlus) || (locPID == KMinus))
1261 // if((locPID == KPlus) || (locPID == KMinus)) //COMPARE: use this instead
1262  {
1263  auto locSystem = locChargedTrackHypothesis->t1_detector();
1264  if((locSystem == SYS_START) || (locSystem == SYS_NULL))
1265  return false; //kaons are rare, and no PID information to find them (swamped with background): just cut these away
1266  }
1267 
1268  //FDC dE/dx
1269  if(locTrackTimeBased->dNumHitsUsedFordEdx_FDC > 0)
1270  {
1271  auto locdEdx = locTrackTimeBased->ddEdx_FDC*1.0E6;
1272  if(!Cut_dEdx(locPID, SYS_FDC, locP, locdEdx))
1273  locPassedCutFlag = false;
1274  }
1275 
1276  //SC dE/dx
1277  auto locSCHitMatchParams = locChargedTrackHypothesis->Get_SCHitMatchParams();
1278  if(locSCHitMatchParams != nullptr)
1279  {
1280  auto locdEdx = locSCHitMatchParams->dEdx*1.0E3;
1281  if(!Cut_dEdx(locPID, SYS_START, locP, locdEdx))
1282  locPassedCutFlag = false;
1283  }
1284 
1285  //TOF dE/dx
1286  auto locTOFHitMatchParams = locChargedTrackHypothesis->Get_TOFHitMatchParams();
1287  if(locTOFHitMatchParams != nullptr)
1288  {
1289  auto locdEdx = locTOFHitMatchParams->dEdx*1.0E3;
1290  if(!Cut_dEdx(locPID, SYS_TOF, locP, locdEdx))
1291  locPassedCutFlag = false;
1292  }
1293 
1294  //BCAL E/p
1295  auto locBCALShowerMatchParams = locChargedTrackHypothesis->Get_BCALShowerMatchParams();
1296  if(locBCALShowerMatchParams != nullptr)
1297  {
1298  const DBCALShower* locBCALShower = locBCALShowerMatchParams->dBCALShower;
1299  double locEOverP = locBCALShower->E/locP;
1300  if(!Cut_EOverP(locPID, SYS_BCAL, locP, locEOverP))
1301  locPassedCutFlag = false;
1302  }
1303 
1304  //FCAL E/p
1305  auto locFCALShowerMatchParams = locChargedTrackHypothesis->Get_FCALShowerMatchParams();
1306  if(locFCALShowerMatchParams != nullptr)
1307  {
1308  const DFCALShower* locFCALShower = locFCALShowerMatchParams->dFCALShower;
1309  double locEOverP = locFCALShower->getEnergy()/locP;
1310  if(!Cut_EOverP(locPID, SYS_FCAL, locP, locEOverP))
1311  locPassedCutFlag = false;
1312  }
1313 
1314  return locPassedCutFlag;
1315 }
1316 
1318 {
1319  japp->WriteLock("DSourceComboer_Cuts");
1320  {
1321  for(auto& locPIDPair : dHistMap_dEdx)
1322  {
1323  for(auto& locSystemPair : locPIDPair.second)
1324  {
1325  auto& locHist = locSystemPair.second;
1326  auto& locVector = ddEdxValueMap[locPIDPair.first][locSystemPair.first];
1327  for(auto& locVectorPair : locVector)
1328  locHist->Fill(locVectorPair.first, locVectorPair.second);
1329  }
1330  }
1331  for(auto& locPIDPair : dHistMap_EOverP)
1332  {
1333  for(auto& locSystemPair : locPIDPair.second)
1334  {
1335  auto& locHist = locSystemPair.second;
1336  auto& locVector = dEOverPValueMap[locPIDPair.first][locSystemPair.first];
1337  for(auto& locVectorPair : locVector)
1338  locHist->Fill(locVectorPair.first, locVectorPair.second);
1339  }
1340  }
1341  }
1342  japp->Unlock("DSourceComboer_Cuts");
1343 
1344  //Reset for next event
1345  for(auto& locPIDPair : ddEdxValueMap)
1346  {
1347  for(auto& locSystemPair : locPIDPair.second)
1348  decltype(locSystemPair.second)().swap(locSystemPair.second);
1349  }
1350  for(auto& locPIDPair : dEOverPValueMap)
1351  {
1352  for(auto& locSystemPair : locPIDPair.second)
1353  decltype(locSystemPair.second)().swap(locSystemPair.second);
1354  }
1355 }
1356 
1357 /********************************************************************* CREATE DSOURCOMBO'S **********************************************************************/
1358 
1360 {
1361  //This builds the combos and creates DParticleCombo & DParticleComboSteps (doing whatever is necessary)
1362  if(dDebugLevel > 0)
1363  cout << "CREATING DSourceCombo's FOR DREACTION " << locReactionVertexInfo->Get_Reaction()->Get_ReactionName() << endl;
1364 
1365  //Initialize results to be returned
1366  DCombosByReaction locOutputComboMap;
1367  auto locReactions = locReactionVertexInfo->Get_Reactions();
1368  for(auto locReaction : locReactions)
1369  {
1370  locOutputComboMap[locReaction] = {};
1371  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::Input] = 1; //is really #events
1372  }
1373 
1374  if(!Check_Reactions(locReactions))
1375  {
1376  if(dDebugLevel > 0)
1377  {
1378  cout << "FINISHED COMBOING:" << endl;
1379  for(auto locComboPair : locOutputComboMap)
1380  cout << "event#, reaction, #combos = " << dEventNumber << ", " << locComboPair.first->Get_ReactionName() << ", " << locComboPair.second.size() << endl;
1381  }
1382  return locOutputComboMap; //no combos!
1383  }
1384 
1385  /******************************************************** COMBOING STEPS *******************************************************
1386  *
1387  * CHARGED STAGE:
1388  *
1389  * OK, we start with charged tracks, because we can't do much with neutrals until we know the vertex to compute the kinematics.
1390  * So, we create our combos, skipping all neutral particles, but filling in all charged tracks.
1391  *
1392  * If mass cuts are needed (e.g. Lambda -> p, pi-), we first create combos of "-> p, pi-", saving them for the USE "X -> p, pi-"
1393  * We then place the invariant mass cut, and those that pass get copied and saved for the USE "Lambda -> p, pi-"
1394  * Thus, storing the decay PID separately from the combo means we can reuse the combo without creating new objects in this case.
1395  *
1396  * Once we have our charged track combos, we can find (most of) the vertices (will discuss exceptions below).
1397  * Once we have the vertices, we can compute the time offsets between the vertices (the amount of time a decaying particle took to decay).
1398  * And we can then place timing cuts on the charged tracks to select which beam bunches are possible.
1399  * Now, you might be thinking that we can cut on the timing of the charged tracks BEFORE we find the vertices, but in some cases we can't.
1400  * For a discussion on this, see the comments in DSourceComboTimeHandler.
1401  *
1402  *
1403  *
1404  * MIXED STAGE: GENERAL
1405  *
1406  * OK, now let's combo some neutrals.
1407  * First, we combo all of the neutrals that are needed with each other, and THEN we combo them with charged tracks.
1408  * (This is how the DSourceComboInfo objects were constructed).
1409  * This is because pi0 comboing will take the longest, and we want to make sure it is done largely independent of any charged tracks.
1410  *
1411  *
1412  * MIXED STAGE: VERTEX-Z
1413  * Now, as discussed earlier, showers can be broken up into z-dependent and z-independent varieties.
1414  * Z-Independent: FCAL photons
1415  * Z-Dependent: BCAL showers or FCAL massive neutrals
1416  * However, since
1417  * Again, for details, see the comments in DSourceComboTimeHandler and DSourceComboP4Handler.
1418  *
1419  * Now, since the z-independent combos can be reused for any vertex-z, they are created first.
1420  * Then, the z-dependent combos are created, and combined with the z-independent ones.
1421  * To do this, it turns out it's easier to just try to create combos with ALL showers, and then skip creating the ones we've already created.
1422  *
1423  * While building combos, mass cuts are placed along the way, EXCEPT on combos with massive neutral particles.
1424  * This is because the exact vertex position is needed to get an accurate massive-neutral momentum.
1425  * While comboing, we want the results to be as re-usable as possible, that's why we use vertex-z bins.
1426  * But vertex-z bins are not sufficient for this, so we will cut on invariant masses with massive neutrals later.
1427  *
1428  *******************************************************************************************************************************/
1429 
1430  //MUST BEWARE DUPLICATE COMBOS
1431  //let's say a combo of charged tracks has 2 valid RF bunches
1432  //and we need to combo 2 pi0s with them
1433  //and the shower timing cuts are loose enough that all 4 showers satisfy both RF bunches
1434  //if we combo the 2 rf bunches separately: WE HAVE DUPLICATE COMBOS
1435  //and doing the duplicate check AFTER the fact takes FOREVER
1436  //therefore, we must take the neutral showers for the 2 rfs, COMBINE THEM, and then COMBO AS A UNIT
1437 
1438  //get step vertex infos (sorted in dependency order)
1439  auto locStepVertexInfos = locReactionVertexInfo->Get_StepVertexInfos();
1440  auto locPrimaryStepVertexInfo = locReactionVertexInfo->Get_StepVertexInfo(0);
1441  auto locPrimaryComboUse = dSourceComboUseReactionMap[locPrimaryStepVertexInfo];
1442  auto locPrimaryComboInfo = std::get<2>(locPrimaryComboUse);
1443 
1444  //handle special case of no charged tracks
1445  if(dDebugLevel > 0)
1446  cout << "Combo charge content: " << dComboInfoChargeContent[std::get<2>(locPrimaryComboUse)] << " (charged/neutral are " << d_Charged << "/" << d_Neutral << ")" << endl;
1447  if(dComboInfoChargeContent[std::get<2>(locPrimaryComboUse)] == d_Neutral)
1448  {
1449  if(dDebugLevel > 0)
1450  cout << "No charged tracks." << endl;
1451  for(auto& locReaction : locReactions)
1452  {
1453  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::Charged_Combos] = 1; //is really #-events
1454  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::Charged_RFBunch] = 1; //is really #-events
1455  }
1456  Combo_WithNeutralsAndBeam(locReactions, locReactionVertexInfo, locPrimaryComboUse, nullptr, {}, locOutputComboMap);
1457 
1458  if(dDebugLevel > 0)
1459  {
1460  cout << "FINISHED COMBOING:" << endl;
1461  for(auto locComboPair : locOutputComboMap)
1462  cout << "event#, reaction, #combos = " << dEventNumber << ", " << locComboPair.first->Get_ReactionName() << ", " << locComboPair.second.size() << endl;
1463  }
1464  return locOutputComboMap;
1465  }
1466 
1467  //Build vertex combos (returns those for the primary vertex, others are stored)
1468  Create_SourceCombos(locPrimaryComboUse, d_ChargedStage, nullptr, 0);
1469  const auto& locReactionChargedCombos = *(Get_CombosSoFar(d_ChargedStage, d_Charged, nullptr)[locPrimaryComboUse]);
1470  for(auto& locReaction : locReactions)
1471  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::Charged_Combos] = locReactionChargedCombos.size();
1472 
1473  if(dDebugLevel > 0)
1474  {
1475  cout << "Charged combos built: " << locReactionChargedCombos.size() << endl;
1476  if(locReactionChargedCombos.empty())
1477  cout << "no combos for event: " << dEventNumber << endl;
1478  }
1479 
1480  //loop over primary vertex combos //each contains decay combos except when dangling
1481  for(const auto& locReactionChargedCombo : locReactionChargedCombos)
1482  {
1483  //Calc all the vertex positions and time offsets for the vertices for these combos (where possible without beam energy)
1484  dSourceComboVertexer->Calc_VertexTimeOffsets_WithCharged(locReactionVertexInfo, locReactionChargedCombo);
1485 
1486  //For the charged tracks, apply timing cuts to determine which RF bunches are possible
1487  vector<int> locBeamBunches_Charged;
1488  if(!dSourceComboTimeHandler->Select_RFBunches_Charged(locReactionVertexInfo, locReactionChargedCombo, locBeamBunches_Charged))
1489  continue; //failed PID timing cuts!
1490  for(auto& locReaction : locReactions)
1492 
1493  //Special case of FULLY charged
1494  auto locChargeContent = dComboInfoChargeContent[locPrimaryComboInfo];
1495  if(locChargeContent == d_Charged)
1496  {
1497  if(dDebugLevel > 0)
1498  cout << "Fully charged." << endl;
1499 
1500  if(false) //COMPARE: Comparison-to-old mode
1501  {
1502  dSourceComboTimeHandler->Vote_OldMethod(locReactionChargedCombo, locBeamBunches_Charged);
1503  if(locBeamBunches_Charged.empty())
1504  continue;
1505  }
1506 
1507  //Select final RF bunch
1508  auto locRFBunch = dSourceComboTimeHandler->Select_RFBunch_Full(locReactionVertexInfo, locReactionChargedCombo, locBeamBunches_Charged);
1509  if(dDebugLevel > 0)
1510  cout << "Selected rf bunch." << endl;
1511 
1512  for(auto& locReaction : locReactions)
1513  {
1518  }
1519 
1520  //combo with beam and save results!!! (if no beam needed, just saves and returns)
1521  Combo_WithBeam(locReactions, locReactionVertexInfo, locPrimaryComboUse, locReactionChargedCombo, locRFBunch, locOutputComboMap);
1522  continue;
1523  }
1524 
1525  //Combo with neutrals and beam
1526  Combo_WithNeutralsAndBeam(locReactions, locReactionVertexInfo, locPrimaryComboUse, locReactionChargedCombo, locBeamBunches_Charged, locOutputComboMap);
1527  }
1528 
1529  if(dDebugLevel > 0)
1530  {
1531  cout << "FINISHED COMBOING:" << endl;
1532  for(auto locComboPair : locOutputComboMap)
1533  cout << "event#, reaction, #combos = " << dEventNumber << ", " << locComboPair.first->Get_ReactionName() << ", " << locComboPair.second.size() << endl;
1534  }
1535 
1536  return locOutputComboMap;
1537 }
1538 
1539 void DSourceComboer::Combo_WithNeutralsAndBeam(const vector<const DReaction*>& locReactions, const DReactionVertexInfo* locReactionVertexInfo, const DSourceComboUse& locPrimaryComboUse, const DSourceCombo* locReactionChargedCombo, const vector<int>& locBeamBunches_Charged, DCombosByReaction& locOutputComboMap)
1540 {
1541  if(dDebugLevel > 0)
1542  {
1543  auto locNumDetectedShowers = dShowersByBeamBunchByZBin[DSourceComboInfo::Get_VertexZIndex_Unknown()][{}].size();
1544  auto locNumFCALShowers = dShowersByBeamBunchByZBin[DSourceComboInfo::Get_VertexZIndex_ZIndependent()][{}].size();
1545  cout << endl << "Comboing neutrals, z-independent, #FCAL/BCAL showers: " << locNumFCALShowers << "/" << locNumDetectedShowers - locNumFCALShowers << endl;
1546  }
1547 
1548  if(dDebugLevel >= 5)
1549  {
1550  cout << "charged combo, vertex zbins: " << locReactionChargedCombo;
1551  auto locVertexZBins = dSourceComboVertexer->Get_VertexZBins(locReactionVertexInfo, locReactionChargedCombo, nullptr, true);
1552  for(auto& locZBin : locVertexZBins)
1553  cout << ", " << int(locZBin);
1554  cout << endl;
1555  }
1556  //if there is a vertex zbin that is out of range, and we need photons: don't allow: will blow up memory due to no invariant mass cuts
1557  for(auto& locStepVertexInfo : locReactionVertexInfo->Get_StepVertexInfos())
1558  {
1559  auto locZBin = dSourceComboVertexer->Get_VertexZBin(locStepVertexInfo, locReactionChargedCombo, nullptr, true);
1561  continue;
1562  if(!locStepVertexInfo->Get_OnlyConstrainTimeParticles().empty())
1563  {
1564  if(dDebugLevel > 0)
1565  cout << "Combo has photons at a vertex that is out of range: don't combo." << endl;
1566  return; //has photons, don't combo
1567  }
1568  }
1569 
1570  //Create full source-particle combos (including neutrals): First using only FCAL showers, then using all showers
1571  Create_SourceCombos(locPrimaryComboUse, d_MixedStage_ZIndependent, locReactionChargedCombo, 0);
1572  auto locZDependentComboUse = Create_ZDependentSourceComboUses(locReactionVertexInfo, locReactionChargedCombo);
1573  if(dDebugLevel > 0)
1574  cout << endl << "Comboing neutrals, z-dependent." << endl;
1575  Create_SourceCombos(locZDependentComboUse, d_MixedStage, locReactionChargedCombo, 0);
1576 
1577  //Then, get the full combos, but only those that satisfy the charged RF bunches
1578  vector<int> locComboRFBunches = locBeamBunches_Charged;
1579  //However, if there are no photons (only massive neutrals), then all of the combos have been saved with RF bunches = {} (empty set)
1580  if(!Get_HasPhotons(std::get<2>(locPrimaryComboUse)))
1581  locComboRFBunches.clear();
1582  const auto& locReactionFullCombos = Get_CombosForComboing(locZDependentComboUse, d_MixedStage, locComboRFBunches, locReactionChargedCombo);
1583  if(dDebugLevel > 0)
1584  cout << endl << "Neutral combos created, # with the charged RF bunches: " << locReactionFullCombos.size() << endl;
1585  for(auto& locReaction : locReactions)
1586  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::Full_Combos] += locReactionFullCombos.size();
1587 
1588  if((dDebugLevel > 0) || (dDebugLevel == -1))
1589  Check_ForDuplicates(locReactionFullCombos);
1590 
1591  //loop over full combos
1592  for(const auto& locReactionFullCombo : locReactionFullCombos)
1593  {
1594  //get common RF bunches between charged & neutral
1595  auto locNeutralRFBunches = dValidRFBunches_ByCombo[std::make_pair(locReactionFullCombo, DSourceComboInfo::Get_VertexZIndex_ZIndependent())];
1596  auto locValidRFBunches = dSourceComboTimeHandler->Get_CommonRFBunches(locBeamBunches_Charged, locNeutralRFBunches);
1597  if(dDebugLevel > 0)
1598  cout << "#charged bunches, #neutral, #common = " << locBeamBunches_Charged.size() << ", " << locNeutralRFBunches.size() << ", " << locValidRFBunches.size() << endl;
1599  if(locValidRFBunches.empty() && (!locNeutralRFBunches.empty() || !locBeamBunches_Charged.empty()))
1600  continue; //fail RF bunch cut
1601 
1602  //if not fully neutral (at least one vertex position is known), do the below
1603  if(locReactionChargedCombo != nullptr)
1604  {
1605  //Calculate vertex positions & time offsets using photons
1606  //not likely to have any effect, but it's necessary sometimes (but rarely)
1607  //E.g. g, p -> K0, Sigma+ K0 -> 3pi: The selected pi0 photons could help define the production vertex
1608  dSourceComboVertexer->Calc_VertexTimeOffsets_WithPhotons(locReactionVertexInfo, locReactionChargedCombo, locReactionFullCombo);
1609 
1610  //Now further select rf bunches, using tracks at the vertices we just found, and BCAL photon showers at any vertex
1611  //this also does PID cuts of photons at charged vertices while we're at it
1612  if(!dSourceComboTimeHandler->Select_RFBunches_PhotonVertices(locReactionVertexInfo, locReactionFullCombo, locValidRFBunches))
1613  {
1614  if(dDebugLevel > 0)
1615  cout << "Failed photon/photon-vertex PID timing cuts" << endl;
1616  continue; //failed PID timing cuts!
1617  }
1618  for(auto& locReaction : locReactions)
1620 
1621  //if no valid RF bunches, but still haven't cut: none of the charged tracks are at known vertices: select RF bunches with charged only
1622  if(locValidRFBunches.empty()) //e.g. g, p -> K0, Sigma+ K0 -> pi+, (pi-)
1623  {
1624  if(!dSourceComboTimeHandler->Select_RFBunches_AllVerticesUnknown(locReactionVertexInfo, locReactionFullCombo, d_Charged, locValidRFBunches))
1625  continue; //failed PID timing cuts!
1626  }
1627  for(auto& locReaction : locReactions)
1629  }
1630  else //fully neutral, so no known vertices, target center was chosen as vertex for comboing showers //e.g. g, p -> pi0, (p)
1631  {
1632  //we will never have a vertex, so do PID cuts for ALL photons using target center to select possible RF bunches
1633  if(!dSourceComboTimeHandler->Select_RFBunches_AllVerticesUnknown(locReactionVertexInfo, locReactionFullCombo, d_Neutral, locValidRFBunches))
1634  continue; //failed PID timing cuts!
1635  for(auto& locReaction : locReactions)
1637  }
1638 
1639  if(false) //COMPARE: Comparison-to-old mode
1640  {
1641  dSourceComboTimeHandler->Vote_OldMethod(locReactionFullCombo, locValidRFBunches);
1642  if(locValidRFBunches.empty())
1643  continue;
1644  }
1645 
1646  //Place mass cuts on massive neutrals: Effectively narrows down RF bunches
1647  //do 2 things at once (where vertex is known) (hence the really long function name):
1648  //calc & cut invariant mass: when massive neutral present
1649  //calc & cut invariant mass: when vertex-z was unknown with only charged tracks, but is known now, and contains BCAL photons (won't happen very often)
1650  if(!dSourceComboP4Handler->Cut_InvariantMass_HasMassiveNeutral_OrPhotonVertex(locReactionVertexInfo, locReactionFullCombo, locValidRFBunches))
1651  continue; //failed cut!
1652  for(auto& locReaction : locReactions)
1654 
1655  //Select final RF bunch //this is not a cut: at least one has passed all cuts (check by the Get_CombosForComboing function & the mass cuts)
1656  auto locRFBunch = dSourceComboTimeHandler->Select_RFBunch_Full(locReactionVertexInfo, locReactionFullCombo, locValidRFBunches);
1657 
1658  //combo with beam and save results!!! (if no beam needed, just saves and returns)
1659  Combo_WithBeam(locReactions, locReactionVertexInfo, locZDependentComboUse, locReactionFullCombo, locRFBunch, locOutputComboMap);
1660  }
1661 }
1662 
1663 void DSourceComboer::Combo_WithBeam(const vector<const DReaction*>& locReactions, const DReactionVertexInfo* locReactionVertexInfo, const DSourceComboUse& locReactionFullComboUse, const DSourceCombo* locReactionFullCombo, int locRFBunch, DCombosByReaction& locOutputComboMap)
1664 {
1665  if(dDebugLevel > 0)
1666  cout << endl << "Comboing beam." << endl;
1667 
1668  //if no beam then we are done!
1669  if(!locReactionVertexInfo->Get_StepVertexInfo(0)->Get_ProductionVertexFlag())
1670  {
1671  if(dDebugLevel > 0)
1672  cout << "No beam particles, we are done!" << endl;
1673 
1674  //place invariant mass cuts using accurate photon kinematics
1675  auto locPassMassCutFlag = dSourceComboP4Handler->Cut_InvariantMass_AccuratePhotonKinematics(locReactionVertexInfo, locReactionFullCombo, nullptr, locRFBunch);
1676  for(const auto& locReaction : locReactions)
1677  {
1681  if(!locPassMassCutFlag)
1682  continue; //FAILED MASS CUTS!
1683 
1687  locOutputComboMap[locReaction].push_back(dParticleComboCreator->Build_ParticleCombo(locReactionVertexInfo, locReactionFullCombo, nullptr, locRFBunch, locReaction->Get_KinFitType()));
1688  }
1689  return;
1690  }
1691 
1692  //Select beam particles
1693  if (abs(locRFBunch) > 2000000000)
1694  return; // proximity to INT_MAX can cause infinite loops, certainly no valid beam particle
1695 
1696  auto locBeamParticles = dSourceComboTimeHandler->Get_BeamParticlesByRFBunch(locRFBunch, dMaxRFBunchCuts[locReactionVertexInfo]);
1697  if(dDebugLevel > 0)
1698  cout << "rf bunch, max #rf bunches, #beams = " << locRFBunch << ", " << dMaxRFBunchCuts[locReactionVertexInfo] << ", " << locBeamParticles.size() << endl;
1699  if(locBeamParticles.empty())
1700  return; //no valid beam particles!!
1701  for(const auto& locReaction : locReactions)
1702  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::Beam_Combos] += locBeamParticles.size();
1703 
1704  //loop over beam particles
1705  for(const auto& locBeamParticle : locBeamParticles)
1706  {
1707  //Calculate remaining vertex positions (that needed to be done via missing mass)
1708  dSourceComboVertexer->Calc_VertexTimeOffsets_WithBeam(locReactionVertexInfo, locReactionFullComboUse, locReactionFullCombo, locBeamParticle);
1709 
1710  //placing timing cuts on the particles at these vertices
1711  if(!dSourceComboTimeHandler->Cut_Timing_MissingMassVertices(locReactionVertexInfo, locReactionFullCombo, locBeamParticle, locRFBunch))
1712  continue; //FAILED TIME CUTS!
1713  for(const auto& locReaction : locReactions)
1715 
1716  //place invariant mass cuts on the particles at these vertices (if they had z-dependent neutral showers (BCAL or massive))
1717  if(!dSourceComboP4Handler->Cut_InvariantMass_MissingMassVertex(locReactionVertexInfo, locReactionFullCombo, locBeamParticle, locRFBunch))
1718  continue; //FAILED MASS CUTS!
1719  for(const auto& locReaction : locReactions)
1721 
1722  //place invariant mass cuts using accurate photon kinematics
1723  if(!dSourceComboP4Handler->Cut_InvariantMass_AccuratePhotonKinematics(locReactionVertexInfo, locReactionFullCombo, locBeamParticle, locRFBunch))
1724  continue; //FAILED MASS CUTS!
1725  for(const auto& locReaction : locReactions)
1727 
1728  //loop over reactions: cut on rf-bunch shift for each reaction, cut on missing mass^2, then save the results
1729  auto locBeamRFBunch = dSourceComboTimeHandler->Calc_RFBunchShift(locBeamParticle->time());
1730  size_t locDeltaRFBunch = abs(locRFBunch - locBeamRFBunch);
1731  for(const auto& locReaction : locReactions)
1732  {
1733  if(dDebugLevel > 0)
1734  cout<< "beam rf bunch, delta rf bunch, reaction, max for reaction = " << locBeamRFBunch << ", " << locDeltaRFBunch << ", " << locReaction->Get_ReactionName() << ", " << dRFBunchCutsByReaction[locReaction] << endl;
1735  if(locDeltaRFBunch > dRFBunchCutsByReaction[locReaction])
1736  continue; //FAILED RF BUNCH CUT
1738 
1739  if(!dSourceComboP4Handler->Cut_MissingMassSquared(locReaction, locReactionVertexInfo, locReactionFullComboUse, locReactionFullCombo, locBeamParticle, locRFBunch))
1740  continue; //FAILED MISSING MASS^2 CUT!
1742 
1743  //build particle combo & save
1744  locOutputComboMap[locReaction].push_back(dParticleComboCreator->Build_ParticleCombo(locReactionVertexInfo, locReactionFullCombo, locBeamParticle, locRFBunch, locReaction->Get_KinFitType()));
1745  if(dDebugLevel >= 10)
1746  {
1747  cout << "Created particle combo, beam energy, combo contents = " << locBeamParticle->energy() << endl;
1748  DAnalysis::Print_SourceCombo(locReactionFullCombo);
1749  }
1750  }
1751  }
1752 }
1753 
1754 /**************************************************************** BUILD SOURCE COMBOS - GENERAL *****************************************************************/
1755 
1756 /*
1757  * suppose reaction is 0) g, p -> omega, p
1758  * 1) omega -> 3pi
1759  * 2) pi0 -> 2g
1760  *
1761  * It will have uses/infos like:
1762  * 0: X -> 1, A (mixed + charged) (both are listed as further decays)
1763  * A: X -> p (charged)
1764  * 1: omega -> B, 2 (mixed) (both are listed as further decays)
1765  * B: X -> pi+, pi- (charged)
1766  * 2: pi0 -> 2g (neutral)
1767  *
1768  * The purpose of passing through the charged combo:
1769  * 1) To retrieve the correct charged combo when comboing it to neutrals to create mixed
1770  * 2) To save the mixed comboing results in a way that they can be reused
1771  *
1772  * The charged combos will be:
1773  * 0: X -> A, 1 //presiding = 0, withnow = A
1774  * A: X -> p //both = nullptr
1775  * 1: omega -> B, 2 //presiding = 1, withnow = B
1776  * B: X -> pi+, pi- //both = nullptr
1777  * 2: pi0 -> 2g //both = nullptr
1778  *
1779  * suppose reaction is 0) g, p -> K0, Sigma+
1780  * 1) K0 -> 3pi
1781  * 2) pi0 -> 2g
1782  * 3) Sigma+ -> pi+, n
1783  *
1784  * It will have uses/infos like:
1785  * 0: X -> A, 1, 3 (mixed -> charged, mixed, mixed) //presiding = 0, withnow = A
1786  * A: X -> p (charged) //both = nullptr
1787  * 1: K0 -> B, 2 (mixed -> charged, neutral) //presiding = 1, withnow = B
1788  * B: X -> pi+, pi- (charged) //both = nullptr
1789  * 2: pi0 -> 2g (neutral) //both = nullptr
1790  * 3: Sigma+ -> C, n (mixed -> charged, n) //presiding = 3, withnow = C
1791  * C: X -> pi+ (charged) //both = nullptr
1792  */
1793 
1794 void DSourceComboer::Create_SourceCombos(const DSourceComboUse& locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, unsigned char locNumTabs)
1795 {
1796  if(dDebugLevel > 0)
1797  {
1798  cout << endl;
1799  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1800  cout << "Creating source combos: Stage, presiding charged combo: " << locComboingStage << ", " << locChargedCombo_Presiding << endl;
1801  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1802  cout << "PRESIDING COMBO:" << endl;
1803  DAnalysis::Print_SourceCombo(locChargedCombo_Presiding, locNumTabs);
1804  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1805  cout << "USE TO CREATE:" << endl;
1806  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
1807  }
1808 
1809  //if on mixed stage, it is impossible for this function to be called with a fully-charged use (already exists!!)
1810  const auto& locDecayPID = std::get<0>(locComboUseToCreate);
1811  const auto& locVertexZBin = std::get<1>(locComboUseToCreate);
1812  const auto& locSourceComboInfo = std::get<2>(locComboUseToCreate);
1813  const auto& locMissingDecayProductFlag = std::get<3>(locComboUseToCreate);
1814 
1815  //Get combos so far
1816  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, dComboInfoChargeContent[locSourceComboInfo], locChargedCombo_Presiding);
1817  if(locSourceCombosByUseSoFar.find(locComboUseToCreate) != locSourceCombosByUseSoFar.end())
1818  {
1819  if(dDebugLevel > 0)
1820  cout << "Already created!" << endl;
1821  return; //we're done!
1822  }
1823 
1824  //we will create these combos for an "Unknown" decay (i.e. no decay, just a direct grouping) (unless already created!)
1825  //then, when we return from this function, we can cut on the invariant mass of the system for any decay we might need it for
1826  DSourceComboUse locUnknownComboUse(Unknown, locVertexZBin, locSourceComboInfo, false, Unknown);
1827  if(locSourceCombosByUseSoFar.find(locUnknownComboUse) == locSourceCombosByUseSoFar.end())
1828  Create_SourceCombos_Unknown(locUnknownComboUse, locComboingStage, locChargedCombo_Presiding, locNumTabs);
1829 
1830  if(dDebugLevel > 0)
1831  {
1832  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1833  cout << "Combos with unknown parent created, desired decay pid = " << locDecayPID << endl;
1834  }
1835 
1836  //if all we want is a direct grouping (unknown), then the combos have already been made: return
1837  if(locDecayPID == Unknown)
1838  return;
1839 
1840  //get the combos that we just created
1841  auto locInfoChargeContent = dComboInfoChargeContent[locSourceComboInfo];
1842  auto locSourceCombos = locSourceCombosByUseSoFar[locUnknownComboUse];
1843 
1844  if((locComboingStage == d_ChargedStage) && (locInfoChargeContent != d_Charged))
1845  {
1846  //don't cut yet! we don't have the neutrals! just copy results and return
1847  if(dDebugLevel > 0)
1848  {
1849  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1850  cout << "On charged stage, need neutrals: done for now" << endl;
1851  }
1852  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, locSourceCombos);
1853 
1854  //Set the resume indices
1855  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
1856  return;
1857  }
1858 
1859  if(locMissingDecayProductFlag)
1860  {
1861  //Don't cut! just copy results and return
1862  if(dDebugLevel > 0)
1863  {
1864  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1865  cout << "Missing decay product: No invariant mass cut." << endl;
1866  }
1867  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, locSourceCombos);
1868 
1869  //Set the resume indices
1870  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
1871  return;
1872  }
1873 
1874  //get combos by beam bunch
1875  auto* locSourceCombosByBeamBunchByUse = (locComboingStage != d_ChargedStage) ? &(Get_SourceCombosByBeamBunchByUse(locInfoChargeContent, locChargedCombo_Presiding)) : nullptr;
1876 
1877  //cannot place an invariant mass cut on massive neutrals yet, because:
1878  //vertex position must first be defined
1879  //although we probably HAVE the vertex position, if it's a fully neutral combo, we don't want to use it:
1880  //results are stored in vertex-z-bins and independent of charged combo: if we cut, we won't be able to reuse the results (because we need PRECISE position, not just a z-bin)
1881  //if it is a mixed combo with known vertex, we can conceivably cut, but there aren't too many of those: Just put off the cuts until later
1882  if(Get_HasMassiveNeutrals(locSourceComboInfo))
1883  {
1884  if(dDebugLevel > 0)
1885  {
1886  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1887  cout << "Massive neutrals, done for now" << endl;
1888  }
1889  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, locSourceCombos);
1890  (*locSourceCombosByBeamBunchByUse)[locComboUseToCreate] = (*locSourceCombosByBeamBunchByUse)[locUnknownComboUse];
1891 
1892  //Set the resume indices
1893  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
1894  return;
1895  }
1896 
1897  //if on the all-showers stage, first copy over ALL fcal-only results
1898  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, Get_SourceComboVectorResource());
1899  if(locComboingStage == d_MixedStage)
1900  Copy_ZIndependentMixedResults(locComboUseToCreate, locChargedCombo_Presiding);
1901 
1902  if(locSourceCombos->empty())
1903  return; //nothing to create
1904 
1905  if((locComboingStage == d_MixedStage) && (locVertexZBin == DSourceComboInfo::Get_VertexZIndex_Unknown()))
1906  {
1907  //we need a zbin for BCAL showers, but it is unknown: can't cut yet!
1908  //However, the FCAL ones were already cut during the z-independent stage, and have already been saved
1909  //so, just copy over the bcal results from the unknown use
1910  for(auto& locCombo : *locSourceCombos)
1911  {
1912  if(!locCombo->Get_IsComboingZIndependent()) //bcal only
1913  locSourceCombosByUseSoFar[locComboUseToCreate]->push_back(locCombo);
1914  }
1915 
1916  //Copy over the combos-by-beam-bunch
1917  auto& locUnknownBothCombosByBeamBunch = (*locSourceCombosByBeamBunchByUse)[locUnknownComboUse];
1918  for(const auto& locComboBeamBunchPair : locUnknownBothCombosByBeamBunch)
1919  {
1920  if(locComboBeamBunchPair.first.size() > 1)
1921  continue; //don't copy the overlap ones: they are not complete & need to be filled on the fly
1922  for(const auto& locCombo : locComboBeamBunchPair.second)
1923  {
1924  if(!locCombo->Get_IsComboingZIndependent()) //bcal only
1925  (*locSourceCombosByBeamBunchByUse)[locComboUseToCreate][locComboBeamBunchPair.first].push_back(locCombo);
1926  }
1927  }
1928 
1929  //Set the resume indices
1930  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
1931  return;
1932  }
1933 
1934  //place an invariant mass cut & save the results
1935  auto locTargetPIDToSubtract = std::get<4>(locComboUseToCreate);
1936  for(const auto& locSourceCombo : *locSourceCombos)
1937  {
1938  //If on all-showers stage, and combo is fcal-only, don't save (combo already created!!)
1939  if((locComboingStage == d_MixedStage) && locSourceCombo->Get_IsComboingZIndependent())
1940  continue; //this combo has already passed the cut & been saved: during the FCAL-only stage
1941  if(!dSourceComboP4Handler->Cut_InvariantMass_NoMassiveNeutrals(locSourceCombo, locDecayPID, locTargetPIDToSubtract, dTargetCenter, locVertexZBin, false))
1942  continue; //vertex not used if accurate-flag is false: can be anything (target center)
1943 
1944  //save the results
1945  locSourceCombosByUseSoFar[locComboUseToCreate]->push_back(locSourceCombo);
1946  if(locComboingStage == d_ChargedStage)
1947  continue;
1948 
1949  //register beam bunches
1950  const auto& locBeamBunches = dValidRFBunches_ByCombo[std::make_pair(locSourceCombo, locVertexZBin)];
1951  Register_ValidRFBunches(locComboUseToCreate, locSourceCombo, locBeamBunches, locComboingStage, locChargedCombo_Presiding);
1952  }
1953 
1954  //Set the resume indices
1955  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
1956 }
1957 
1958 void DSourceComboer::Create_SourceCombos_Unknown(const DSourceComboUse& locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, unsigned char locNumTabs)
1959 {
1960  /****************************************************** COMBOING PARTICLES *****************************************************
1961  *
1962  * First combo VERTICALLY, and then HORIZONTALLY
1963  * What does this mean?
1964  * Vertically: Make combos of size N of each PID needed (e.g. 3 pi0s)
1965  * Horizontally: Make combos of different PIDs (e.g. 2pi0, pi+, pi-, p)
1966  *
1967  * Why start with vertical comboing?
1968  * because the thing that takes the most time is when someone decides to analyze (e.g.) 2pi0, 3pi0, then 3pi0 eta, 3pi0 something else, 4pi0, etc.
1969  * we want to make the Npi0 combos as needed, then reuse the Npi0s when making combos of other types
1970  * thus we want to build vertically (pi0s together, then etas together), and THEN horizontally (combine pi0s & etas, etc)
1971  * plus, when building vertically, it's easier to keep track of things since the PID / decay-parent is the same
1972  *
1973  * Build all possible combos for all NEEDED GROUPINGS for each of the FURTHER DECAYS (if not done already)
1974  * this becomes a series of recursive calls
1975  * e.g. if need 3 pi0s, call for 2pi0s, which calls for 1pi0, which calls for 2g
1976  * then do the actual pi0 groupings on the return
1977  *
1978  * Note, if we combo vertically (e.g. 3pi0s, 2pi+'s, etc.), they are created with a use that is strictly that content.
1979  * Then, when we combo them horizontally, they are promoted out of the vertical combo, at the same level as everything else in the new horizontal combo.
1980  * This reduces the depth-complexity of the combos.
1981  *
1982  *******************************************************************************************************************************/
1983 
1984  if(dDebugLevel > 0)
1985  {
1986  cout << endl;
1987  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1988  cout << "Create_SourceCombos_Unknown: Stage, presiding charged combo, use to create = " << locComboingStage << ", " << locChargedCombo_Presiding << endl;
1989  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1990  cout << "PRESIDING COMBO:" << endl;
1991  DAnalysis::Print_SourceCombo(locChargedCombo_Presiding, locNumTabs);
1992  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
1993  cout << "USE TO CREATE:" << endl;
1994  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
1995  }
1996 
1997  //get use info, combos
1998  auto locComboInfoToCreate = std::get<2>(locComboUseToCreate);
1999  auto locChargeContent = dComboInfoChargeContent[locComboInfoToCreate];
2000  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locChargeContent, locChargedCombo_Presiding);
2001 
2002  Combo_Vertically_AllDecays(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding, locNumTabs);
2003  if(locSourceCombosByUseSoFar.find(locComboUseToCreate) != locSourceCombosByUseSoFar.end())
2004  {
2005  if(dDebugLevel > 0)
2006  cout << "We're done!" << endl;
2007  return; //we're done!
2008  }
2009 
2010  if((locComboingStage == d_ChargedStage) || (locChargeContent == d_Neutral))
2011  {
2012  Combo_Vertically_AllParticles(locComboUseToCreate, locComboingStage, locNumTabs); //no such thing as a "mixed" particle
2013  if(locSourceCombosByUseSoFar.find(locComboUseToCreate) != locSourceCombosByUseSoFar.end())
2014  {
2015  if(dDebugLevel > 0)
2016  cout << "We're done!" << endl;
2017  return; //we're done!
2018  }
2019  }
2020 
2021  //OK, now build horizontally!! //group particles with different PIDs
2022  Combo_Horizontally_All(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding, locNumTabs);
2023 }
2024 
2025 /************************************************************** BUILD PHOTON COMBOS - VERTICALLY ****************************************************************/
2026 
2027 void DSourceComboer::Combo_Vertically_AllDecays(const DSourceComboUse& locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, unsigned char locNumTabs)
2028 {
2029  if(dDebugLevel >= 5)
2030  {
2031  cout << endl;
2032  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2033  cout << "Combo_Vertically_AllDecays: Stage, presiding charged combo = " << locComboingStage << ", " << locChargedCombo_Presiding << endl;
2034  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2035  cout << "PRESIDING COMBO:" << endl;
2036  DAnalysis::Print_SourceCombo(locChargedCombo_Presiding, locNumTabs);
2037  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2038  cout << "USE TO CREATE:" << endl;
2039  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
2040  }
2041 
2042  //get combo use contents
2043  auto locComboInfo = std::get<2>(locComboUseToCreate);
2044  auto locVertexZBin = std::get<1>(locComboUseToCreate);
2045  auto locNumParticlesNeeded = locComboInfo->Get_NumParticles();
2046  auto locFurtherDecays = locComboInfo->Get_FurtherDecays();
2047 
2048  //for each further decay map entry (e.g. pi0, 3), this is a collection of the uses representing those groupings //e.g. Unknown -> 3pi0
2049  for(const auto& locFurtherDecayPair : locFurtherDecays)
2050  {
2051  auto& locSourceComboDecayUse = locFurtherDecayPair.first; //e.g. pi0, -> 2g
2052  auto& locNumDecaysNeeded = locFurtherDecayPair.second; //N of the above decay (e.g. pi0s)
2053  auto locSourceComboDecayInfo = std::get<2>(locSourceComboDecayUse);
2054  auto locDecayChargeContent = dComboInfoChargeContent[locSourceComboDecayInfo];
2055 
2056  if((locComboingStage == d_ChargedStage) && (locDecayChargeContent == d_Neutral))
2057  continue; //skip for now!!
2058  //if on a mixed stage, and the to-build combo info is fully charged, skip it: it's already been done
2059  if((locComboingStage != d_ChargedStage) && (locDecayChargeContent == d_Charged))
2060  continue;
2061 
2062  if(locNumDecaysNeeded == 1)
2063  {
2064  //must dive down to get the next charged combo
2065  //building for the first time: the first one (later ones will be grabbed when building these combos vertically (in Combo_Vertically_NDecays))
2066  auto locChargedCombo_NextPresiding = Get_NextChargedCombo(locChargedCombo_Presiding, locSourceComboDecayUse, locComboingStage, true, 1);
2067  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locDecayChargeContent, locChargedCombo_NextPresiding);
2068 
2069  //build the decay combo directly
2070  if(locSourceCombosByUseSoFar.find(locSourceComboDecayUse) == locSourceCombosByUseSoFar.end()) //if not done already!
2071  {
2072  //must return to top-level combo function to build this decay, as this may have any structure
2073  Create_SourceCombos(locSourceComboDecayUse, locComboingStage, locChargedCombo_NextPresiding, locNumTabs + 1);
2074  }
2075  else if(dDebugLevel > 0)
2076  cout << "This decay already created!" << endl;
2077 
2078  continue; //no actual comboing needs to be done with this just yet, will do so in horizontal stage
2079  }
2080 
2081  //OK, so we need a grouping of N > 1 decays (e.g. pi0s)
2082  //so, let's create a use of Unknown -> N pi0s (e.g.)
2083  //if we can just utilize the use from the input combo-info, then we will. if not, we'll make a new one
2084  auto locNeededGroupingUse = locComboUseToCreate;
2085  if((locFurtherDecays.size() > 1) || !locNumParticlesNeeded.empty()) //if true: can't use the input
2086  {
2087  auto locGroupingComboInfo = GetOrMake_SourceComboInfo({}, {std::make_pair(locSourceComboDecayUse, locNumDecaysNeeded)}, locNumTabs); // -> N pi0s (e.g.)
2088  locNeededGroupingUse = std::make_tuple(Unknown, locVertexZBin, locGroupingComboInfo, false, Unknown); // Unknown -> Npi0s (e.g.)
2089  }
2090 
2091  // Now, see whether the combos for this grouping have already been done
2092  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locDecayChargeContent, locChargedCombo_Presiding);
2093  if(locSourceCombosByUseSoFar.find(locNeededGroupingUse) != locSourceCombosByUseSoFar.end())
2094  {
2095  if(dDebugLevel > 0)
2096  cout << "This decay already created!" << endl;
2097  continue; //it's already done!!
2098  }
2099 
2100  //it's not already done. darn it.
2101  //build an info and a use for a direct grouping of N - 1 decays //e.g. 2pi0s
2102  auto locNMinus1ComboUse = locSourceComboDecayUse; //initialize (is valid if #needed == 2, otherwise will create it)
2103  if(locNumDecaysNeeded > 2)
2104  {
2105  auto locNMinus1Info = GetOrMake_SourceComboInfo({}, {std::make_pair(locSourceComboDecayUse, locNumDecaysNeeded - 1)}, locNumTabs); // 0 detected particles, N - 1 pi0s (e.g.)
2106  locNMinus1ComboUse = std::make_tuple(Unknown, locVertexZBin, locNMinus1Info, false, Unknown); // Unknown -> N - 1 pi0s (e.g.)
2107  }
2108 
2109  // Now, see whether the combos for the direct N - 1 grouping have already been done. If not, create them
2110  if(locNumDecaysNeeded == 2) //(so N - 1 = 1)
2111  {
2112  auto locChargedCombo_NextPresiding = Get_NextChargedCombo(locChargedCombo_Presiding, locSourceComboDecayUse, locComboingStage, true, 1);
2113  locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locDecayChargeContent, locChargedCombo_NextPresiding);
2114  if(locSourceCombosByUseSoFar.find(locNMinus1ComboUse) == locSourceCombosByUseSoFar.end())
2115  Create_SourceCombos(locNMinus1ComboUse, locComboingStage, locChargedCombo_NextPresiding, locNumTabs + 1);
2116  }
2117  else
2118  {
2119  //if not done yet, no need to go to top-level combo function since just N - 1: can re-call this one
2120  if(locSourceCombosByUseSoFar.find(locNMinus1ComboUse) == locSourceCombosByUseSoFar.end())
2121  Combo_Vertically_AllDecays(locNMinus1ComboUse, locComboingStage, locChargedCombo_Presiding, locNumTabs + 1);
2122  }
2123 
2124  //Finally, we can actually DO the grouping, between the N - 1 combos and the one-off combos
2125  Combo_Vertically_NDecays(locNeededGroupingUse, locNMinus1ComboUse, locSourceComboDecayUse, locComboingStage, locChargedCombo_Presiding, locNumTabs);
2126  }
2127 }
2128 
2129 void DSourceComboer::Combo_Vertically_NDecays(const DSourceComboUse& locComboUseToCreate, const DSourceComboUse& locNMinus1ComboUse, const DSourceComboUse& locSourceComboDecayUse, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, unsigned char locNumTabs)
2130 {
2131  if(dDebugLevel >= 5)
2132  {
2133  cout << endl;
2134  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2135  cout << "Combo_Vertically_NDecays: Stage, presiding charged combo = " << locComboingStage << ", " << locChargedCombo_Presiding << endl;
2136  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2137  cout << "PRESIDING COMBO:" << endl;
2138  DAnalysis::Print_SourceCombo(locChargedCombo_Presiding, locNumTabs);
2139  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2140  cout << "USE TO CREATE:" << endl;
2141  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
2142  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2143  cout << "N - 1 USE:" << endl;
2144  DAnalysis::Print_SourceComboUse(locNMinus1ComboUse, locNumTabs);
2145  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2146  cout << "SINGLE USE:" << endl;
2147  DAnalysis::Print_SourceComboUse(locSourceComboDecayUse, locNumTabs);
2148  }
2149 
2150  auto locNIs2Flag = (locNMinus1ComboUse == locSourceComboDecayUse); //true if need exactly 2 decaying particles
2151 
2152  //Get combos so far
2153  auto locComboInfoToCreate = std::get<2>(locComboUseToCreate);
2154  auto locChargeContent = dComboInfoChargeContent[locComboInfoToCreate];
2155  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locChargeContent, locChargedCombo_Presiding);
2156 
2157  //e.g. we are grouping 1 pi0 with N - 1 pi0s to make a combo of N pi0s
2158  //so, let's get the combos for (e.g.) 1 pi0 and for N - 1 pi0s
2159  const auto& locCombos_NMinus1 = *locSourceCombosByUseSoFar[locNMinus1ComboUse]; //Combos are a vector of (e.g.): -> N - 1 pi0s
2160 
2161  //if on the all-showers stage, first copy over ALL fcal-only results
2162  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, Get_SourceComboVectorResource());
2163  if(locComboingStage == d_MixedStage)
2164  Copy_ZIndependentMixedResults(locComboUseToCreate, locChargedCombo_Presiding);
2165 
2166  if(dDebugLevel >= 20)
2167  cout << "n-1 size: " << locCombos_NMinus1.size() << endl;
2168  if(locCombos_NMinus1.empty())
2169  return; //nothing to create
2170 
2171  //if comboing N mixed combos (locComboUseToCreate) (which are thus all used in the same step), do this:
2172  //locChargedCombo_WithNow corresponds to N mixed combos
2173  auto locZIndependentDecayUse = Get_ZIndependentUse(locSourceComboDecayUse);
2174  auto locInstanceUse = locCombos_NMinus1.front()->Get_IsComboingZIndependent() ? locZIndependentDecayUse : locSourceComboDecayUse;
2175  auto locFirstNMinus1FurtherDecayCombos = locCombos_NMinus1.front()->Get_FurtherDecayCombos();
2176  size_t locInstance = 2; //changed below if needed
2177  if(!locNIs2Flag)
2178  {
2179  auto locIteratorPair = std::equal_range(locFirstNMinus1FurtherDecayCombos.begin(), locFirstNMinus1FurtherDecayCombos.end(), locInstanceUse, DSourceCombo::DCompare_FurtherDecays());
2180  locInstance = (*locIteratorPair.first).second.size() + 1; //numbering starts with 1, not 0
2181  }
2182  auto locPreviousPresidingCombo = Get_NextChargedCombo(locChargedCombo_Presiding, locSourceComboDecayUse, locComboingStage, true, locInstance);
2183  auto locChargedCombo_WithPrevious = Get_ChargedCombo_WithNow(locPreviousPresidingCombo, locComboInfoToCreate, locComboingStage);
2184 
2185  if(dDebugLevel >= 20)
2186  cout << "instance = " << locInstance << ", combos: presiding, previous-presiding, with-previous: " << locChargedCombo_Presiding << ", " << locPreviousPresidingCombo << ", " << locChargedCombo_WithPrevious << endl;
2187 
2188  //now, for each combo of N - 1 (e.g.) pi0s, see which of the single-decay combos are a valid grouping
2189  //valid grouping:
2190  //TEST 1: If (e.g.) pi0s have names "A", "B", "C", don't include the grouping "ABA", and don't include "ACB" if we already have "ABC"
2191  //TEST 2: Also, don't re-use a shower we've already used (e.g. if A & C each contain the same photon, don't group them together)
2192  //Technically, if we pass Test 2 we automatically pass Test 1.
2193  //However, validating for Test 1 is much faster, as discussed below.
2194  auto& locComboResultsVector = *(locSourceCombosByUseSoFar[locComboUseToCreate]);
2195  for(const auto& locCombo_NMinus1 : locCombos_NMinus1)
2196  {
2197  //loop over potential combos to add to the group, creating a new combo for each valid (non-duplicate) grouping
2198  //however, we don't have to loop over all of the combos!!
2199 
2200  //first of all, get the potential combos that satisfy the RF bunches for the N - 1 combo
2201  const auto& locValidRFBunches_NMinus1 = dValidRFBunches_ByCombo[std::make_pair(locCombo_NMinus1, std::get<1>(locNMinus1ComboUse))];
2202  const auto& locDecayCombos_1 = Get_CombosForComboing(locSourceComboDecayUse, locComboingStage, locValidRFBunches_NMinus1, locPreviousPresidingCombo);
2203  if(dDebugLevel >= 20)
2204  cout << "decay combos vector address, size: " << &locDecayCombos_1 << ", " << locDecayCombos_1.size() << endl;
2205  if(dDebugLevel >= 100)
2206  {
2207  cout << "Vector combos: ";
2208  for(auto& locCombo : locDecayCombos_1)
2209  cout << locCombo << ", ";
2210  cout << endl;
2211  }
2212 
2213  //now, note that all of the combos are stored in the order in which they were created (e.g. A, B, C, D)
2214  //so (e.g.), groupings of 2 will be created and saved in the order: AB, AC, AD, BC, BD, CD
2215  //above, on the B-loop, we start the search at "C," not at A, because this was already tested on an earlier pass
2216  //therefore, start the search one AFTER the LAST (e.g. -> 2 photon) combo of the N - 1 group
2217  //this will guarantee we pass "TEST 1" without ever checking
2218 
2219  //actually, we already saved the iterator to the first (e.g.) pi0 to test when we saved the N - 1 combo, so just retrieve it
2220  auto locNMinus1ComboDecayUse = locCombo_NMinus1->Get_IsComboingZIndependent() ? locZIndependentDecayUse : locSourceComboDecayUse;
2221  auto locNMinus1FurtherDecayCombos = locCombo_NMinus1->Get_FurtherDecayCombos();
2222  auto locNMinus1LastCombo = locCombo_NMinus1;
2223  if(!locNIs2Flag)
2224  {
2225  auto locIteratorPair = std::equal_range(locNMinus1FurtherDecayCombos.begin(), locNMinus1FurtherDecayCombos.end(), locNMinus1ComboDecayUse, DSourceCombo::DCompare_FurtherDecays());
2226  locNMinus1LastCombo = (*locIteratorPair.first).second.back();
2227  }
2228 
2229  auto locComboSearchIndex = Get_ResumeAtIndex_Combos(locSourceComboDecayUse, locNMinus1LastCombo, locValidRFBunches_NMinus1, locComboingStage);
2230  if(dDebugLevel >= 20)
2231  cout << "n-1 last combo, begin search index = : " << locNMinus1LastCombo << ", " << locComboSearchIndex << endl;
2232  if(locComboSearchIndex == locDecayCombos_1.size())
2233  continue; //e.g. this combo is "AD" and there are only 4 reconstructed combos (ABCD): no potential matches! move on to the next N - 1 combo
2234 
2235  //before we loop, first get all of the showers used to make the N - 1 grouping, and sort it so that we can quickly search it
2236  auto locUsedParticles_NMinus1 = DAnalysis::Get_SourceParticles(locCombo_NMinus1->Get_SourceParticles(true)); //true: entire chain
2237  std::sort(locUsedParticles_NMinus1.begin(), locUsedParticles_NMinus1.end()); //must sort, because when retrieving entire chain is unsorted
2238 
2239  //this function will do our "TEST 2"
2240  auto Search_Duplicates = [&locUsedParticles_NMinus1](const JObject* locParticle) -> bool
2241  {return std::binary_search(locUsedParticles_NMinus1.begin(), locUsedParticles_NMinus1.end(), locParticle);};
2242 
2243  auto locIsZIndependent_NMinus1 = locCombo_NMinus1->Get_IsComboingZIndependent();
2244 
2245  //now loop over the potential combos
2246  for(; locComboSearchIndex != locDecayCombos_1.size(); ++locComboSearchIndex)
2247  {
2248  const auto locDecayCombo_1 = locDecayCombos_1[locComboSearchIndex];
2249 
2250  //If on all-showers stage, and combo is fcal-only, don't save (combo already created!!)
2251  auto locIsZIndependent = locIsZIndependent_NMinus1 && locDecayCombo_1->Get_IsComboingZIndependent();
2252  if((locComboingStage == d_MixedStage) && locIsZIndependent)
2253  continue; //this combo has already been created (assuming it was valid): during the FCAL-only stage
2254 
2255  //conduct "TEST 2" search: search the N - 1 shower vector to see if any of the showers in this combo are duplicated
2256  auto locUsedParticles_1 = DAnalysis::Get_SourceParticles(locDecayCombo_1->Get_SourceParticles(true)); //true: entire chain
2257  if(std::any_of(locUsedParticles_1.begin(), locUsedParticles_1.end(), Search_Duplicates))
2258  continue; //at least one photon was a duplicate, this combo won't work
2259 
2260  //no duplicates: this combo is unique. build a new combo!
2261 
2262  //See which RF bunches match up //guaranteed to be at least one, due to selection in Get_ParticlesForComboing() function
2263  auto locValidRFBunches = dSourceComboTimeHandler->Get_CommonRFBunches(locValidRFBunches_NMinus1, dValidRFBunches_ByCombo[std::make_pair(locDecayCombo_1, std::get<1>(locSourceComboDecayUse))]);
2264 
2265  //Combine the decay combos
2266  vector<const DSourceCombo*> locAllDecayCombos;
2267  if(locNIs2Flag) //N = 2 Two identical combos (e.g. 2 of pi0 -> 2g)
2268  locAllDecayCombos = {locCombo_NMinus1, locDecayCombo_1};
2269  else //combine a combo of N - 1 (e.g. pi0) decays to this new one
2270  {
2271  auto locIteratorPair = std::equal_range(locNMinus1FurtherDecayCombos.begin(), locNMinus1FurtherDecayCombos.end(), locNMinus1ComboDecayUse, DSourceCombo::DCompare_FurtherDecays());
2272  locAllDecayCombos = (*locIteratorPair.first).second;
2273  locAllDecayCombos.push_back(locDecayCombo_1);
2274  }
2275 
2276  //then create the new combo
2277  DSourceCombosByUse_Small locFurtherDecayCombos = {std::make_pair(locSourceComboDecayUse, locAllDecayCombos)}; //arguments (e.g.): (pi0, -> 2g), N combos of: -> 2g
2278  auto locCombo = Get_SourceComboResource();
2279  locCombo->Set_Members({}, locFurtherDecayCombos, locIsZIndependent); // 1 combo of N (e.g.) pi0s
2280  if(dDebugLevel >= 10)
2281  {
2282  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2283  cout << "CREATED COMBO:" << endl;
2284  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
2285  }
2286 
2287  //save it!
2288  locComboResultsVector.push_back(locCombo);
2289  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, locChargedCombo_Presiding);
2290  }
2291  }
2292  if((dDebugLevel > 0) || (dDebugLevel == -1))
2293  Check_ForDuplicates(locComboResultsVector);
2294 
2295  //Set the resume indices
2296  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
2297  if(dDebugLevel >= 5)
2298  {
2299  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2300  cout << "Combo_Vertically_NDecays: NUM SOURCE COMBOS CREATED: " << locComboResultsVector.size() << endl;
2301  }
2302 }
2303 
2304 void DSourceComboer::Combo_Vertically_AllParticles(const DSourceComboUse& locComboUseToCreate, ComboingStage_t locComboingStage, unsigned char locNumTabs)
2305 {
2306  if(dDebugLevel >= 5)
2307  {
2308  cout << endl;
2309  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2310  cout << "Combo_Vertically_AllParticles: Stage = " << locComboingStage << endl;
2311  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2312  cout << "USE TO CREATE:" << endl;
2313  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
2314  }
2315 
2316  //get combo use contents
2317  auto locVertexZBin = std::get<1>(locComboUseToCreate);
2318  auto locNumParticlesNeeded = std::get<2>(locComboUseToCreate)->Get_NumParticles();
2319  auto locFurtherDecays = std::get<2>(locComboUseToCreate)->Get_FurtherDecays();
2320 
2321  //Get combos so far //guaranteed not to be mixed
2322  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, d_Neutral); //if not neutral then is on charged stage: argument doesn't matter
2323 
2324  //for each further decay map entry (e.g. pi0, 3), this is a collection of the uses representing those groupings //e.g. Unknown -> 3pi0
2325  for(const auto& locParticlePair : locNumParticlesNeeded)
2326  {
2327  //get PID information
2328  auto& locPID = locParticlePair.first; //e.g. pi0, -> 2g
2329  auto& locNumPIDNeeded = locParticlePair.second; //N of the above decay (e.g. pi0s)
2330 
2331  if(locNumPIDNeeded == 1)
2332  continue; //nothing to do vertically; we will combo this horizontally later
2333 
2334  if((locComboingStage == d_ChargedStage) && (ParticleCharge(locPID) == 0))
2335  continue; //skip for now!!
2336  if((locComboingStage != d_ChargedStage) && (ParticleCharge(locPID) != 0))
2337  continue; //already done!
2338 
2339  //OK, so we need a grouping of N > 1 particles with the same PID (e.g. g's)
2340  //so, let's create a use of Unknown -> N g's (e.g.)
2341  //if we can just utilize the use from the input combo-info, then we will. if not, we'll make a new one
2342  DSourceComboUse locNeededGroupingUse = locComboUseToCreate;
2343  if((locNumParticlesNeeded.size() > 1) || !locFurtherDecays.empty()) //if true: can't use the input
2344  {
2345  auto locGroupingComboInfo = GetOrMake_SourceComboInfo({std::make_pair(locPID, locNumPIDNeeded)}, {}, locNumTabs); // -> N g's (e.g.)
2346  locNeededGroupingUse = std::make_tuple(Unknown, locVertexZBin, locGroupingComboInfo, false, Unknown); // Unknown -> N g's (e.g.)
2347  }
2348 
2349  //See whether the combos for this grouping have already been done
2350  if(locSourceCombosByUseSoFar.find(locNeededGroupingUse) != locSourceCombosByUseSoFar.end())
2351  {
2352  if(dDebugLevel > 0)
2353  cout << "This group already created!" << endl;
2354  continue; //it's already done!!
2355  }
2356 
2357  //it's not already done. darn it.
2358  //if it's a direct combo of 2 particles, just make it and continue
2359  if(locNumPIDNeeded == 2)
2360  {
2361  Combo_Vertically_NParticles(locNeededGroupingUse, DSourceComboUse(), locComboingStage, locNumTabs);
2362  continue;
2363  }
2364 
2365  //build an info and a use for a direct grouping of N - 1 particles //e.g. 3 g's
2366  auto locNMinus1Info = GetOrMake_SourceComboInfo({std::make_pair(locPID, locNumPIDNeeded - 1)}, {}, locNumTabs); // N - 1 g's (e.g.), no decaying particles
2367  DSourceComboUse locNMinus1ComboUse(Unknown, locVertexZBin, locNMinus1Info, false, Unknown); // Unknown -> N - 1 g's (e.g.)
2368 
2369  // Now, see whether the combos for the direct N - 1 grouping have already been done. If not, create them
2370  if(locSourceCombosByUseSoFar.find(locNMinus1ComboUse) == locSourceCombosByUseSoFar.end())
2371  Combo_Vertically_AllParticles(locNMinus1ComboUse, locComboingStage, locNumTabs + 1); //no need to go to top-level combo function since just N - 1: can re-call this one
2372 
2373  //Finally, we can actually DO the grouping, between the N - 1 particles and one more particle
2374  Combo_Vertically_NParticles(locNeededGroupingUse, locNMinus1ComboUse, locComboingStage, locNumTabs);
2375  }
2376 }
2377 
2378 void DSourceComboer::Combo_Vertically_NParticles(const DSourceComboUse& locComboUseToCreate, const DSourceComboUse& locNMinus1ComboUse, ComboingStage_t locComboingStage, unsigned char locNumTabs)
2379 {
2380  if(dDebugLevel >= 5)
2381  {
2382  cout << endl;
2383  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2384  cout << "Combo_Vertically_NParticles: Stage = " << locComboingStage << endl;
2385  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2386  cout << "USE TO CREATE:" << endl;
2387  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
2388  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2389  cout << "N - 1 USE:" << endl;
2390  DAnalysis::Print_SourceComboUse(locNMinus1ComboUse, locNumTabs);
2391  }
2392 
2393  //either: combining two particles with the same PID to create a new combo, or combining a combo of N particles (with the same PID) with one more particle
2394  auto locComboInfo = std::get<2>(locComboUseToCreate);
2395  auto locParticlePair = locComboInfo->Get_NumParticles().back(); //is guaranteed to be size 1
2396  auto locPID = locParticlePair.first;
2397  auto locNumParticles = locParticlePair.second;
2398  auto locVertexZBin = std::get<1>(locComboUseToCreate);
2399 
2400  //Get combos so far //guaranteed not to be mixed
2401  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, d_Neutral); //if not neutral then is on charged stage: argument doesn't matter
2402 
2403  //check if comboing massive neutrals: comboing results are independent of z-bin: see if this has already been done for a different (non-z-independent) zbin
2404  if((ParticleCharge(locPID) == 0) && (ParticleMass(locPID) > 0.0) && (locComboingStage == d_MixedStage))
2405  {
2406  //if so, just copy the results into this zbin!
2407  //note that this can only the case for comboing particles, not comboing combos: decays contain uses, which have the zbin specified
2408  for(signed char locZBin = 0; locZBin < (signed char)dSourceComboTimeHandler->Get_NumVertexZBins(); ++locZBin)
2409  {
2410  if(locZBin == locVertexZBin)
2411  continue;
2412  auto locZBinUse = DSourceComboUse{Unknown, locZBin, locComboInfo, false, Unknown};
2413  if(locSourceCombosByUseSoFar.find(locZBinUse) == locSourceCombosByUseSoFar.end())
2414  continue;
2415 
2416  //just copy the results!
2417  locSourceCombosByUseSoFar[locComboUseToCreate] = locSourceCombosByUseSoFar[locZBinUse];
2418  auto& locSourceCombosByBeamBunchByUse = Get_SourceCombosByBeamBunchByUse(dComboInfoChargeContent[locComboInfo], nullptr);
2419  locSourceCombosByBeamBunchByUse.emplace(locComboUseToCreate, locSourceCombosByBeamBunchByUse[locZBinUse]);
2420 
2421  //Set the resume indices
2422  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
2423  if(dDebugLevel >= 5)
2424  {
2425  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2426  cout << "Combo_Vertically_NParticles: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseSoFar[locComboUseToCreate]->size() << endl;
2427  }
2428  return;
2429  }
2430  }
2431 
2432  //if on the all-showers stage, first copy over ALL fcal-only results
2433  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, Get_SourceComboVectorResource());
2434  if(locComboingStage == d_MixedStage)
2435  Copy_ZIndependentMixedResults(locComboUseToCreate, nullptr);
2436 
2437  if(locNumParticles == 2)
2438  {
2439  //Get particles for comboing
2440  const auto& locParticles = Get_ParticlesForComboing(locPID, locComboingStage, {}, locVertexZBin);
2441  if(locParticles.size() < 2)
2442  return; //not enough to create combos
2443 
2444  auto locLastIteratorToCheck = std::prev(locParticles.end());
2445  for(auto locFirstIterator = locParticles.begin(); locFirstIterator != locLastIteratorToCheck; ++locFirstIterator)
2446  {
2447  auto locRFBunches_First = (locPID == Gamma) ? dSourceComboTimeHandler->Get_ValidRFBunches(*locFirstIterator, locVertexZBin) : vector<int>{};
2448  for(auto locSecondIterator = std::next(locFirstIterator); locSecondIterator != locParticles.end(); ++locSecondIterator)
2449  {
2450  auto locIsZIndependent = (locComboingStage == d_MixedStage_ZIndependent) || (Get_IsComboingZIndependent(*locFirstIterator, locPID) && Get_IsComboingZIndependent(*locSecondIterator, locPID));
2451  if((locComboingStage == d_MixedStage) && locIsZIndependent)
2452  continue; //this combo has already been created (assuming it was valid): during the FCAL-only stage
2453 
2454  //See which RF bunches match up, if any //if charged or massive neutrals, ignore (they don't choose at this stage)
2455  auto locValidRFBunches = (locPID != Gamma) ? vector<int>{} : dSourceComboTimeHandler->Get_CommonRFBunches(locRFBunches_First, *locSecondIterator, locVertexZBin);
2456  if((locPID == Gamma) && locValidRFBunches.empty())
2457  continue;
2458 
2459  //check if this combo is unique!
2460  //it will not be unique if: comboing photons, on mixed stage, and this combo has already been created for a different zbin
2461  //if so, don't duplicate memory
2462  //note that this can only the case for comboing particles, not comboing combos: decays contain uses, which have the zbin specified
2463  vector<const JObject*> locPhotons;
2464  if((locPID == Gamma) && (locComboingStage == d_MixedStage))
2465  {
2466  locPhotons = {*locFirstIterator, *locSecondIterator};
2467  auto locIterator = dNPhotonsToComboMap.find(locPhotons);
2468  if(locIterator != dNPhotonsToComboMap.end()) //if true, already exists, copy it and continue;
2469  {
2470  auto locCombo = locIterator->second;
2471  locSourceCombosByUseSoFar[locComboUseToCreate]->push_back(locCombo);
2472  if(dDebugLevel >= 15)
2473  {
2474  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2475  cout << "COPIED COMBO:" << endl;
2476  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
2477  }
2478  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, nullptr);
2479  continue;
2480  }
2481  }
2482 
2483  auto locCombo = Get_SourceComboResource();
2484  locCombo->Set_Members({std::make_pair(locPID, *locFirstIterator), std::make_pair(locPID, *locSecondIterator)}, {}, locIsZIndependent);
2485  locSourceCombosByUseSoFar[locComboUseToCreate]->push_back(locCombo); //save it //in creation order
2486  if(locPID == Gamma)
2487  dNPhotonsToComboMap.emplace(locPhotons, locCombo);
2488  if(dDebugLevel >= 10)
2489  {
2490  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2491  cout << "CREATED COMBO:" << endl;
2492  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
2493  }
2494 
2495  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, nullptr);
2496  }
2497  }
2498  if((dDebugLevel > 0) || (dDebugLevel == -1))
2499  Check_ForDuplicates(*(locSourceCombosByUseSoFar[locComboUseToCreate]));
2500 
2501  //Set the resume indices
2502  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
2503  if(dDebugLevel >= 5)
2504  {
2505  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2506  cout << "Combo_Vertically_NParticles: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseSoFar[locComboUseToCreate]->size() << endl;
2507  }
2508  return;
2509  }
2510 
2511  //create combo of N same-PID-particles by adding one particle to previously-created combos of N - 1 same-PID-particles
2512  const auto& locCombos_NMinus1 = *locSourceCombosByUseSoFar[locNMinus1ComboUse]; //Each combo contains a vector of N - 1 same-PID-particles
2513  for(const auto& locCombo_NMinus1 : locCombos_NMinus1)
2514  {
2515  //Get particles for comboing
2516  const auto& locValidRFBunches_NMinus1 = dValidRFBunches_ByCombo[std::make_pair(locCombo_NMinus1, std::get<1>(locNMinus1ComboUse))];
2517  const auto& locParticles = Get_ParticlesForComboing(locPID, locComboingStage, locValidRFBunches_NMinus1, locVertexZBin);
2518 
2519  //retrieve where to begin the search
2520  auto locLastParticleInCombo = locCombo_NMinus1->Get_SourceParticles(false).back().second;
2521  auto locParticleSearchIndex = Get_ResumeAtIndex_Particles(locPID, locLastParticleInCombo, locValidRFBunches_NMinus1, locVertexZBin);
2522  if(dDebugLevel >= 20)
2523  cout << "particle index, #particles = " << locParticleSearchIndex << ", " << locParticles.size() << endl;
2524  if(locParticleSearchIndex == locParticles.size())
2525  continue; //e.g. this combo is "AD" and there are only 4 reconstructed combos (ABCD): no potential matches! move on to the next N - 1 combo
2526 
2527  auto locIsZIndependent_NMinus1 = locCombo_NMinus1->Get_IsComboingZIndependent();
2528 
2529  for(; locParticleSearchIndex < locParticles.size(); ++locParticleSearchIndex)
2530  {
2531  auto& locParticle = locParticles[locParticleSearchIndex];
2532  auto locIsZIndependent = (locComboingStage == d_MixedStage_ZIndependent) || (locIsZIndependent_NMinus1 && Get_IsComboingZIndependent(locParticle, locPID));
2533  if((locComboingStage == d_MixedStage) && locIsZIndependent)
2534  continue; //this combo has already been created (assuming it was valid): during the FCAL-only stage
2535 
2536  //See which RF bunches match up //guaranteed to be at least one, due to selection in Get_ParticlesForComboing() function
2537  //if charged or massive neutrals, ignore (they don't choose at this stage)
2538  auto locValidRFBunches = (locPID != Gamma) ? vector<int>{} : dSourceComboTimeHandler->Get_CommonRFBunches(locValidRFBunches_NMinus1, locParticle, locVertexZBin);
2539 
2540  auto locComboParticlePairs = locCombo_NMinus1->Get_SourceParticles();
2541 
2542  //check if this combo is unique!
2543  //it will not be unique if: comboing photons, on mixed stage, and this combo has already been created for a different zbin
2544  //if so, don't duplicate memory
2545  //note that this can only the case for comboing particles, not comboing combos: decays contain uses, which have the zbin specified
2546  vector<const JObject*> locPhotons;
2547  if((locPID == Gamma) && (locComboingStage == d_MixedStage))
2548  {
2549  auto GetParticle = [](const pair<Particle_t, const JObject*>& locPair) -> const JObject* {return locPair.second;};
2550  std::transform(locComboParticlePairs.begin(), locComboParticlePairs.end(), std::back_inserter(locPhotons), GetParticle);
2551  locPhotons.push_back(locParticle);
2552  auto locIterator = dNPhotonsToComboMap.find(locPhotons);
2553  if(locIterator != dNPhotonsToComboMap.end()) //if true, already exists, save it and continue;
2554  {
2555  auto locCombo = locIterator->second;
2556  locSourceCombosByUseSoFar[locComboUseToCreate]->push_back(locCombo);
2557  if(dDebugLevel >= 15)
2558  {
2559  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2560  cout << "COPIED COMBO:" << endl;
2561  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
2562  }
2563  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, nullptr);
2564  continue;
2565  }
2566  }
2567 
2568  locComboParticlePairs.emplace_back(locPID, locParticle);
2569  auto locCombo = Get_SourceComboResource();
2570  locCombo->Set_Members(locComboParticlePairs, {}, locIsZIndependent);
2571  locSourceCombosByUseSoFar[locComboUseToCreate]->push_back(locCombo); //save it //in creation order
2572  if(locPID == Gamma)
2573  dNPhotonsToComboMap.emplace(locPhotons, locCombo);
2574  if(dDebugLevel >= 10)
2575  {
2576  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2577  cout << "CREATED COMBO:" << endl;
2578  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
2579  }
2580 
2581  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, nullptr);
2582  }
2583  }
2584  if((dDebugLevel > 0) || (dDebugLevel == -1))
2585  Check_ForDuplicates(*(locSourceCombosByUseSoFar[locComboUseToCreate]));
2586 
2587  //Set the resume indices
2588  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
2589  if(dDebugLevel >= 5)
2590  {
2591  cout << "Combo_Vertically_NParticles: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseSoFar[locComboUseToCreate]->size() << endl;
2592  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2593  }
2594 }
2595 
2596 /************************************************************* BUILD PHOTON COMBOS - HORIZONTALLY ***************************************************************/
2597 
2598 void DSourceComboer::Combo_Horizontally_All(const DSourceComboUse& locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, unsigned char locNumTabs)
2599 {
2600  if(dDebugLevel >= 5)
2601  {
2602  cout << endl;
2603  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2604  cout << "Combo_Horizontally_All: Stage, presiding charged combo = " << locComboingStage << ", " << locChargedCombo_Presiding << endl;
2605  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2606  cout << "PRESIDING COMBO:" << endl;
2607  DAnalysis::Print_SourceCombo(locChargedCombo_Presiding, locNumTabs);
2608  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2609  cout << "USE TO CREATE:" << endl;
2610  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
2611  }
2612 
2613  //get combo use contents
2614  auto locVertexZBin = std::get<1>(locComboUseToCreate);
2615  const auto& locComboInfoToCreate = std::get<2>(locComboUseToCreate);
2616  auto locNumParticlesNeeded = locComboInfoToCreate->Get_NumParticles();
2617  auto locFurtherDecays = locComboInfoToCreate->Get_FurtherDecays();
2618 
2619  //first handle special cases:
2620  if(locNumParticlesNeeded.empty() && (locFurtherDecays.size() == 1))
2621  {
2622  Create_Combo_OneDecay(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding, locNumTabs);
2623  return;
2624  }
2625  if(locFurtherDecays.empty() && (locNumParticlesNeeded.size() == 1))
2626  {
2627  //we just need N (e.g.) photons together
2628  auto& locParticlePair = locNumParticlesNeeded.front();
2629  if(locParticlePair.second > 1)
2630  return; //already done when comboing vertically!!
2631 
2632  //not much of a combo if there's only 1, is it? //e.g. 1 charged track at a vertex
2633  if((locComboingStage == d_ChargedStage) && (ParticleCharge(locParticlePair.first) == 0))
2634  return; //skip for now!!
2635 
2636  Create_Combo_OneParticle(locComboUseToCreate, locComboingStage, locNumTabs);
2637  return;
2638  }
2639 
2640  //see if there is another combo that already exists that is a subset of what we requested
2641  //e.g. if we need a charged combo, a neutral combo, and a mixed: search for:
2642  //charged + neutral (no mixed)
2643  //charged + mixed (no neutral)
2644  //neutral + mixed (no charged)
2645  //e.g. if we need 2pi0s, one omega, and 1g: search for:
2646  //2pi0s, one omega: if exists, just combo that with 1g
2647  //2pi0s, one photon: if exists, just combo with one omega
2648  //etc.
2649 
2650  DSourceComboUse locComboUse_SubsetToBuild(Unknown, locVertexZBin, nullptr, false, Unknown);
2651  DSourceComboUse locComboUse_SubsetToAdd(Unknown, locVertexZBin, nullptr, false, Unknown);
2652  auto locChargedCombo_SubsetToBuildPresiding = locChargedCombo_Presiding;
2653 
2654  //First test the case: 1 set of particles, 1 decay
2655  bool locMissingSubsetIsDecayFlag = true; //set false if otherwise
2656  if((locFurtherDecays.size() == 1) && (locNumParticlesNeeded.size() == 1))
2657  {
2658  if(dDebugLevel >= 5)
2659  cout << "1 decay, 1 type of particle needed" << endl;
2660 
2661  //build the particles first, then we'll add the decay horizontally
2662  //unless the particles are fully neutral and we're on the charged stage: then we'll do the decay first
2663  if((locComboingStage == d_ChargedStage) && (ParticleCharge(locNumParticlesNeeded[0].first) == 0))
2664  {
2665  if(dDebugLevel >= 5)
2666  cout << "build decay first" << endl;
2667  auto locAllBut1ComboUse = locFurtherDecays[0].first; //do decay first
2668  locMissingSubsetIsDecayFlag = false;
2669 
2670  //Get combos so far
2671  auto locAllBut1ComboInfo = std::get<2>(locAllBut1ComboUse);
2672  auto locChargedCombo_NextPresiding = Get_NextChargedCombo(locChargedCombo_Presiding, locAllBut1ComboUse, locComboingStage, true, 1);
2673  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, dComboInfoChargeContent[locAllBut1ComboInfo], locChargedCombo_NextPresiding);
2674 
2675  // Now, see whether the combos for this grouping have already been done
2676  if(locSourceCombosByUseSoFar.find(locAllBut1ComboUse) == locSourceCombosByUseSoFar.end()) //if true: not yet
2677  locComboUse_SubsetToBuild = locAllBut1ComboUse; //will build below before the end of the function
2678  else //yes, it's already been done!
2679  {
2680  //and, since we are in the charged stage and the remaining particles are neutral, we are done for now
2681  //just copy the all-but-1 as the desired combos
2682  auto& locAllBut1Combos = locSourceCombosByUseSoFar[locAllBut1ComboUse];
2683  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, locAllBut1Combos);
2684  //Set the resume indices
2685  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_NextPresiding);
2686  return;
2687  }
2688  }
2689  else //particles first
2690  {
2691  if(dDebugLevel >= 5)
2692  cout << "build particles first" << endl;
2693 
2694  auto locAllBut1ComboInfo = GetOrMake_SourceComboInfo(locNumParticlesNeeded, {}, locNumTabs);
2695  auto locAllBut1ZBin = (Get_ChargeContent(locAllBut1ComboInfo) != d_Charged) ? locVertexZBin : DSourceComboInfo::Get_VertexZIndex_ZIndependent();
2696  DSourceComboUse locAllBut1ComboUse{Unknown, locAllBut1ZBin, locAllBut1ComboInfo, false, Unknown}; //Unknown -> particles
2697 
2698  //Get combos so far //not mixed charge: with-now is nullptr
2699  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, dComboInfoChargeContent[locAllBut1ComboInfo], nullptr);
2700 
2701  //if we can directly build the to-add decay use, then we will (only 1 combo requested)
2702  //else we must build a new use first, and then promote the to-add use when comboing horizontally
2703  auto locToAddComboInfo = (locFurtherDecays[0].second == 1) ? std::get<2>(locFurtherDecays[0].first) : GetOrMake_SourceComboInfo({}, {std::make_pair(locFurtherDecays[0].first, locFurtherDecays[0].second)}, locNumTabs);
2704  auto locToAddChargeContent = dComboInfoChargeContent[locToAddComboInfo];
2705  auto locToAddZBin = (locToAddChargeContent != d_Charged) ? locVertexZBin : DSourceComboInfo::Get_VertexZIndex_ZIndependent();
2706  auto locToAddComboUse = (locFurtherDecays[0].second == 1) ? locFurtherDecays[0].first : DSourceComboUse{Unknown, locToAddZBin, locToAddComboInfo, false, Unknown};
2707 
2708  // Now, see whether the combos for this grouping have already been done
2709  if(locSourceCombosByUseSoFar.find(locAllBut1ComboUse) == locSourceCombosByUseSoFar.end()) //if true: not yet
2710  {
2711  locComboUse_SubsetToBuild = locAllBut1ComboUse; //will build below before the end of the function
2712  locComboUse_SubsetToAdd = locToAddComboUse;
2713  locChargedCombo_SubsetToBuildPresiding = nullptr; //not needed when building particles
2714  }
2715  else //yes, it's already been done!
2716  {
2717  //just combo the All-but-1 combos to those from this decay and save the results
2718  Combo_Horizontally_AddCombo(locComboUseToCreate, locAllBut1ComboUse, locToAddComboUse, locComboingStage, locChargedCombo_Presiding, false, locNumTabs);
2719  return;
2720  }
2721  }
2722  }
2723  else //at least 2 of one type (decays / particles) needed:
2724  {
2725  //for each further decay map entry (e.g. pi0, 3), this is a collection of the uses representing those groupings //e.g. Unknown -> 3pi0
2726  //decays are sorted by: mixed-charge first, then fully-neutral, then fully-charged
2727  //within a charge: loop from heaviest-mass to least (most likely to be missing)
2728  for(auto locDecayIterator = locFurtherDecays.begin(); locDecayIterator != locFurtherDecays.end(); ++locDecayIterator)
2729  {
2730  if(locFurtherDecays.size() == 1)
2731  break; //will work with the particles instead
2732  if(dDebugLevel >= 5)
2733  cout << "2+ decays: try to build subsets with a decay missing" << endl;
2734 
2735  //build a DSourceComboUse with everything EXCEPT this set of decays, and see if it already exists
2736  //build the further-decays, removing this decay
2737  auto locFurtherDecaysToSearchFor = locFurtherDecays;
2738  const auto& locSourceComboUse_ThisDecay = locDecayIterator->first;
2739  locFurtherDecaysToSearchFor.erase(locFurtherDecaysToSearchFor.begin() + std::distance(locFurtherDecays.begin(), locDecayIterator));
2740 
2741  //if we can directly build the further-decays-to-search-for use, then we will (e.g. only one use & only 1 combo requested)
2742  //else we must build a new use first, and then expand the all-but-1 when comboing horizontally
2743  auto locAllBut1ComboInfo = std::get<2>((*locFurtherDecaysToSearchFor.begin()).first); //changed below if needed
2744 
2745  //if we can directly build the to-add decay use, then we will (only 1 combo requested)
2746  //else we must build a new use first, and then promote the to-add use when comboing horizontally
2747  auto locToAddComboInfo = (locDecayIterator->second == 1) ? std::get<2>(locSourceComboUse_ThisDecay) : GetOrMake_SourceComboInfo({}, {std::make_pair(locSourceComboUse_ThisDecay, locDecayIterator->second)}, locNumTabs);
2748  auto locToAddChargeContent = dComboInfoChargeContent[locToAddComboInfo];
2749  auto locToAddZBin = (locToAddChargeContent != d_Charged) ? locVertexZBin : DSourceComboInfo::Get_VertexZIndex_ZIndependent();
2750  auto locToAddComboUse = (locDecayIterator->second == 1) ? locSourceComboUse_ThisDecay : DSourceComboUse{Unknown, locToAddZBin, locToAddComboInfo, false, Unknown};
2751 
2752  //guard against special cases for the all-but-1 combo use //must be after check on whether all-but-1 is charged (it itself is special case)
2753  auto locAllBut1ComboUse = DSourceComboUse{Unknown, locVertexZBin, locAllBut1ComboInfo, false, Unknown}; //may change below
2754  auto locChargedCombo_PresidingToUse = locChargedCombo_Presiding; //may change below
2755  if((locFurtherDecaysToSearchFor.size() > 1) || !locNumParticlesNeeded.empty())
2756  {
2757  locAllBut1ComboInfo = GetOrMake_SourceComboInfo(locNumParticlesNeeded, locFurtherDecaysToSearchFor, locNumTabs);
2758  auto locAllBut1ZBin = (Get_ChargeContent(locAllBut1ComboInfo) != d_Charged) ? locVertexZBin : DSourceComboInfo::Get_VertexZIndex_ZIndependent();
2759  locAllBut1ComboUse = DSourceComboUse{Unknown, locAllBut1ZBin, locAllBut1ComboInfo, false, Unknown};
2760  }
2761  else if((locFurtherDecaysToSearchFor.size() == 1) && (locFurtherDecaysToSearchFor[0].second > 1))
2762  {
2763  locAllBut1ComboInfo = GetOrMake_SourceComboInfo({}, {std::make_pair(locSourceComboUse_ThisDecay, locDecayIterator->second)}, locNumTabs);
2764  auto locAllBut1ZBin = (Get_ChargeContent(locAllBut1ComboInfo) != d_Charged) ? locVertexZBin : DSourceComboInfo::Get_VertexZIndex_ZIndependent();
2765  locAllBut1ComboUse = DSourceComboUse{Unknown, locAllBut1ZBin, locAllBut1ComboInfo, false, Unknown};
2766  }
2767  else if(locFurtherDecaysToSearchFor.size() == 1)
2768  {
2769  locAllBut1ComboUse = (*locFurtherDecaysToSearchFor.begin()).first;
2770  locChargedCombo_PresidingToUse = Get_NextChargedCombo(locChargedCombo_Presiding, locAllBut1ComboUse, locComboingStage, true, 1);
2771  }
2772 
2773  if(dDebugLevel >= 20)
2774  {
2775  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2776  cout << "Comboing decays together. All-But-1 Use:" << endl;
2777  DAnalysis::Print_SourceComboUse(locAllBut1ComboUse, locNumTabs);
2778  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
2779  cout << "Comboing decays together. To-Add Use:" << endl;
2780  DAnalysis::Print_SourceComboUse(locToAddComboUse, locNumTabs);
2781  }
2782 
2783  //if on charged stage and in this loop, at least one decay must be charged
2784  //(if fully neutral wouldn't be trying to build, and if comboing with neutral the charged are represented by the decay e.g. X -> pi+, pi-, p ... etc.)
2785  //we will guard against the case of to-add == neutral in Combo_Horizontally_AddDecay()
2786  //here, we guard against the case of all-but-1 == neutral
2787  auto locAllBut1ChargeContent = dComboInfoChargeContent[locAllBut1ComboInfo];
2788  if((locComboingStage == d_ChargedStage) && (locAllBut1ChargeContent == d_Neutral))
2789  {
2790  //just copy the to-add (which is either charged or mixed) as the desired combos
2791  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locToAddChargeContent);
2792  auto& locToAddCombos = locSourceCombosByUseSoFar[locToAddComboUse];
2793  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, locToAddCombos);
2794  if(dDebugLevel > 0)
2795  cout << "Save for later!" << endl;
2796 
2797  //Set the resume indices
2798  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
2799  return;
2800  }
2801 
2802  if((locComboingStage != d_ChargedStage) && (locAllBut1ChargeContent == d_Charged))
2803  {
2804  //yes, it's already been done!
2805  //just combo the All-but-1 combos to those from this decay and return the results
2806  //don't promote particles or expand all-but-1: create new combo ABOVE all-but-1, that will contain all-but-1 and to-add side-by-side
2807  Combo_Horizontally_AddDecay(locComboUseToCreate, locAllBut1ComboUse, locToAddComboUse, locComboingStage, locChargedCombo_Presiding, false, locNumTabs);
2808  return;
2809  }
2810 
2811  //Get combos so far
2812  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locAllBut1ChargeContent, locChargedCombo_PresidingToUse);
2813 
2814  // Now, see whether the combos for this grouping have already been done
2815  if(locSourceCombosByUseSoFar.find(locAllBut1ComboUse) == locSourceCombosByUseSoFar.end()) //if true: not yet
2816  {
2817  //if on the first one (heaviest mass), save this subset in case we need to create it (if nothing else already done)
2818  if(locDecayIterator == locFurtherDecays.begin())
2819  {
2820  locComboUse_SubsetToBuild = locAllBut1ComboUse;
2821  locChargedCombo_SubsetToBuildPresiding = locChargedCombo_PresidingToUse;
2822  locComboUse_SubsetToAdd = locToAddComboUse;
2823  }
2824  continue; // try the next decay
2825  }
2826 
2827  //yes, it's already been done!
2828  //just combo the All-but-1 combos to those from this decay and save the results
2829  auto locExpandAllBut1Flag = Get_ExpandAllBut1Flag(locComboingStage, locAllBut1ComboUse, Get_ChargeContent(std::get<2>(locToAddComboUse)));
2830  Combo_Horizontally_AddDecay(locComboUseToCreate, locAllBut1ComboUse, locToAddComboUse, locComboingStage, locChargedCombo_Presiding, locExpandAllBut1Flag, locNumTabs);
2831  return;
2832  }
2833 
2834  //ok, none of the subsets without a decay has yet been created. let's try subsets without detected particles
2835  if((locComboingStage == d_ChargedStage) || (dComboInfoChargeContent[locComboInfoToCreate] == d_Neutral)) //no loose particles when mixing charged & neutral
2836  {
2837  if(dDebugLevel >= 5)
2838  cout << "try to build subsets with some detected particles missing" << endl;
2839  for(auto locParticleIterator = locNumParticlesNeeded.begin(); locParticleIterator != locNumParticlesNeeded.end(); ++locParticleIterator)
2840  {
2841  //build a DSourceComboUse with everything EXCEPT this set of particles, and see if it already exists
2842  //combo the particle horizontally, removing this PID
2843  auto locNumParticlesToSearchFor = locNumParticlesNeeded;
2844  const auto& locParticlePair = *locParticleIterator;
2845  locNumParticlesToSearchFor.erase(locNumParticlesToSearchFor.begin() + std::distance(locNumParticlesNeeded.begin(), locParticleIterator));
2846  if(dDebugLevel >= 5)
2847  cout << "particle pair: " << locParticlePair.first << ", " << int(locParticlePair.second) << endl;
2848 
2849  //build the DSourceComboUse
2850  auto locAllBut1ComboInfo = GetOrMake_SourceComboInfo(locNumParticlesToSearchFor, locFurtherDecays, locNumTabs);
2851  if((locComboingStage == d_ChargedStage) && (dComboInfoChargeContent[locAllBut1ComboInfo] == d_Neutral))
2852  continue; //this won't be done yet!
2853  auto locChargedContentAllBut1 = Get_ChargeContent(locAllBut1ComboInfo);
2854  auto locAllBut1ZBin = (locChargedContentAllBut1 != d_Charged) ? locVertexZBin : DSourceComboInfo::Get_VertexZIndex_ZIndependent();
2855  DSourceComboUse locAllBut1ComboUse(Unknown, locAllBut1ZBin, locAllBut1ComboInfo, false, Unknown); // Unknown -> everything but these particles
2856 
2857  //Get combos so far
2858  auto locChargedCombo_PresidingToUse = locChargedCombo_Presiding;
2859  if(locFurtherDecays.empty() && (locNumParticlesToSearchFor.size() == 1))
2860  locChargedCombo_PresidingToUse = Get_NextChargedCombo(locChargedCombo_Presiding, locAllBut1ComboUse, locComboingStage, true, 1);
2861  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locChargedContentAllBut1, locChargedCombo_PresidingToUse); //if not neutral then is on charged stage: argument doesn't matter
2862 
2863  // Now, see whether the combos for this grouping have already been done
2864  if(locSourceCombosByUseSoFar.find(locAllBut1ComboUse) == locSourceCombosByUseSoFar.end()) //if true: not yet
2865  {
2866  //if on the first one and there's no decays, save this subset in case we need to create it (if nothing else already done)
2867  if((locParticleIterator == locNumParticlesNeeded.begin()) && (locFurtherDecays.size() < 2))
2868  {
2869  locComboUse_SubsetToBuild = locAllBut1ComboUse;
2870  locChargedCombo_SubsetToBuildPresiding = locChargedCombo_PresidingToUse;
2871  locMissingSubsetIsDecayFlag = false;
2872  }
2873  continue; // try the next PID
2874  }
2875 
2876  //yes, it's already been done!
2877  //just combo the All-but-1 combos to those from this particle and return the results
2878  bool locExpandAllBut1Flag = false; //changed if following conditions hold
2879  if(std::get<0>(locAllBut1ComboUse) == Unknown) //check other conditions
2880  locExpandAllBut1Flag = (locAllBut1ComboInfo->Get_NumParticles().size() + locAllBut1ComboInfo->Get_FurtherDecays().size()) > 1; //true: has already been comboed horizontally once
2881  Combo_Horizontally_AddParticles(locComboUseToCreate, locAllBut1ComboUse, locParticlePair, locComboingStage, locChargedCombo_Presiding, locExpandAllBut1Flag, locNumTabs);
2882  return;
2883  }
2884  }
2885  }
2886 
2887  //none of the possible immediate subsets have been created
2888  //therefore, create one of them (the one without the heaviest particle), and then do the remaining combo
2889  if(dDebugLevel >= 5)
2890  cout << "build subset: dive down" << endl;
2891  Create_SourceCombos(locComboUse_SubsetToBuild, locComboingStage, locChargedCombo_SubsetToBuildPresiding, locNumTabs + 1); //this may be something we want to combo vertically: go back to (unknown) mother function
2892 
2893  //do the final combo!
2894  if(locMissingSubsetIsDecayFlag) //subset was missing a decay PID
2895  {
2896  if(dDebugLevel >= 5)
2897  cout << "do final add: add decay" << endl;
2898  auto locExpandAllBut1Flag = Get_ExpandAllBut1Flag(locComboingStage, locComboUse_SubsetToBuild, Get_ChargeContent(std::get<2>(locComboUse_SubsetToAdd)));
2899  Combo_Horizontally_AddDecay(locComboUseToCreate, locComboUse_SubsetToBuild, locComboUse_SubsetToAdd, locComboingStage, locChargedCombo_Presiding, locExpandAllBut1Flag, locNumTabs);
2900  }
2901  else //subset was missing a detected PID
2902  {
2903  if(dDebugLevel >= 5)
2904  cout << "do final add: add particles" << endl;
2905  auto locToAddChargeContent = (ParticleCharge(locNumParticlesNeeded.front().first) == 0) ? d_Neutral : d_Charged;
2906  auto locExpandAllBut1Flag = Get_ExpandAllBut1Flag(locComboingStage, locComboUse_SubsetToBuild, locToAddChargeContent);
2907  Combo_Horizontally_AddParticles(locComboUseToCreate, locComboUse_SubsetToBuild, locNumParticlesNeeded.front(), locComboingStage, locChargedCombo_Presiding, locExpandAllBut1Flag, locNumTabs);
2908  }
2909 }
2910 
2911 bool DSourceComboer::Get_ExpandAllBut1Flag(ComboingStage_t locComboingStage, const DSourceComboUse& locAllBut1ComboUse, Charge_t locToAddChargeContent)
2912 {
2913  if(std::get<0>(locAllBut1ComboUse) != Unknown)
2914  return false;
2915 
2916  auto locAllBut1ComboInfo = std::get<2>(locAllBut1ComboUse);
2917  auto locAllBut1ChargeContent = Get_ChargeContent(locAllBut1ComboInfo);
2918 
2919  //special case: if one is charged and the other is not: do not expand: must be side-by-side
2920  if((locAllBut1ChargeContent == d_Charged) != (locToAddChargeContent == d_Charged))
2921  return false;
2922  //don't expand if all-but-1 is mixed and merely contains a promoted charged combo
2923  if((locComboingStage == d_ChargedStage) && (locAllBut1ChargeContent == d_AllCharges))
2924  {
2925  size_t locNumNeutralUses = locAllBut1ComboInfo->Get_NumParticles().size();
2926  size_t locNumNonNeutralUses = 0;
2927  for(auto& locAllBut1DecayPair : locAllBut1ComboInfo->Get_FurtherDecays())
2928  {
2929  if(dComboInfoChargeContent[std::get<2>(locAllBut1DecayPair.first)] == d_Neutral)
2930  ++locNumNeutralUses;
2931  else
2932  ++locNumNonNeutralUses;
2933  }
2934  if((locNumNeutralUses >= 1) && (locNumNonNeutralUses == 1))
2935  return false; //merely a promoted charged combo
2936  }
2937 
2938  auto locExpandAllBut1Flag = ((locAllBut1ComboInfo->Get_NumParticles().size() + locAllBut1ComboInfo->Get_FurtherDecays().size()) > 1); //if true: has already been comboed horizontally once
2939  //special case: if only content is a single decay use, but > 1 combo of that use
2940  if(!locExpandAllBut1Flag && locAllBut1ComboInfo->Get_NumParticles().empty() && (locAllBut1ComboInfo->Get_FurtherDecays()[0].second > 1)) //
2941  locExpandAllBut1Flag = true;
2942  return locExpandAllBut1Flag;
2943 }
2944 
2945 void DSourceComboer::Combo_Horizontally_AddDecay(const DSourceComboUse& locComboUseToCreate, const DSourceComboUse& locComboUseAllBut1, const DSourceComboUse& locComboUseToAdd, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, bool locExpandAllBut1Flag, unsigned char locNumTabs)
2946 {
2947  auto locChargeContentUseToAdd = Get_ChargeContent(std::get<2>(locComboUseToAdd));
2948  if((locComboingStage == d_ChargedStage) && (locChargeContentUseToAdd == d_Neutral))
2949  {
2950  //this won't be done yet! just copy the all-but-1 as the desired combos
2951  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, d_Charged);
2952  auto& locAllBut1Combos = locSourceCombosByUseSoFar[locComboUseAllBut1];
2953  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, locAllBut1Combos);
2954  if(dDebugLevel > 0)
2955  cout << "Save for later!" << endl;
2956 
2957  //Set the resume indices
2958  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
2959  return;
2960  }
2961 
2962  //create the combos for the use-to-add if they haven't been created yet
2963  auto locChargedCombo_NextPresiding = Get_NextChargedCombo(locChargedCombo_Presiding, locComboUseToAdd, locComboingStage, true, 1);
2964  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locChargeContentUseToAdd, locChargedCombo_NextPresiding);
2965  if(locSourceCombosByUseSoFar.find(locComboUseToAdd) == locSourceCombosByUseSoFar.end()) //if true: not yet
2966  Create_SourceCombos(locComboUseToAdd, locComboingStage, locChargedCombo_Presiding, locNumTabs + 1); //this may be something we want to combo vertically: go back to (unknown) mother function
2967 
2968  //combo it horizontally with the rest
2969  Combo_Horizontally_AddCombo(locComboUseToCreate, locComboUseAllBut1, locComboUseToAdd, locComboingStage, locChargedCombo_Presiding, locExpandAllBut1Flag, locNumTabs);
2970 }
2971 
2972 void DSourceComboer::Combo_Horizontally_AddParticles(const DSourceComboUse& locComboUseToCreate, const DSourceComboUse& locComboUseAllBut1, const pair<Particle_t, unsigned char>& locParticlePairToAdd, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, bool locExpandAllBut1Flag, unsigned char locNumTabs)
2973 {
2974  if((locComboingStage == d_ChargedStage) && (ParticleCharge(locParticlePairToAdd.first) == 0))
2975  {
2976  //this won't be done yet! just copy the all-but-1 as the desired combos
2977  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, d_Charged);
2978  auto& locAllBut1Combos = locSourceCombosByUseSoFar[locComboUseAllBut1];
2979  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, locAllBut1Combos);
2980  if(dDebugLevel > 0)
2981  cout << "Save for later!" << endl;
2982 
2983  //Set the resume indices
2984  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
2985  return;
2986  }
2987  if(locParticlePairToAdd.second > 1)
2988  {
2989  //create a combo use for X -> N particles of this type
2990  auto locSourceInfoToAdd = GetOrMake_SourceComboInfo({locParticlePairToAdd}, {}, locNumTabs);
2991  auto locChargeContentUseToAdd = Get_ChargeContent(locSourceInfoToAdd);
2992  auto locToAddZBin = (locChargeContentUseToAdd != d_Charged) ? std::get<1>(locComboUseToCreate) : DSourceComboInfo::Get_VertexZIndex_ZIndependent();
2993  DSourceComboUse locComboUseToAdd(Unknown, locToAddZBin, locSourceInfoToAdd, false, Unknown);
2994 
2995  //create the combos for the use-to-add if they haven't been created yet
2996  auto locChargedCombo_NextPresiding = Get_NextChargedCombo(locChargedCombo_Presiding, locComboUseToAdd, locComboingStage, true, 1);
2997  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locChargeContentUseToAdd, locChargedCombo_NextPresiding);
2998  if(locSourceCombosByUseSoFar.find(locComboUseToAdd) == locSourceCombosByUseSoFar.end()) //if true: not yet
2999  Create_SourceCombos(locComboUseToAdd, locComboingStage, locChargedCombo_Presiding, locNumTabs + 1); //this may be something we want to combo vertically: go back to (unknown) mother function
3000 
3001  //combo it horizontally with the rest
3002  Combo_Horizontally_AddCombo(locComboUseToCreate, locComboUseAllBut1, locComboUseToAdd, locComboingStage, locChargedCombo_Presiding, locExpandAllBut1Flag, locNumTabs);
3003  }
3004  else
3005  Combo_Horizontally_AddParticle(locComboUseToCreate, locComboUseAllBut1, locParticlePairToAdd.first, locComboingStage, locChargedCombo_Presiding, locNumTabs);
3006 }
3007 
3008 void DSourceComboer::Create_Combo_OneParticle(const DSourceComboUse& locComboUseToCreate, ComboingStage_t locComboingStage, unsigned char locNumTabs)
3009 {
3010  if(dDebugLevel >= 5)
3011  {
3012  cout << endl;
3013  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3014  cout << "Create_Combo_OneParticle: Stage = " << locComboingStage << endl;
3015  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3016  cout << "USE TO CREATE:" << endl;
3017  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
3018  }
3019 
3020  //not much of a combo if there's only 1, is it? //e.g. 1 charged track at a vertex
3021 
3022  //Get combos so far
3023  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, d_Neutral); //if not neutral then is on charged stage: argument doesn't matter
3024 
3025  //get combo use contents
3026  auto locVertexZBin = std::get<1>(locComboUseToCreate);
3027  auto locComboInfo = std::get<2>(locComboUseToCreate);
3028  auto locParticlePair = locComboInfo->Get_NumParticles().front();
3029  auto locPID = locParticlePair.first;
3030 
3031  //check if comboing massive neutrals: comboing results are independent of z-bin: see if this has already been done for a different (non-z-independent) zbin
3032  if((ParticleCharge(locPID) == 0) && (ParticleMass(locPID) > 0.0) && (locComboingStage == d_MixedStage))
3033  {
3034  //if so, just copy the results into this zbin!
3035  //note that this is only the case for comboing particles, not comboing combos: decays contain uses, which have the zbin specified
3036  for(signed char locZBin = 0; locZBin < (signed char)dSourceComboTimeHandler->Get_NumVertexZBins(); ++locZBin)
3037  {
3038  if(locZBin == locVertexZBin)
3039  continue;
3040  auto locZBinUse = DSourceComboUse{Unknown, locZBin, locComboInfo, false, Unknown};
3041  if(locSourceCombosByUseSoFar.find(locZBinUse) == locSourceCombosByUseSoFar.end())
3042  continue;
3043 
3044  //just copy the results!
3045  locSourceCombosByUseSoFar[locComboUseToCreate] = locSourceCombosByUseSoFar[locZBinUse];
3046  auto& locSourceCombosByBeamBunchByUse = Get_SourceCombosByBeamBunchByUse(dComboInfoChargeContent[locComboInfo], nullptr);
3047  locSourceCombosByBeamBunchByUse.emplace(locComboUseToCreate, locSourceCombosByBeamBunchByUse[locZBinUse]);
3048 
3049  //Set the resume indices
3050  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
3051  if(dDebugLevel >= 5)
3052  {
3053  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3054  cout << "Create_Combo_OneParticle: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseSoFar[locComboUseToCreate]->size() << endl;
3055  }
3056  return;
3057  }
3058  }
3059 
3060  //if on the mixed stage, must be doing all neutrals: first copy over ALL fcal-only results
3061  locSourceCombosByUseSoFar.emplace(locComboUseToCreate, Get_SourceComboVectorResource());
3062  if(locComboingStage == d_MixedStage)
3063  Copy_ZIndependentMixedResults(locComboUseToCreate, nullptr);
3064 
3065  //Get particles for comboing
3066  const auto& locParticles = Get_ParticlesForComboing(locPID, locComboingStage, {}, locVertexZBin);
3067  for(const auto& locParticle : locParticles)
3068  {
3069  auto locIsZIndependent = Get_IsComboingZIndependent(locParticle, locPID);
3070  if((locComboingStage == d_MixedStage) && locIsZIndependent)
3071  continue; //this combo has already been created (assuming it was valid): during the FCAL-only stage
3072 
3073  //check if this combo is unique!
3074  //it will not be unique if: comboing photons, on mixed stage, and this combo has already been created for a different zbin
3075  //if so, don't duplicate memory
3076  //note that this can only the case for comboing particles, not comboing combos: decays contain uses, which have the zbin specified
3077  if((locPID == Gamma) && (locComboingStage == d_MixedStage))
3078  {
3079  auto locIterator = dNPhotonsToComboMap.find({locParticle});
3080  if(locIterator != dNPhotonsToComboMap.end()) //if true, already exists, save it and continue;
3081  {
3082  auto locCombo = locIterator->second;
3083  locSourceCombosByUseSoFar[locComboUseToCreate]->push_back(locCombo);
3084  if(dDebugLevel >= 15)
3085  {
3086  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3087  cout << "COPIED COMBO:" << endl;
3088  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
3089  }
3090  continue;
3091  }
3092  }
3093 
3094  auto locCombo = Get_SourceComboResource();
3095  locCombo->Set_Members({std::make_pair(locPID, locParticle)}, {}, locIsZIndependent);
3096  if(locPID == Gamma)
3097  dNPhotonsToComboMap.emplace(vector<const JObject*>{locParticle}, locCombo);
3098  if(dDebugLevel >= 10)
3099  {
3100  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3101  cout << "CREATED COMBO:" << endl;
3102  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
3103  }
3104 
3105  locSourceCombosByUseSoFar[locComboUseToCreate]->push_back(locCombo); //save it //in creation order
3106  if(locPID == Gamma)
3107  Register_ValidRFBunches(locComboUseToCreate, locCombo, dSourceComboTimeHandler->Get_ValidRFBunches(locParticle, locVertexZBin), locComboingStage, nullptr);
3108  else
3109  Register_ValidRFBunches(locComboUseToCreate, locCombo, {}, locComboingStage, nullptr);
3110  }
3111  if((dDebugLevel > 0) || (dDebugLevel == -1))
3112  Check_ForDuplicates(*(locSourceCombosByUseSoFar[locComboUseToCreate]));
3113 
3114  //Set the resume indices
3115  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
3116  if(dDebugLevel >= 5)
3117  {
3118  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3119  cout << "Create_Combo_OneParticle: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseSoFar[locComboUseToCreate]->size() << endl;
3120  }
3121 }
3122 
3123 void DSourceComboer::Create_Combo_OneDecay(const DSourceComboUse& locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, unsigned char locNumTabs)
3124 {
3125  //here you're literally just promoting something you created earlier. e.g. creating a X -> Lambda combo
3126  //this can happen if your channel is something like g, p -> (K+), Lambda
3127  if(dDebugLevel >= 5)
3128  {
3129  cout << endl;
3130  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3131  cout << "Create_Combo_OneDecay: Stage = " << locComboingStage << endl;
3132  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3133  cout << "USE TO CREATE:" << endl;
3134  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
3135  }
3136 
3137  auto locComboInfoToCreate = std::get<2>(locComboUseToCreate);
3138 
3139  auto locFurtherDecays = locComboInfoToCreate->Get_FurtherDecays();
3140  auto& locDecayUse = locFurtherDecays[0].first;
3141  auto locChargedCombo_PreviousPresiding = Get_NextChargedCombo(locChargedCombo_Presiding, locDecayUse, locComboingStage, true, 1);
3142 
3143  //if on the all-showers stage, first copy over ALL fcal-only results
3144  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, dComboInfoChargeContent[locComboInfoToCreate], locChargedCombo_PreviousPresiding);
3145  auto& locSourceCombosByUseToSaveTo = Get_CombosSoFar(locComboingStage, dComboInfoChargeContent[locComboInfoToCreate], locChargedCombo_Presiding);
3146  locSourceCombosByUseToSaveTo.emplace(locComboUseToCreate, Get_SourceComboVectorResource());
3147 
3148  if(locComboingStage == d_MixedStage)
3149  Copy_ZIndependentMixedResults(locComboUseToCreate, locChargedCombo_Presiding);
3150 
3151  //The decay itself has already been created (when comboing vertically)
3152  auto& locDecayCombos = locSourceCombosByUseSoFar[locDecayUse];
3153  for(auto locDecayCombo : *locDecayCombos)
3154  {
3155  auto locIsZIndependent = locDecayCombo->Get_IsComboingZIndependent();
3156  if((locComboingStage == d_MixedStage) && locIsZIndependent)
3157  continue; //this combo has already been created (assuming it was valid): during the FCAL-only stage
3158 
3159  auto& locValidRFBunches = dValidRFBunches_ByCombo[std::make_pair(locDecayCombo, std::get<1>(locDecayUse))];
3160 
3161  //create new combo!
3162  auto locCombo = Get_SourceComboResource();
3163  locCombo->Set_Members({}, {std::make_pair(locDecayUse, vector<const DSourceCombo*>{locDecayCombo})}, locIsZIndependent);
3164 
3165  //save it!
3166  locSourceCombosByUseToSaveTo[locComboUseToCreate]->push_back(locCombo);
3167  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, locChargedCombo_Presiding);
3168  }
3169 
3170  //Set the resume indices
3171  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
3172  if(dDebugLevel >= 5)
3173  {
3174  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3175  cout << "Create_Combo_OneDecay: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseToSaveTo[locComboUseToCreate]->size() << endl;
3176  }
3177 }
3178 
3179 void DSourceComboer::Combo_Horizontally_AddCombo(const DSourceComboUse& locComboUseToCreate, const DSourceComboUse& locAllBut1ComboUse, const DSourceComboUse& locSourceComboUseToAdd, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, bool locExpandAllBut1Flag, unsigned char locNumTabs)
3180 {
3181  if(dDebugLevel >= 5)
3182  {
3183  cout << endl;
3184  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3185  cout << "Combo_Horizontally_AddCombo: Stage, presiding charged combo = " << locComboingStage << ", " << locChargedCombo_Presiding << endl;
3186  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3187  cout << "PRESIDING COMBO:" << endl;
3188  DAnalysis::Print_SourceCombo(locChargedCombo_Presiding, locNumTabs);
3189  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3190  cout << "USE TO CREATE:" << endl;
3191  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
3192  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3193  cout << "All-But-1 USE:" << endl;
3194  DAnalysis::Print_SourceComboUse(locAllBut1ComboUse, locNumTabs);
3195  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3196  cout << "USE TO ADD:" << endl;
3197  DAnalysis::Print_SourceComboUse(locSourceComboUseToAdd, locNumTabs);
3198  }
3199 
3200  //e.g. we are grouping N pi0s and M photons (> 1) with L etas (>= 1), etc. to make combos
3201  //so, let's get the combos for the main grouping
3202 
3203  //Get combos so far
3204  auto locComboInfo_AllBut1 = std::get<2>(locAllBut1ComboUse);
3205  auto locComboInfoToCreate = std::get<2>(locComboUseToCreate);
3206  auto locChargeContent_AllBut1 = dComboInfoChargeContent[locComboInfo_AllBut1];
3207  auto locChargedCombo_WithNow = Get_ChargedCombo_WithNow(locChargedCombo_Presiding, locComboInfoToCreate, locComboingStage);
3208 //cout << "save to: stage, charge, withnow = " << locComboingStage << ", " << dComboInfoChargeContent[locComboInfoToCreate] << ", " << locChargedCombo_WithNow << endl;
3209  auto& locSourceCombosByUseToSaveTo = Get_CombosSoFar(locComboingStage, dComboInfoChargeContent[locComboInfoToCreate], locChargedCombo_Presiding);
3210 
3211  bool locGetFromSoFarFlag = (locComboingStage == d_ChargedStage) || (locChargeContent_AllBut1 != d_Charged);
3212  auto locChargedCombo_PresidingAllBut1 = locChargedCombo_Presiding;
3213  if((locComboInfo_AllBut1->Get_FurtherDecays().size() == 1) && locComboInfo_AllBut1->Get_NumParticles().empty())
3214  locChargedCombo_PresidingAllBut1 = Get_NextChargedCombo(locChargedCombo_Presiding, locAllBut1ComboUse, locComboingStage, true, 1);
3215  auto& locSourceCombosByUseAllBut1 = Get_CombosSoFar(locComboingStage, locChargeContent_AllBut1, locChargedCombo_PresidingAllBut1);
3216 
3217  vector<const DSourceCombo*> locChargedComboVector = {locChargedCombo_WithNow}; //ugh
3218  auto locCombos_AllBut1 = locGetFromSoFarFlag ? locSourceCombosByUseAllBut1[locAllBut1ComboUse] : &locChargedComboVector; //Combos are a vector of (e.g.): -> N pi0s
3219 
3220  auto locChargeContent_ToAdd = dComboInfoChargeContent[std::get<2>(locSourceComboUseToAdd)];
3221  if((locComboingStage == d_ChargedStage) && (locChargeContent_ToAdd == d_Neutral))
3222  {
3223  //can't add neutrals, so we are already done! just copy the results to the new vector
3224  locSourceCombosByUseToSaveTo.emplace(locComboUseToCreate, locCombos_AllBut1);
3225  //Set the resume indices
3226  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
3227  if(dDebugLevel > 0)
3228  cout << "Save for later!" << endl;
3229  return;
3230  }
3231 
3232  //if on the all-showers stage, first copy over ALL fcal-only results
3233  locSourceCombosByUseToSaveTo.emplace(locComboUseToCreate, Get_SourceComboVectorResource());
3234  if(locComboingStage == d_MixedStage)
3235  Copy_ZIndependentMixedResults(locComboUseToCreate, locChargedCombo_Presiding);
3236 
3237  if(locCombos_AllBut1->empty())
3238  return; //nothing to do
3239 
3240  auto locDecayPID_UseToAdd = std::get<0>(locSourceComboUseToAdd);
3241  auto locComboInfo_UseToAdd = std::get<2>(locSourceComboUseToAdd);
3242 
3243  //determine whether we should promote the contents of the combos we are combining up to the new combo (else set combo as decay of new combo)
3244  auto locComboInfo_UseToCreate = std::get<2>(locComboUseToCreate);
3245  DSourceComboUse locNonNeutralUse{Unknown, 0, nullptr, false, Unknown};
3246  bool locPromoteToAddFlag = Get_PromoteFlag(locComboingStage, locDecayPID_UseToAdd, locComboInfo_UseToCreate, locComboInfo_UseToAdd, locNonNeutralUse); //is ignored if charged
3247  bool locPromoteAllBut1Flag = Get_PromoteFlag(locComboingStage, std::get<0>(locAllBut1ComboUse), locComboInfo_UseToCreate, locComboInfo_AllBut1, locNonNeutralUse);
3248  if(dDebugLevel >= 20)
3249  cout << "flags: expand all-but-1, promote to-add, promote all-but-1: " << locExpandAllBut1Flag << ", " << locPromoteToAddFlag << ", " << locPromoteAllBut1Flag << endl;
3250 
3251  //check if on mixed stage but comboing to charged
3252  if((locComboingStage != d_ChargedStage) && (locChargeContent_ToAdd == d_Charged))
3253  {
3254  //only one valid option for to-add: locChargedCombo_WithNow: create all combos immediately
3255  for(const auto& locCombo_AllBut1 : *locCombos_AllBut1)
3256  {
3257  auto locIsZIndependent = locCombo_AllBut1->Get_IsComboingZIndependent();
3258  if((locComboingStage == d_MixedStage) && locIsZIndependent)
3259  continue; //this combo has already been created (assuming it was valid): during the FCAL-only stage
3260 
3261  //get the valid RF bunches (those for the all-but-1, because we are comboing with charged which is "all")
3262  const auto& locValidRFBunches = dValidRFBunches_ByCombo[std::make_pair(locCombo_AllBut1, std::get<1>(locAllBut1ComboUse))];
3263 
3264  //create new combo!
3265  auto locCombo = Get_SourceComboResource();
3266 
3267  //get contents of the all-but-1 so that we can add to them
3268  auto locFurtherDecayCombos_AllBut1 = locCombo_AllBut1->Get_FurtherDecayCombos(); //the all-but-1 combo contents by use
3269  auto locComboParticles_AllBut1 = locCombo_AllBut1->Get_SourceParticles();
3270 
3271  if(locExpandAllBut1Flag)
3272  {
3273  locFurtherDecayCombos_AllBut1.emplace_back(locSourceComboUseToAdd, vector<const DSourceCombo*>{locChargedCombo_WithNow});
3274  locCombo->Set_Members(locComboParticles_AllBut1, locFurtherDecayCombos_AllBut1, locIsZIndependent); // create combo with all PIDs
3275  }
3276  else
3277  {
3278  if(locPromoteAllBut1Flag)
3279  {
3280  //promote contents of all-but-1 above the to-add level
3281  //so, really, use the all-but-1 as the basis, and put the to-add as a another decay in the all-but-1
3282  locFurtherDecayCombos_AllBut1.emplace_back(locSourceComboUseToAdd, vector<const DSourceCombo*>{locChargedCombo_WithNow});
3283  locCombo->Set_Members(locComboParticles_AllBut1, locFurtherDecayCombos_AllBut1, locIsZIndependent);
3284  }
3285  else //no promotions: side by side in a new combo
3286  {
3287  DSourceCombosByUse_Small locFurtherDecayCombos_Needed;
3288  locFurtherDecayCombos_Needed.emplace_back(locAllBut1ComboUse, vector<const DSourceCombo*>{locCombo_AllBut1});
3289  locFurtherDecayCombos_Needed.emplace_back(locSourceComboUseToAdd, vector<const DSourceCombo*>{locChargedCombo_WithNow});
3290  locCombo->Set_Members({}, locFurtherDecayCombos_Needed, locIsZIndependent); // create combo with all PIDs
3291  }
3292  }
3293  if(dDebugLevel >= 10)
3294  {
3295  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3296  cout << "CREATED COMBO:" << endl;
3297  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
3298  }
3299 
3300  //save it!
3301  locSourceCombosByUseToSaveTo[locComboUseToCreate]->push_back(locCombo);
3302  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, locChargedCombo_Presiding);
3303  }
3304  if((dDebugLevel > 0) || (dDebugLevel == -1))
3305  Check_ForDuplicates(*(locSourceCombosByUseToSaveTo[locComboUseToCreate]));
3306 
3307  //Set the resume indices
3308  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
3309  if(dDebugLevel >= 5)
3310  {
3311  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3312  cout << "Combo_Horizontally_AddCombo: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseToSaveTo[locComboUseToCreate]->size() << endl;
3313  }
3314  return;
3315  }
3316 
3317  //get the previous charged combo (if needed)
3318  auto locNextPresidingCombo = Get_NextChargedCombo(locChargedCombo_Presiding, locSourceComboUseToAdd, locComboingStage, true, 1);
3319  auto locChargedCombo_WithPrevious = Get_ChargedCombo_WithNow(locNextPresidingCombo, std::get<2>(locSourceComboUseToAdd), locComboingStage);
3320  if(dDebugLevel >= 20)
3321  cout << "combos: presiding, next-presiding, with-previous: " << locChargedCombo_Presiding << ", " << locNextPresidingCombo << ", " << locChargedCombo_WithPrevious << endl;
3322 
3323  //now, for each combo of all-but-1-PIDs, see which of the to-add combos we can group to it
3324  //valid grouping: Don't re-use a shower we've already used
3325  for(const auto& locCombo_AllBut1 : *locCombos_AllBut1)
3326  {
3327  //first of all, get the potential combos that satisfy the RF bunches for the all-but-1 combo
3328  const auto& locValidRFBunches_AllBut1 = dValidRFBunches_ByCombo[std::make_pair(locCombo_AllBut1, std::get<1>(locAllBut1ComboUse))];
3329  const auto& locDecayCombos_ToAdd = Get_CombosForComboing(locSourceComboUseToAdd, locComboingStage, locValidRFBunches_AllBut1, locNextPresidingCombo);
3330 
3331  //before we loop, first get all of the showers used to make the all-but-1 grouping, and sort it so that we can quickly search it
3332  auto locUsedParticles_AllBut1 = DAnalysis::Get_SourceParticles(locCombo_AllBut1->Get_SourceParticles(true)); //true: entire chain
3333  std::sort(locUsedParticles_AllBut1.begin(), locUsedParticles_AllBut1.end()); //must sort, because when retrieving entire chain is unsorted
3334 
3335  //this function will do our validity test
3336  auto Search_Duplicates = [&locUsedParticles_AllBut1](const JObject* locParticle) -> bool
3337  {return std::binary_search(locUsedParticles_AllBut1.begin(), locUsedParticles_AllBut1.end(), locParticle);};
3338 
3339  auto locIsZIndependent_AllBut1 = locCombo_AllBut1->Get_IsComboingZIndependent();
3340 
3341  //loop over potential combos to add to the group, creating a new combo for each valid (non-duplicate) grouping
3342  for(const auto& locDecayCombo_ToAdd : locDecayCombos_ToAdd)
3343  {
3344  auto locIsZIndependent = (locIsZIndependent_AllBut1 && locDecayCombo_ToAdd->Get_IsComboingZIndependent());
3345  if((locComboingStage == d_MixedStage) && locIsZIndependent)
3346  continue; //this combo has already been created (assuming it was valid): during the FCAL-only stage
3347 
3348  //search the all-but-1 shower vector to see if any of the showers in this combo are duplicated
3349  auto locUsedParticles_ToAdd = DAnalysis::Get_SourceParticles(locDecayCombo_ToAdd->Get_SourceParticles(true)); //true: entire chain
3350 
3351  //conduct search
3352  if(std::any_of(locUsedParticles_ToAdd.begin(), locUsedParticles_ToAdd.end(), Search_Duplicates))
3353  continue; //at least one photon was a duplicate, this combo won't work
3354 
3355  //no duplicates: this combo is unique. build a new combo
3356 
3357  //See which RF bunches match up //guaranteed to be at least one, due to selection in Get_CombosForComboing() function
3358  vector<int> locValidRFBunches = {}; //if charged or massive neutrals, ignore (they don't choose at this stage)
3359  if(locComboingStage != d_ChargedStage)
3360  locValidRFBunches = dSourceComboTimeHandler->Get_CommonRFBunches(locValidRFBunches_AllBut1, dValidRFBunches_ByCombo[std::make_pair(locDecayCombo_ToAdd, std::get<1>(locSourceComboUseToAdd))]);
3361 
3362  //create new combo!
3363  auto locCombo = Get_SourceComboResource();
3364 
3365  //get contents of the all-but-1 so that we can add to them
3366  auto locFurtherDecayCombos_AllBut1 = locCombo_AllBut1->Get_FurtherDecayCombos(); //the all-but-1 combo contents by use
3367  auto locComboParticles_AllBut1 = locCombo_AllBut1->Get_SourceParticles(false);
3368  if(locExpandAllBut1Flag)
3369  {
3370  if(locPromoteToAddFlag)
3371  {
3372  //promote all contents of to-add to the all-but-1 level
3373  auto locUsedParticlePairs_ToAdd = locDecayCombo_ToAdd->Get_SourceParticles(false);
3374  locComboParticles_AllBut1.insert(locComboParticles_AllBut1.end(), locUsedParticlePairs_ToAdd.begin(), locUsedParticlePairs_ToAdd.end());
3375  auto locFurtherDecayCombos_ToAdd = locDecayCombo_ToAdd->Get_FurtherDecayCombos();
3376  locFurtherDecayCombos_AllBut1.insert(locFurtherDecayCombos_AllBut1.end(), locFurtherDecayCombos_ToAdd.begin(), locFurtherDecayCombos_ToAdd.end());
3377  }
3378  else
3379  locFurtherDecayCombos_AllBut1.emplace_back(locSourceComboUseToAdd, vector<const DSourceCombo*>{locDecayCombo_ToAdd});
3380  locCombo->Set_Members(locComboParticles_AllBut1, locFurtherDecayCombos_AllBut1, locIsZIndependent); // create combo with all PIDs
3381  }
3382  else //side by side in a new combo
3383  {
3384  auto locComboParticlePairs_ToAdd = locDecayCombo_ToAdd->Get_SourceParticles(false);
3385  auto locFurtherDecayCombos_ToAdd = locDecayCombo_ToAdd->Get_FurtherDecayCombos();
3386  if(locPromoteAllBut1Flag && locPromoteToAddFlag)
3387  {
3388  //union of particles & decays from each
3389  //so, use the all-but-1 as a basis, and merge the to-add content in at the same level
3390  locFurtherDecayCombos_AllBut1.insert(locFurtherDecayCombos_AllBut1.end(), locFurtherDecayCombos_ToAdd.begin(), locFurtherDecayCombos_ToAdd.end());
3391  locComboParticles_AllBut1.insert(locComboParticles_AllBut1.end(), locComboParticlePairs_ToAdd.begin(), locComboParticlePairs_ToAdd.end());
3392  locCombo->Set_Members(locComboParticles_AllBut1, locFurtherDecayCombos_AllBut1, locIsZIndependent);
3393  }
3394  else if(locPromoteAllBut1Flag)
3395  {
3396  //promote contents of all-but-1 above the to-add level
3397  //so, really, use the all-but-1 as the basis, and put the to-add as a another decay in the all-but-1
3398  locFurtherDecayCombos_AllBut1.emplace_back(locSourceComboUseToAdd, vector<const DSourceCombo*>{locDecayCombo_ToAdd});
3399  locCombo->Set_Members(locComboParticles_AllBut1, locFurtherDecayCombos_AllBut1, locIsZIndependent);
3400  }
3401  else if(locPromoteToAddFlag)
3402  {
3403  //promote contents of to-add above the all-but-1 level
3404  //so, really, use the to-add as the basis, and put the all-but-1 as a another decay in the to-add
3405  //if non-neutral-use info is not nullptr, then the all-but-1 is a charged combo that has been promoted: keep the combo, but change the use
3406  auto locAllBut1UseToUse = (std::get<2>(locNonNeutralUse) == nullptr) ? locAllBut1ComboUse : locNonNeutralUse;
3407  locFurtherDecayCombos_ToAdd.emplace_back(locAllBut1UseToUse, vector<const DSourceCombo*>{locCombo_AllBut1});
3408  locCombo->Set_Members(locComboParticlePairs_ToAdd, locFurtherDecayCombos_ToAdd, locIsZIndependent);
3409  }
3410  else //promote nothing
3411  {
3412  DSourceCombosByUse_Small locFurtherDecayCombos_Needed;
3413  //if non-neutral-use info is not nullptr, then the all-but-1 is a charged combo that has been promoted: keep the combo, but change the use
3414  auto locAllBut1UseToUse = (std::get<2>(locNonNeutralUse) == nullptr) ? locAllBut1ComboUse : locNonNeutralUse;
3415  locFurtherDecayCombos_Needed.emplace_back(locAllBut1UseToUse, vector<const DSourceCombo*>{locCombo_AllBut1});
3416  locFurtherDecayCombos_Needed.emplace_back(locSourceComboUseToAdd, vector<const DSourceCombo*>{locDecayCombo_ToAdd});
3417  locCombo->Set_Members({}, locFurtherDecayCombos_Needed, locIsZIndependent);
3418  }
3419  }
3420  if(dDebugLevel >= 10)
3421  {
3422  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3423  cout << "CREATED COMBO:" << endl;
3424  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
3425  }
3426 
3427  //save it! //in creation order!
3428  locSourceCombosByUseToSaveTo[locComboUseToCreate]->push_back(locCombo);
3429  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, locChargedCombo_Presiding);
3430  }
3431  }
3432  if((dDebugLevel > 0) || (dDebugLevel == -1))
3433  Check_ForDuplicates(*(locSourceCombosByUseToSaveTo[locComboUseToCreate]));
3434 
3435  //Set the resume indices
3436  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
3437  if(dDebugLevel >= 5)
3438  {
3439  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3440  cout << "Combo_Horizontally_AddCombo: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseToSaveTo[locComboUseToCreate]->size() << endl;
3441  }
3442 }
3443 
3444 void DSourceComboer::Combo_Horizontally_AddParticle(const DSourceComboUse& locComboUseToCreate, const DSourceComboUse& locAllBut1ComboUse, Particle_t locPID, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding, unsigned char locNumTabs)
3445 {
3446  if(dDebugLevel >= 5)
3447  {
3448  cout << endl;
3449  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3450  cout << "Combo_Horizontally_AddParticle: Stage, PID-to-add, presiding charged combo = " << locComboingStage << ", " << locPID << ", " << locChargedCombo_Presiding << endl;
3451  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3452  cout << "PRESIDING COMBO:" << endl;
3453  DAnalysis::Print_SourceCombo(locChargedCombo_Presiding, locNumTabs);
3454  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3455  cout << "USE TO CREATE:" << endl;
3456  DAnalysis::Print_SourceComboUse(locComboUseToCreate, locNumTabs);
3457  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3458  cout << "All-But-1 USE:" << endl;
3459  DAnalysis::Print_SourceComboUse(locAllBut1ComboUse, locNumTabs);
3460  }
3461 
3462  //Get combos so far
3463  auto locComboInfo_AllBut1 = std::get<2>(locAllBut1ComboUse);
3464  auto locChargeContent_AllBut1 = dComboInfoChargeContent[locComboInfo_AllBut1];
3465  auto locComboInfoToCreate = std::get<2>(locAllBut1ComboUse);
3466  auto locChargedCombo_WithNow = Get_ChargedCombo_WithNow(locChargedCombo_Presiding, locComboInfoToCreate, locComboingStage);
3467  auto locChargedCombo_PresidingAllBut1 = locChargedCombo_Presiding;
3468  if((locComboInfo_AllBut1->Get_FurtherDecays().size() == 1) && locComboInfo_AllBut1->Get_NumParticles().empty())
3469  locChargedCombo_PresidingAllBut1 = Get_NextChargedCombo(locChargedCombo_Presiding, locAllBut1ComboUse, locComboingStage, true, 1);
3470  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(locComboingStage, locChargeContent_AllBut1, locChargedCombo_PresidingAllBut1);
3471 
3472  //e.g. we are grouping a whole bunch of particles and decays with a lone particle to make new combos
3473  auto& locSourceCombosByUseToSaveTo = Get_CombosSoFar(locComboingStage, dComboInfoChargeContent[locComboInfoToCreate], locChargedCombo_Presiding);
3474 
3475  vector<const DSourceCombo*> locChargedComboVector = {locChargedCombo_WithNow}; //ugh
3476  bool locGetFromSoFarFlag = (locComboingStage == d_ChargedStage) || (locChargeContent_AllBut1 != d_Charged);
3477  auto locCombos_AllBut1 = locGetFromSoFarFlag ? locSourceCombosByUseSoFar[locAllBut1ComboUse] : &locChargedComboVector; //Combos are a vector of (e.g.): -> N pi0s
3478 
3479  if((locComboingStage == d_ChargedStage) && (ParticleCharge(locPID) == 0))
3480  {
3481  //can't add neutrals, so we are already done! just copy the results to the new vector
3482  locSourceCombosByUseToSaveTo[locComboUseToCreate] = locCombos_AllBut1;
3483  //Set the resume indices
3484  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, nullptr);
3485  if(dDebugLevel > 0)
3486  cout << "Save for later!" << endl;
3487  return;
3488  }
3489 
3490  //if on the all-showers stage, first copy over ALL fcal-only results
3491  locSourceCombosByUseToSaveTo.emplace(locComboUseToCreate, Get_SourceComboVectorResource());
3492  if(locComboingStage == d_MixedStage)
3493  Copy_ZIndependentMixedResults(locComboUseToCreate, nullptr);
3494 
3495  auto locVertexZBin = std::get<1>(locComboUseToCreate);
3496 
3497  //loop over the combos
3498  for(const auto& locCombo_AllBut1 : *locCombos_AllBut1)
3499  {
3500  //now, for each combo of all-but-1-PIDs, see which of the particles can group to it
3501  //valid grouping: Don't re-use a particle we've already used
3502 
3503  //before we loop, first get all of the particles of the given PID used to make the all-but-1 grouping, and sort it so that we can quickly search it
3504  auto locUsedParticlePairs_AllBut1 = locCombo_AllBut1->Get_SourceParticles(true);
3505 
3506  auto locUsedParticles_AllBut1 = DAnalysis::Get_SourceParticles(locUsedParticlePairs_AllBut1, ParticleCharge(locPID)); //true: entire chain
3507  std::sort(locUsedParticles_AllBut1.begin(), locUsedParticles_AllBut1.end()); //necessary: may be out of order due to comboing of different decays
3508 
3509  //also, pre-get the further decays & FCAL-only flag, as we'll need them to build new combos
3510  auto locFurtherDecays = locCombo_AllBut1->Get_FurtherDecayCombos(); //the all-but-1 combo contents by use
3511  auto locIsZIndependent_AllBut1 = locCombo_AllBut1->Get_IsComboingZIndependent();
3512 
3513  //Get potential particles for comboing
3514  const auto& locValidRFBunches_AllBut1 = dValidRFBunches_ByCombo[std::make_pair(locCombo_AllBut1, std::get<1>(locAllBut1ComboUse))];
3515  const auto& locParticles = Get_ParticlesForComboing(locPID, locComboingStage, locValidRFBunches_AllBut1, locVertexZBin);
3516 
3517  //loop over potential showers to add to the group, creating a new combo for each valid (non-duplicate) grouping
3518  for(const auto& locParticle : locParticles)
3519  {
3520  auto locIsZIndependent = (locComboingStage == d_MixedStage_ZIndependent) || (locIsZIndependent_AllBut1 && Get_IsComboingZIndependent(locParticle, locPID));
3521  if((locComboingStage == d_MixedStage) && locIsZIndependent)
3522  continue; //this combo has already been created (assuming it was valid): during the FCAL-only stage
3523 
3524  //conduct search
3525  if(std::binary_search(locUsedParticles_AllBut1.begin(), locUsedParticles_AllBut1.end(), locParticle))
3526  continue; //this particle has already been used, this combo won't work
3527 
3528  //See which RF bunches match up //guaranteed to be at least one, due to selection in Get_ParticlesForComboing() function
3529  //if charged or massive neutrals, ignore (they don't choose at this stage)
3530  vector<int> locValidRFBunches = (locPID != Gamma) ? locValidRFBunches_AllBut1 : dSourceComboTimeHandler->Get_CommonRFBunches(locValidRFBunches_AllBut1, locParticle, locVertexZBin);
3531 
3532  //no duplicates: this combo is unique. build a new combo
3533  auto locComboParticles = locCombo_AllBut1->Get_SourceParticles(false);
3534  locComboParticles.emplace_back(locPID, locParticle);
3535  auto locCombo = Get_SourceComboResource();
3536  locCombo->Set_Members(locComboParticles, locFurtherDecays, locIsZIndependent); // create combo with all PIDs
3537  if(dDebugLevel >= 10)
3538  {
3539  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3540  cout << "CREATED COMBO:" << endl;
3541  DAnalysis::Print_SourceCombo(locCombo, locNumTabs);
3542  }
3543 
3544  //save it! //in creation order!
3545  locSourceCombosByUseToSaveTo[locComboUseToCreate]->push_back(locCombo);
3546  Register_ValidRFBunches(locComboUseToCreate, locCombo, locValidRFBunches, locComboingStage, locChargedCombo_Presiding);
3547  }
3548  }
3549  if((dDebugLevel > 0) || (dDebugLevel == -1))
3550  Check_ForDuplicates(*(locSourceCombosByUseToSaveTo[locComboUseToCreate]));
3551 
3552  //Set the resume indices
3553  Build_ComboResumeIndices(locComboUseToCreate, locComboingStage, locChargedCombo_Presiding);
3554  if(dDebugLevel >= 5)
3555  {
3556  for(decltype(locNumTabs) locTabNum = 0; locTabNum < locNumTabs; ++locTabNum) cout << "\t";
3557  cout << "Combo_Horizontally_AddParticle: NUM SOURCE COMBOS CREATED: " << locSourceCombosByUseToSaveTo[locComboUseToCreate]->size() << endl;
3558  }
3559 }
3560 
3561 /***************************************************************** PARTICLE UTILITY FUNCTIONS *****************************************************************/
3562 
3563 const vector<const JObject*>& DSourceComboer::Get_ParticlesForComboing(Particle_t locPID, ComboingStage_t locComboingStage, const vector<int>& locBeamBunches, signed char locVertexZBin)
3564 {
3565  //find all particles that have an overlapping beam bunch with the input
3566 
3567  //SPECIAL CASES FOR NEUTRALS:
3568  //massive neutral: all showers
3569  //unknown RF: All showers
3570  //unknown vertex, known RF: from each zbin, all showers that were valid for that rf bunch (already setup)
3571 
3572  if(ParticleCharge(locPID) != 0) //charged tracks
3573  return dTracksByPID[locPID]; //rf bunch & vertex-z are irrelevant
3574  else if(locPID != Gamma) //massive neutrals
3575  return dShowersByBeamBunchByZBin[DSourceComboInfo::Get_VertexZIndex_Unknown()][{}]; //all neutrals: cannot do PID at all, and cannot do mass cuts until a specific vertex is chosen, so vertex-z doesn't matter
3576 
3577  if(locComboingStage == d_MixedStage_ZIndependent) //fcal
3578  {
3580  auto locGroupBunchIterator = dShowersByBeamBunchByZBin[locVertexZBin].find(locBeamBunches);
3581  if(locGroupBunchIterator != dShowersByBeamBunchByZBin[locVertexZBin].end())
3582  return locGroupBunchIterator->second;
3583  return Get_ShowersByBeamBunch(locBeamBunches, dShowersByBeamBunchByZBin[locVertexZBin], locVertexZBin);
3584  }
3585 
3586  if(locBeamBunches.empty())
3587  return dShowersByBeamBunchByZBin[DSourceComboInfo::Get_VertexZIndex_Unknown()][{}]; //all showers, regardless of vertex-z
3588 
3589  auto locGroupBunchIterator = dShowersByBeamBunchByZBin[locVertexZBin].find(locBeamBunches);
3590  if(locGroupBunchIterator != dShowersByBeamBunchByZBin[locVertexZBin].end())
3591  return locGroupBunchIterator->second;
3592  return Get_ShowersByBeamBunch(locBeamBunches, dShowersByBeamBunchByZBin[locVertexZBin], locVertexZBin);
3593 }
3594 
3595 const vector<const JObject*>& DSourceComboer::Get_ShowersByBeamBunch(const vector<int>& locBeamBunches, DPhotonShowersByBeamBunch& locShowersByBunch, signed char locVertexZBin)
3596 {
3597  if(locBeamBunches.empty())
3598  return locShowersByBunch[{}];
3599 
3600  //find all particles that have an overlapping beam bunch with the input
3601  //this won't happen often (max probably tens of times each event), so we can be a little inefficient
3602  vector<int> locBunchesSoFar = {*locBeamBunches.begin()};
3603  for(auto locBunchIterator = std::next(locBeamBunches.begin()); locBunchIterator != locBeamBunches.end(); ++locBunchIterator)
3604  {
3605  const auto& locComboShowers = locShowersByBunch[locBunchesSoFar];
3606  const auto& locBunchShowers = locShowersByBunch[{*locBunchIterator}];
3607 
3608  locBunchesSoFar.push_back(*locBunchIterator);
3609  if(locShowersByBunch.find(locBunchesSoFar) != locShowersByBunch.end())
3610  continue; //this subset already created and indexed
3611 
3612  if(locBunchShowers.empty())
3613  {
3614  locShowersByBunch.emplace(locBunchesSoFar, locComboShowers);
3615  Build_ParticleIndices(Gamma, locBeamBunches, locShowersByBunch[locBunchesSoFar], locVertexZBin);
3616  continue;
3617  }
3618 
3619  //merge and move-emplace
3620  vector<const JObject*> locMergeResult;
3621  locMergeResult.reserve(locComboShowers.size() + locBunchShowers.size());
3622  std::set_union(locComboShowers.begin(), locComboShowers.end(), locBunchShowers.begin(), locBunchShowers.end(), std::back_inserter(locMergeResult));
3623  locShowersByBunch.emplace(locBunchesSoFar, std::move(locMergeResult));
3624  Build_ParticleIndices(Gamma, locBunchesSoFar, locShowersByBunch[locBunchesSoFar], locVertexZBin);
3625  }
3626  return locShowersByBunch[locBeamBunches];
3627 }
3628 
3629 /******************************************************************* COMBO UTILITY FUNCTIONS ******************************************************************/
3630 
3631 void DSourceComboer::Register_ValidRFBunches(const DSourceComboUse& locSourceComboUse, const DSourceCombo* locSourceCombo, const vector<int>& locRFBunches, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding)
3632 {
3633  //search and register
3634  auto locComboInfo = std::get<2>(locSourceComboUse);
3635  dValidRFBunches_ByCombo.emplace(std::make_pair(locSourceCombo, std::get<1>(locSourceComboUse)), locRFBunches);
3636 
3637  //also, register for each individual bunch: so that we can get valid combos for some input rf bunches later
3638  if(locComboingStage != d_ChargedStage)
3639  {
3640  auto& locSourceCombosByBeamBunchByUse = Get_SourceCombosByBeamBunchByUse(dComboInfoChargeContent[locComboInfo], locChargedCombo_Presiding);
3641  auto& locCombosByBeamBunch = locSourceCombosByBeamBunchByUse[locSourceComboUse];
3642  for(const auto& locBeamBunch : locRFBunches)
3643  locCombosByBeamBunch[{locBeamBunch}].push_back(locSourceCombo);
3644  }
3645 }
3646 
3647 void DSourceComboer::Build_ComboResumeIndices(const DSourceComboUse& locSourceComboUse, ComboingStage_t locComboingStage, const DSourceCombo* locChargedCombo_Presiding)
3648 {
3649  auto locComboInfo = std::get<2>(locSourceComboUse);
3650 
3651  auto& locComboVector = *(Get_CombosSoFar(locComboingStage, dComboInfoChargeContent[locComboInfo], locChargedCombo_Presiding)[locSourceComboUse]);
3652  std::sort(locComboVector.begin(), locComboVector.end());
3653  Build_ComboIndices(locSourceComboUse, {}, locComboVector, locComboingStage);
3654 
3655  //register for each individual bunch: so that we can get valid combos for some input rf bunches later
3656  if(locComboingStage != d_ChargedStage)
3657  {
3658  auto& locSourceCombosByBeamBunchByUse = Get_SourceCombosByBeamBunchByUse(dComboInfoChargeContent[locComboInfo], locChargedCombo_Presiding);
3659  auto& locCombosByBeamBunch = locSourceCombosByBeamBunchByUse[locSourceComboUse];
3660  for(auto& locRFPair : locCombosByBeamBunch)
3661  {
3662  std::sort(locRFPair.second.begin(), locRFPair.second.end());
3663  Build_ComboIndices(locSourceComboUse, locRFPair.first, locRFPair.second, locComboingStage);
3664  }
3665  }
3666 }
3667 
3668 const vector<const DSourceCombo*>& DSourceComboer::Get_CombosForComboing(const DSourceComboUse& locComboUse, ComboingStage_t locComboingStage, const vector<int>& locBeamBunches, const DSourceCombo* locChargedCombo_PresidingPrevious)
3669 {
3670  if(dDebugLevel >= 20)
3671  {
3672  cout << "Get_CombosForComboing: stage, #bunches, charged combo, bunches " << locComboingStage << ", " << locBeamBunches.size() << ", " << locChargedCombo_PresidingPrevious << ", ";
3673  for(auto& locBunch : locBeamBunches)
3674  cout << locBunch << ", ";
3675  cout << endl;
3676  cout << "GET-COMBOS USE:" << endl;
3677  Print_SourceComboUse(locComboUse);
3678  }
3679 
3680  //THE INPUT locChargedCombo MUST BE:
3681  //Whatever charged combo PREVIOUSLY presided when creating the combos you're trying to get
3682  //find all combos for the given use that have an overlapping beam bunch with the input
3683  auto locChargeContent = dComboInfoChargeContent[std::get<2>(locComboUse)];
3684  if(locBeamBunches.empty() || (locChargeContent == d_Charged)) //e.g. fully charged, or a combo of 2 KLongs (RF bunches not saved for massive neutrals)
3685  return *((Get_CombosSoFar(locComboingStage, locChargeContent, locChargedCombo_PresidingPrevious))[locComboUse]);
3686 
3687  auto& locSourceCombosByBeamBunchByUse = Get_SourceCombosByBeamBunchByUse(locChargeContent, locChargedCombo_PresidingPrevious);
3688  auto locGroupBunchIterator = locSourceCombosByBeamBunchByUse[locComboUse].find(locBeamBunches);
3689  if(locGroupBunchIterator != locSourceCombosByBeamBunchByUse[locComboUse].end())
3690  return locGroupBunchIterator->second;
3691 
3692  return Get_CombosByBeamBunch(locComboUse, locSourceCombosByBeamBunchByUse[locComboUse], locBeamBunches, locComboingStage);
3693 }
3694 
3695 const vector<const DSourceCombo*>& DSourceComboer::Get_CombosByBeamBunch(const DSourceComboUse& locComboUse, DCombosByBeamBunch& locCombosByBunch, const vector<int>& locBeamBunches, ComboingStage_t locComboingStage)
3696 {
3697  if(dDebugLevel >= 20)
3698  {
3699  cout << "Get_CombosByBeamBunch: stage, # bunches, bunches: " << locComboingStage << ", " << locBeamBunches.size() << ", ";
3700  for(auto& locBunch : locBeamBunches)
3701  cout << locBunch << ", ";
3702  cout << endl;
3703  }
3704  if(locBeamBunches.empty())
3705  {
3706  Build_ComboIndices(locComboUse, locBeamBunches, locCombosByBunch[locBeamBunches], locComboingStage);
3707  return locCombosByBunch[{}];
3708  }
3709 
3710  //find all combos for the given use that have an overlapping beam bunch with the input
3711  //this shouldn't be called very many times per event, so we can be a little inefficient
3712  vector<int> locBunchesSoFar = {*locBeamBunches.begin()}; //0
3713  for(auto locBunchIterator = std::next(locBeamBunches.begin()); locBunchIterator != locBeamBunches.end(); ++locBunchIterator)
3714  {
3715  //get vectors and sort them: sort needed for union below
3716  auto& locCombosSoFar = locCombosByBunch[locBunchesSoFar]; //{0}
3717  auto& locBunchCombos = locCombosByBunch[{*locBunchIterator}]; //{1}
3718 
3719  locBunchesSoFar.push_back(*locBunchIterator); //{0, 1}
3720  if(locCombosByBunch.find(locBunchesSoFar) != locCombosByBunch.end())
3721  continue; //this subset already created and indexed
3722 
3723  if(locBunchCombos.empty())
3724  {
3725  locCombosByBunch.emplace(locBunchesSoFar, locCombosSoFar);
3726  Build_ComboIndices(locComboUse, locBeamBunches, locCombosByBunch[locBunchesSoFar], locComboingStage);
3727  continue;
3728  }
3729 
3730  //merge and move-emplace
3731  vector<const DSourceCombo*> locMergeResult;
3732  locMergeResult.reserve(locCombosSoFar.size() + locBunchCombos.size());
3733  std::set_union(locCombosSoFar.begin(), locCombosSoFar.end(), locBunchCombos.begin(), locBunchCombos.end(), std::back_inserter(locMergeResult));
3734  locCombosByBunch.emplace(locBunchesSoFar, std::move(locMergeResult)); //when building for 0, 1, 2 this replaces 0,1
3735  Build_ComboIndices(locComboUse, locBunchesSoFar, locCombosByBunch[locBunchesSoFar], locComboingStage);
3736  }
3737 
3738  return locCombosByBunch[locBeamBunches];
3739 }
3740 
3741 void DSourceComboer::Copy_ZIndependentMixedResults(const DSourceComboUse& locComboUseToCreate, const DSourceCombo* locChargedCombo_Presiding)
3742 {
3743  //Copy the results from the FCAL-only stage through to the both stage (that way we don't have to repeat them)
3744 
3745  //THE INPUT locChargedCombo MUST BE:
3746  //Whatever charged combo you are about to combo horizontally with to make this new, mixed combo
3747 
3748  //Get combos so far
3749  auto locChargeContent = dComboInfoChargeContent[std::get<2>(locComboUseToCreate)];
3750  auto& locSourceCombosByUseSoFar = Get_CombosSoFar(d_MixedStage_ZIndependent, locChargeContent, locChargedCombo_Presiding);
3751 
3752  //Get FCAL results
3753  auto locComboUseFCAL = Get_ZIndependentUse(locComboUseToCreate);
3754  if(dDebugLevel >= 20)
3755  {
3756  cout << "FCAL USE: " << endl;
3757  Print_SourceComboUse(locComboUseFCAL);
3758  }
3759  if(locSourceCombosByUseSoFar.find(locComboUseFCAL) == locSourceCombosByUseSoFar.end())
3760  return; //no results to copy, just return
3761  const auto& locFCALComboVector = *(locSourceCombosByUseSoFar[locComboUseFCAL]);
3762  if(dDebugLevel >= 20)
3763  cout << "copying " << locFCALComboVector.size() << " from the fcal vector" << endl;
3764  if(locFCALComboVector.empty())
3765  return;
3766 
3767  //Copy over the combos
3768  auto& locBothComboVector = *(locSourceCombosByUseSoFar[locComboUseToCreate]);
3769  locBothComboVector.reserve(locFCALComboVector.size() + dInitialComboVectorCapacity);
3770  locBothComboVector.assign(locFCALComboVector.begin(), locFCALComboVector.end());
3771 
3772  //Copy over the combos-by-beam-bunch
3773  auto& locSourceCombosByBeamBunchByUse = Get_SourceCombosByBeamBunchByUse(locChargeContent, locChargedCombo_Presiding);
3774  const auto& locCombosByBeamBunch = locSourceCombosByBeamBunchByUse[locComboUseFCAL];
3775  for(const auto& locComboBeamBunchPair : locCombosByBeamBunch)
3776  {
3777  if(locComboBeamBunchPair.first.size() == 1) //don't copy the overlap ones: they are not complete & need to be filled on the fly
3778  locSourceCombosByBeamBunchByUse[locComboUseToCreate].emplace(locComboBeamBunchPair);
3779  }
3780 }
3781 
3782 const DSourceCombo* DSourceComboer::Get_VertexPrimaryCombo(const DSourceCombo* locReactionCombo, const DReactionStepVertexInfo* locStepVertexInfo)
3783 {
3784  //if it's the production vertex, just return the input
3785  if(locStepVertexInfo->Get_ProductionVertexFlag())
3786  return locReactionCombo;
3787 
3788  //see if it's already been determined before: if so, just return it
3789  auto locCreationPair = std::make_pair(locReactionCombo, locStepVertexInfo);
3790  auto locIterator = dVertexPrimaryComboMap.find(locCreationPair);
3791  if(locIterator != dVertexPrimaryComboMap.end())
3792  return locIterator->second;
3793 
3794  //find it
3795  auto locReaction = locStepVertexInfo->Get_Reaction();
3796  auto locDesiredStepIndex = locStepVertexInfo->Get_StepIndices().front();
3797  auto locVertexPrimaryCombo = Get_StepSourceCombo(locReaction, locDesiredStepIndex, locReactionCombo, 0);
3798 
3799  //save it and return it
3800  dVertexPrimaryComboMap.emplace(locCreationPair, locVertexPrimaryCombo);
3801  return locVertexPrimaryCombo;
3802 }
3803 
3804 const DSourceCombo* DSourceComboer::Get_VertexPrimaryCombo(const DSourceCombo* locReactionCombo, const DReactionStepVertexInfo* locStepVertexInfo) const
3805 {
3806  //if it's the production vertex, just return the input
3807  if(locStepVertexInfo->Get_ProductionVertexFlag())
3808  return locReactionCombo;
3809 
3810  //see if it's already been determined before: if so, just return it
3811  auto locCreationPair = std::make_pair(locReactionCombo, locStepVertexInfo);
3812  auto locIterator = dVertexPrimaryComboMap.find(locCreationPair);
3813  if(locIterator != dVertexPrimaryComboMap.end())
3814  return locIterator->second;
3815 
3816  //find it
3817  auto locReaction = locStepVertexInfo->Get_Reaction();
3818  auto locDesiredStepIndex = locStepVertexInfo->Get_StepIndices().front();
3819  auto locVertexPrimaryCombo = Get_StepSourceCombo(locReaction, locDesiredStepIndex, locReactionCombo, 0);
3820 
3821  //return it
3822  return locVertexPrimaryCombo;
3823 }
3824 
3825 const DSourceCombo* DSourceComboer::Get_StepSourceCombo(const DReaction* locReaction, size_t locDesiredStepIndex, const DSourceCombo* locVertexPrimaryCombo, size_t locVertexPrimaryStepIndex) const
3826 {
3827  if(dDebugLevel >= 100)
3828  cout << "reaction, desired step index, current step index: " << locReaction->Get_ReactionName() << ", " << locDesiredStepIndex << ", " << locVertexPrimaryStepIndex << endl;
3829  if(locDesiredStepIndex == locVertexPrimaryStepIndex)
3830  return locVertexPrimaryCombo;
3831 
3832  //Get the list of steps we need to traverse //particle pair: step index, particle instance index
3833  vector<pair<size_t, int>> locParticleIndices = {std::make_pair(locDesiredStepIndex, DReactionStep::Get_ParticleIndex_Initial())};
3834  while(locParticleIndices.back().first != locVertexPrimaryStepIndex)
3835  {
3836  auto locParticlePair = DAnalysis::Get_InitialParticleDecayFromIndices(locReaction, locParticleIndices.back().first);
3837  if(dDebugLevel >= 100)
3838  cout << "decay from pair: " << locParticlePair.first << ", " << locParticlePair.second << endl;
3839  auto locStep = locReaction->Get_ReactionStep(locParticlePair.first);
3840  auto locInstanceIndex = DAnalysis::Get_ParticleInstanceIndex(locStep, locParticlePair.second);
3841  locParticleIndices.emplace_back(locParticlePair.first, locInstanceIndex);
3842  if(dDebugLevel >= 100)
3843  cout << "save indices: " << locParticlePair.first << ", " << locInstanceIndex << endl;
3844  }
3845 
3846  //start from back of locParticleIndices, searching
3847  while(true)
3848  {
3849  auto locNextStep = locParticleIndices[locParticleIndices.size() - 2].first;
3850  auto locInstanceIndexToFind = locParticleIndices.back().second;
3851  const auto& locUseToFind = dSourceComboUseReactionStepMap.find(locReaction)->second.find(locNextStep)->second;
3852  if(dDebugLevel >= 100)
3853  cout << "next step, instance to find, use to find: " << locNextStep << ", " << locInstanceIndexToFind << endl;
3854  if(dDebugLevel >= 100)
3855  Print_SourceComboUse(locUseToFind);
3856  locVertexPrimaryCombo = Find_Combo_AtThisStep(locVertexPrimaryCombo, locUseToFind, locInstanceIndexToFind);
3857  if(dDebugLevel >= 100)
3858  cout << "pointer = " << locVertexPrimaryCombo << endl;
3859  if(locVertexPrimaryCombo == nullptr)
3860  return nullptr; //e.g. entirely neutral step when input is charged
3861  if(locNextStep == locDesiredStepIndex)
3862  return locVertexPrimaryCombo;
3863  locParticleIndices.pop_back();
3864  }
3865 
3866  return nullptr;
3867 }
3868 
3869 const DSourceCombo* DSourceComboer::Find_Combo_AtThisStep(const DSourceCombo* locSourceCombo, DSourceComboUse locUseToFind, size_t locDecayInstanceIndex) const
3870 {
3871  //ignores z-bin when comparing
3872  //if z-dependent, go to z-independent use
3873  if(std::get<1>(locUseToFind) != DSourceComboInfo::Get_VertexZIndex_ZIndependent())
3874  locUseToFind = dZDependentUseToIndependentMap.find(locUseToFind)->second;
3875  if(dDebugLevel >= 100)
3876  {
3877  cout << "Find_Combo_AtThisStep: USE TO FIND:" << endl;
3878  DAnalysis::Print_SourceComboUse(locUseToFind);
3879  }
3880  for(const auto& locDecayPair : locSourceCombo->Get_FurtherDecayCombos())
3881  {
3882  //if z-dependent, go to z-independent
3883  auto locDecayUse = locDecayPair.first;
3884  if(std::get<1>(locDecayUse) != DSourceComboInfo::Get_VertexZIndex_ZIndependent())
3885  locDecayUse = dZDependentUseToIndependentMap.find(locDecayUse)->second;
3886  if(dDebugLevel >= 100)
3887  {
3888  cout << "USE TO CHECK:" << endl;
3889  DAnalysis::Print_SourceComboUse(locDecayUse);
3890  }
3891 
3892  if(locDecayUse == locUseToFind) //good, do stuff
3893  return locDecayPair.second[locDecayInstanceIndex];
3894  if(std::get<0>(locDecayUse) != Unknown)
3895  continue; //is another step!
3896 
3897  //vector of combos is guaranteed to be size 1, and it's guaranteed that none of ITS further decays are unknown
3898  auto locComboToSearch = locDecayPair.second[0];
3899  if(dDebugLevel >= 100)
3900  cout << "#to-check decay uses: " << locComboToSearch->Get_FurtherDecayCombos().size() << endl;
3901  for(const auto& locNestedDecayPair : locComboToSearch->Get_FurtherDecayCombos())
3902  {
3903  //if z-dependent, go to z-independent
3904  auto locNestedDecayUse = locNestedDecayPair.first;
3905  if(std::get<1>(locNestedDecayUse) != DSourceComboInfo::Get_VertexZIndex_ZIndependent())
3906  locNestedDecayUse = dZDependentUseToIndependentMap.find(locNestedDecayUse)->second;
3907  if(dDebugLevel >= 100)
3908  {
3909  cout << "NESTED USE TO CHECK:" << endl;
3910  DAnalysis::Print_SourceComboUse(locNestedDecayUse);
3911  }
3912  if(locNestedDecayUse == locUseToFind) //good, do stuff
3913  return locNestedDecayPair.second[locDecayInstanceIndex];
3914  }
3915  }
3916 
3917  //Not found: Either invalid request, OR the input is a fully-charged combo being used for a Use that contains neutrals (created during charged-only stage): Return the input, it is already what you want
3918  return locSourceCombo;
3919 }
3920 
3921 pair<DSourceComboUse, size_t> DSourceComboer::Get_StepSourceComboUse(const DReaction* locReaction, size_t locDesiredStepIndex, DSourceComboUse locVertexPrimaryComboUse, size_t locVertexPrimaryStepIndex) const
3922 {
3923  //size_t: combo instance
3924  if(dDebugLevel >= 100)
3925  cout << "reaction, desired step index, current step index: " << locReaction->Get_ReactionName() << ", " << locDesiredStepIndex << ", " << locVertexPrimaryStepIndex << endl;
3926  if(locDesiredStepIndex == locVertexPrimaryStepIndex)
3927  return std::make_pair(locVertexPrimaryComboUse, size_t(1));
3928 
3929  //Get the list of steps we need to traverse //particle pair: step index, particle instance index
3930  vector<pair<size_t, int>> locParticleIndices = {std::make_pair(locDesiredStepIndex, DReactionStep::Get_ParticleIndex_Initial())};
3931  while(locParticleIndices.back().first != locVertexPrimaryStepIndex)
3932  {
3933  auto locParticlePair = DAnalysis::Get_InitialParticleDecayFromIndices(locReaction, locParticleIndices.back().first);
3934  if(dDebugLevel >= 100)
3935  cout << "decay from pair: " << locParticlePair.first << ", " << locParticlePair.second << endl;
3936  auto locStep = locReaction->Get_ReactionStep(locParticlePair.first);
3937  auto locInstanceIndex = DAnalysis::Get_ParticleInstanceIndex(locStep, locParticlePair.second);
3938  locParticleIndices.emplace_back(locParticlePair.first, locInstanceIndex);
3939  if(dDebugLevel >= 100)
3940  cout << "save indices: " << locParticlePair.first << ", " << locInstanceIndex << endl;
3941  }
3942 
3943  //start from back of locParticleIndices, searching
3944  while(true)
3945  {
3946  auto locNextStep = locParticleIndices[locParticleIndices.size() - 2].first;
3947  auto locInstanceIndexToFind = locParticleIndices.back().second;
3948  const auto& locUseToFind = dSourceComboUseReactionStepMap.find(locReaction)->second.find(locNextStep)->second;
3949  if(dDebugLevel >= 100)
3950  cout << "next step, instance to find, use to find: " << locNextStep << ", " << locInstanceIndexToFind << endl;
3951  if(dDebugLevel >= 100)
3952  Print_SourceComboUse(locUseToFind);
3953  locVertexPrimaryComboUse = Find_ZDependentUse_AtThisStep(locVertexPrimaryComboUse, locUseToFind, locInstanceIndexToFind);
3954  if(std::get<2>(locVertexPrimaryComboUse) == nullptr)
3955  return std::make_pair(locVertexPrimaryComboUse, size_t(locInstanceIndexToFind + 1)); //e.g. entirely neutral step when input is charged
3956  if(locNextStep == locDesiredStepIndex)
3957  return std::make_pair(locVertexPrimaryComboUse, size_t(locInstanceIndexToFind + 1));
3958  locParticleIndices.pop_back();
3959  }
3960  return std::make_pair(DSourceComboUse(Unknown, 0, nullptr, 0, Unknown), size_t(1));
3961 }
3962 
3963 DSourceComboUse DSourceComboer::Find_ZDependentUse_AtThisStep(const DSourceComboUse& locSourceComboUse, DSourceComboUse locUseToFind, size_t locDecayInstanceIndex) const
3964 {
3965  //ignores z-bin when comparing
3966  //if z-dependent, go to z-independent use
3967  if(std::get<1>(locUseToFind) != DSourceComboInfo::Get_VertexZIndex_ZIndependent())
3968  locUseToFind = dZDependentUseToIndependentMap.find(locUseToFind)->second;
3969  if(dDebugLevel >= 100)
3970  {
3971  cout << "Find_Combo_AtThisStep: USE TO FIND:" << endl;
3972  DAnalysis::Print_SourceComboUse(locUseToFind);
3973  }
3974  for(const auto& locDecayPair : std::get<2>(locSourceComboUse)->Get_FurtherDecays())
3975  {
3976  //if z-dependent, go to z-independent
3977  auto locDecayUse = locDecayPair.first;
3978  auto locZIndependentDecayUse = locDecayUse;
3979  if(std::get<1>(locDecayUse) != DSourceComboInfo::Get_VertexZIndex_ZIndependent())
3980  locZIndependentDecayUse = dZDependentUseToIndependentMap.find(locDecayUse)->second;
3981  if(dDebugLevel >= 100)
3982  {
3983  cout << "USE TO CHECK:" << endl;
3984  DAnalysis::Print_SourceComboUse(locDecayUse);
3985  }
3986 
3987  if(locZIndependentDecayUse == locUseToFind) //good, do stuff
3988  return locDecayUse;
3989  if(std::get<0>(locDecayUse) != Unknown)
3990  continue; //is another step!
3991 
3992  //check other uses at this step (further depth guaranteed to be only 1)
3993  if(dDebugLevel >= 100)
3994  cout << "#to-check decay uses: " << std::get<2>(locDecayUse)->Get_FurtherDecays().size() << endl;
3995  for(const auto& locNestedDecayPair : std::get<2>(locDecayUse)->Get_FurtherDecays())
3996  {
3997  //if z-dependent, go to z-independent
3998  auto locNestedDecayUse = locNestedDecayPair.first;
3999  auto locZIndependentNestedDecayUse = locNestedDecayUse;
4000  if(std::get<1>(locNestedDecayUse) != DSourceComboInfo::Get_VertexZIndex_ZIndependent())
4001  locZIndependentNestedDecayUse = dZDependentUseToIndependentMap.find(locNestedDecayUse)->second;
4002  if(dDebugLevel >= 100)
4003  {
4004  cout << "NESTED USE TO CHECK:" << endl;
4005  DAnalysis::Print_SourceComboUse(locZIndependentNestedDecayUse);
4006  }
4007  if(locZIndependentNestedDecayUse == locUseToFind) //good, do stuff
4008  return locNestedDecayUse;
4009  }
4010  }
4011 
4012  //Not found: Either invalid request, OR the input is a fully-charged combo being used for a Use that contains neutrals (created during charged-only stage): Return the input, it is already what you want
4013  return DSourceComboUse(Unknown, 0, nullptr, 0, Unknown);
4014 }
4015 /*
4016  * For K0, Sigma+, p the full combos will be:
4017  * 0: X -> A, 1, 3 (mixed -> charged, mixed, mixed)
4018  * A: X -> p (charged)
4019  * 1: K0 -> B, 2 (mixed -> charged, neutral)
4020  * B: X -> pi+, pi- (charged)
4021  * 2: pi0 -> 2g (neutral)
4022  * 3: Sigma+ -> C, n (mixed -> charged, n)
4023  * C: X -> pi+ (charged)
4024  *
4025  * For XXX, the charged combos will be:
4026  * 0: X -> A, B, C
4027  * A: X -> p
4028  * B: X -> pi+, pi-
4029  * C: X -> pi+
4030  *
4031  * For XXX, the presiding/withnow combos will be:
4032  * 0: X -> A, 1, 3 (mixed -> charged, mixed, mixed) //presiding = 0, withnow = A
4033  * A: X -> p (charged) //both = nullptr
4034  * 1: K0 -> B, 2 (mixed -> charged, neutral) //presiding = 0, withnow = B
4035  * B: X -> pi+, pi- (charged) //both = nullptr
4036  * 2: pi0 -> 2g (neutral) //both = nullptr
4037  * 3: Sigma+ -> C, n (mixed -> charged, n) //presiding = 0, withnow = C
4038  * C: X -> pi+ (charged) //both = nullptr
4039  *
4040 */
4041 
4042 const DSourceCombo* DSourceComboer::Get_ChargedCombo_WithNow(const DSourceCombo* locChargedCombo_Presiding, const DSourceComboInfo* locToCreateComboInfo, ComboingStage_t locComboingStage) const
4043 {
4044  if(locChargedCombo_Presiding == nullptr)
4045  return nullptr;
4046 
4047  //find the charged use what use we want
4048  DSourceComboUse locWithNowComboUse{Unknown, 0, nullptr, false, Unknown};
4049  for(const auto& locDecayComboPair : locToCreateComboInfo->Get_FurtherDecays())
4050  {
4051  if(Get_ChargeContent(std::get<2>(locDecayComboPair.first)) != d_Charged)
4052  continue;
4053  locWithNowComboUse = locDecayComboPair.first;
4054  break;
4055  }
4056 
4057  if(std::get<2>(locWithNowComboUse) == nullptr)
4058  {
4059  if(dDebugLevel >= 20)
4060  cout << "CHARGED COMBO WITH NOW: Same as presiding." << endl;
4061  return locChargedCombo_Presiding; //the info we are trying to create will use the presiding itself
4062  }
4063  return Get_NextChargedCombo(locChargedCombo_Presiding, locWithNowComboUse, locComboingStage, false, 0);
4064 }
4065 
4066 const DSourceCombo* DSourceComboer::Get_NextChargedCombo(const DSourceCombo* locChargedCombo_Presiding, const DSourceComboUse& locNextComboUse, ComboingStage_t locComboingStage, bool locGetPresidingFlag, size_t locInstance) const
4067 {
4068  //locInstance starts from ONE!! (and is not used if locGetPresidingFlag = false (getting withnow))
4069  if(locComboingStage == d_ChargedStage)
4070  return nullptr;
4071  if(locChargedCombo_Presiding == nullptr)
4072  return nullptr;
4073  if(Get_ChargeContent(std::get<2>(locNextComboUse)) == d_Neutral)
4074  return nullptr; //not needed
4075 
4076  auto locFurtherDecayCombos = locChargedCombo_Presiding->Get_FurtherDecayCombos();
4077 
4078  auto locUseToFind = (locComboingStage == d_MixedStage_ZIndependent) ? locNextComboUse : dZDependentUseToIndependentMap.find(locNextComboUse)->second;
4079  auto locIteratorPair = std::equal_range(locFurtherDecayCombos.begin(), locFurtherDecayCombos.end(), locUseToFind, DSourceCombo::DCompare_FurtherDecays());
4080 
4081  if(dDebugLevel >= 20)
4082  {
4083  cout << "Get_NextChargedCombo: get-presiding-flag, instance, stage, find result, use-to-find: " << locGetPresidingFlag << ", " << locInstance << ", " << locComboingStage << ", " << (locIteratorPair.first != locIteratorPair.second) << endl;
4084  DAnalysis::Print_SourceComboUse(locUseToFind);
4085  cout << "Presiding combo:" << endl;
4086  DAnalysis::Print_SourceCombo(locChargedCombo_Presiding);
4087  }
4088 
4089  //check if the use you are looking for is a temporary (e.g. vertical grouping of 2KShorts when comboing horizontally)
4090  //or if the charged combos were supposed to be comboed with neutrals, but were instead promoted: no intermediary charged combo, just retuern current
4091  if(locIteratorPair.first == locIteratorPair.second)
4092  return locChargedCombo_Presiding; //temporary: the presiding is still the same!
4093 
4094  //get the vector of potential charged combos
4095  auto locNextChargedComboVector = (*locIteratorPair.first).second;
4096  if(dDebugLevel >= 20)
4097  cout << "next charged combo vector size: = " << locNextChargedComboVector.size() << endl;
4098 
4099  //if getting with-now, size is guaranteed to be 1, just get the first one
4100  if(!locGetPresidingFlag)
4101  {
4102  if(dDebugLevel >= 20)
4103  {
4104  cout << "CHARGED COMBO WITH NOW:" << endl;
4105  Print_SourceCombo(locNextChargedComboVector[0]);
4106  }
4107  return locNextChargedComboVector[0];
4108  }
4109 
4110  //if on z-independent, don't need to do anything fancy, just return the requested instance
4111  if(locComboingStage == d_MixedStage_ZIndependent)
4112  return locNextChargedComboVector[locInstance - 1];
4113 
4114  //there might be multiple combos (e.g. K0 decays), each at a different vertex-z
4115  //so, we must retrieve the N'th charged combo with the correct vertex-z bin
4116  size_t locCount = 0;
4117  auto locDesiredVertexZBin = std::get<1>(locNextComboUse);
4118  for(const auto& locNextPotentialCombo : locNextChargedComboVector)
4119  {
4120  //either locNextPotentialCombo is the primary combo of a detached vertex (and we need to check the zbin), or it's not (returns -1) and we don't
4121  if(IsDetachedVertex(std::get<0>(locUseToFind)))
4122  {
4123  auto locVertexChargeContent = DAnalysis::Get_ChargeContent_ThisVertex(std::get<2>(locUseToFind));
4124  auto locIsCombo2ndVertex = (locVertexChargeContent == d_Neutral);
4125  auto locIsVertexKnown = dSourceComboVertexer->Get_IsVertexKnown_NoBeam(false, locNextPotentialCombo, locIsCombo2ndVertex);
4126  auto locNextVertexZBin = dSourceComboVertexer->Get_VertexZBin_NoBeam(false, locNextPotentialCombo, locIsCombo2ndVertex); //defaults to center of target if not known
4127 
4128  if(dDebugLevel >= 20)
4129  cout << "detached next potential combo, next zbin, desired zbin = " << locNextPotentialCombo << ", " << int(locNextVertexZBin) << ", " << int(locDesiredVertexZBin) << endl;
4130  //if desired = independent we don't care, or if unknown we don't need to check
4131  if(locIsVertexKnown && (locNextVertexZBin != locDesiredVertexZBin) && (locDesiredVertexZBin != DSourceComboInfo::Get_VertexZIndex_ZIndependent()))
4132  continue;
4133  }
4134  if(dDebugLevel >= 20)
4135  cout << "pre-count, instance = " << locCount << ", " << locInstance << endl;
4136  if(++locCount == locInstance)
4137  return locNextPotentialCombo;
4138  }
4139 
4140  if(dDebugLevel >= 20)
4141  cout << "uh oh" << endl;
4142  return nullptr; //uh oh ...
4143 }
4144 
4145 bool DSourceComboer::Get_PromoteFlag(ComboingStage_t locComboingStage, Particle_t locDecayPID_UseToCheck, const DSourceComboInfo* locComboInfo_UseToCreate, const DSourceComboInfo* locComboInfo_UseToCheck, DSourceComboUse& locNonNeutralUse) const
4146 {
4147  locNonNeutralUse = DSourceComboUse{Unknown, 0, nullptr, false, Unknown};
4148  if(locDecayPID_UseToCheck != Unknown)
4149  return false;
4150 
4151  auto locFurtherDecayInfo_UseToCheck = locComboInfo_UseToCheck->Get_FurtherDecays();
4152 
4153  //don't promote if is mixed and merely contains a promoted charged combo
4154  if((locComboingStage == d_ChargedStage) && (Get_ChargeContent(locComboInfo_UseToCheck) == d_AllCharges))
4155  {
4156  size_t locNumNeutralUses = locComboInfo_UseToCheck->Get_NumParticles().size();
4157  size_t locNumNonNeutralUses = 0;
4158  for(auto& locAllBut1DecayPair : locFurtherDecayInfo_UseToCheck)
4159  {
4160  if(Get_ChargeContent(std::get<2>(locAllBut1DecayPair.first)) == d_Neutral)
4161  ++locNumNeutralUses;
4162  else
4163  {
4164  ++locNumNonNeutralUses;
4165  locNonNeutralUse = locAllBut1DecayPair.first;
4166  }
4167  }
4168  if((locNumNeutralUses >= 1) && (locNumNonNeutralUses == 1))
4169  return false; //merely a promoted charged combo
4170  locNonNeutralUse = DSourceComboUse{Unknown, 0, nullptr, false, Unknown}; //reset in case > 1
4171  }
4172 
4173 //we must: ungroup all-but-1 use: save the existing combo under the charged/mixed use & ditch the neutral decay uses
4174 //in this case: it becomes a no-promote, but a different use
4175 
4176  if(!locFurtherDecayInfo_UseToCheck.empty())
4177  {
4178  auto locFurtherDecayInfo_UseToCreate = locComboInfo_UseToCreate->Get_FurtherDecays();
4179  return std::binary_search(locFurtherDecayInfo_UseToCreate.begin(), locFurtherDecayInfo_UseToCreate.end(), locFurtherDecayInfo_UseToCheck.front(), DSourceComboInfo::DCompare_FurtherDecays());
4180  }
4181  else
4182  {
4183  auto locNumParticles_ToAdd = locComboInfo_UseToCheck->Get_NumParticles();
4184  auto locNumParticles_UseToCreate = locComboInfo_UseToCreate->Get_NumParticles();
4185  return std::binary_search(locNumParticles_UseToCreate.begin(), locNumParticles_UseToCreate.end(), locNumParticles_ToAdd.front(), DSourceComboInfo::DCompare_ParticlePairPIDs());
4186  }
4187 }
4188 
4189 bool DSourceComboer::Check_Reactions(vector<const DReaction*>& locReactions)
4190 {
4191  //All of the reactions in the vertex-info are guaranteed to have the same channel content
4192  //They just may differ in actions, or skims
4193  //So, we can check #particles for just one reaction, but must check skims for all reactions
4194  if(!Check_NumParticles(locReactions.front()))
4195  {
4196  if(dDebugLevel > 0)
4197  cout << "Not enough particles: No combos." << endl;
4198  return false;
4199  }
4200  for(auto& locReaction : locReactions)
4201  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::Min_Particles] = 1; //is really #-events
4202 
4203  //Check Max neutrals
4204  auto locNumNeutralNeeded = locReactions.front()->Get_FinalPIDs(-1, false, false, d_Neutral, true).size(); //no missing, no decaying, include duplicates
4205  auto locNumDetectedShowers = dShowersByBeamBunchByZBin[DSourceComboInfo::Get_VertexZIndex_Unknown()][{}].size();
4206  if((locNumNeutralNeeded > 0) && (locNumDetectedShowers > dMaxNumNeutrals))
4207  {
4208  if(dDebugLevel > 0)
4209  cout << "Too many neutrals: No combos." << endl;
4210  return false;
4211  }
4212  //Check additional showers
4213  auto NumShower_Checker = [&](const DReaction* locReaction) -> bool
4214  {
4215  auto locCutPair = locReaction->Get_MaxExtraShowers();
4216  if(!locCutPair.first)
4217  return false;
4218  return ((locNumDetectedShowers - locNumNeutralNeeded) > locCutPair.second);
4219  };
4220  locReactions.erase(std::remove_if(locReactions.begin(), locReactions.end(), NumShower_Checker), locReactions.end());
4221  if(locReactions.empty())
4222  {
4223  if(dDebugLevel > 0)
4224  cout << "Too many showers (" << locNumDetectedShowers << "): No combos." << endl;
4225  return false;
4226  }
4227 
4228  //Check Max charged tracks
4229  auto locNumTracksNeeded = locReactions.front()->Get_FinalPIDs(-1, false, false, d_Charged, true).size(); //no missing, no decaying, include duplicates
4230  auto NumExtra_Checker = [&](const DReaction* locReaction) -> bool
4231  {
4232  auto locCutPair = locReaction->Get_MaxExtraGoodTracks();
4233  if(!locCutPair.first)
4234  return false;
4235  return ((dNumChargedTracks - locNumTracksNeeded) > locCutPair.second);
4236  };
4237  locReactions.erase(std::remove_if(locReactions.begin(), locReactions.end(), NumExtra_Checker), locReactions.end());
4238  if(locReactions.empty())
4239  {
4240  if(dDebugLevel > 0)
4241  cout << "Too many tracks (" << dNumChargedTracks << "): No combos." << endl;
4242  return false;
4243  }
4244  for(auto& locReaction : locReactions)
4245  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::Max_Particles] = 1; //is really #-events
4246 
4247  //Check skims
4248  auto Skim_Checker = [this](const DReaction* locReaction) -> bool{return !Check_Skims(locReaction);};
4249  locReactions.erase(std::remove_if(locReactions.begin(), locReactions.end(), Skim_Checker), locReactions.end());
4250  if(locReactions.empty())
4251  {
4252  if(dDebugLevel > 0)
4253  cout << "Event not in skim: No combos." << endl;
4254  return false;
4255  }
4256  for(auto& locReaction : locReactions)
4257  dNumCombosSurvivedStageTracker[locReaction][DConstructionStage::In_Skim] = 1; //is really #-events
4258 
4259  return true;
4260 }
4261 
4263 {
4264  if(dDebugLevel > 0)
4265  cout << "Checking #particles" << endl;
4266 
4267  //see if enough particles were detected to build this reaction
4268  auto locReactionPIDs = locReaction->Get_FinalPIDs(-1, false, false, d_AllCharges, true); //no missing, no decaying, include duplicates
4269  auto locPIDMap = DAnalysis::Convert_VectorToCountMap<Particle_t>(locReactionPIDs);
4270  size_t locNumPositiveNeeded = 0, locNumNegativeNeeded = 0, locNumNeutralNeeded = 0;
4271  for(const auto& locPIDPair : locPIDMap)
4272  {
4273  if(ParticleCharge(locPIDPair.first) > 0)
4274  locNumPositiveNeeded += locPIDPair.second;
4275  else if(ParticleCharge(locPIDPair.first) < 0)
4276  locNumNegativeNeeded += locPIDPair.second;
4277  else
4278  locNumNeutralNeeded += locPIDPair.second;
4279  }
4280  auto locNumDetectedShowers = dShowersByBeamBunchByZBin[DSourceComboInfo::Get_VertexZIndex_Unknown()][{}].size();
4281 
4282  //check by charge
4283  if(dDebugLevel > 0)
4284  cout << "q+: Need " << locNumPositiveNeeded << ", Have " << dTracksByCharge[true].size() << endl;
4285  if(dTracksByCharge[true].size() < locNumPositiveNeeded)
4286  return false;
4287  if(dDebugLevel > 0)
4288  cout << "q-: Need " << locNumNegativeNeeded << ", Have " << dTracksByCharge[false].size() << endl;
4289  if(dTracksByCharge[false].size() < locNumNegativeNeeded)
4290  return false;
4291  if(dDebugLevel > 0)
4292  cout << "q+/-: Need " << locNumNegativeNeeded + locNumPositiveNeeded << ", Have " << dNumChargedTracks << endl;
4293  if(dNumChargedTracks < (locNumNegativeNeeded + locNumPositiveNeeded))
4294  return false;
4295  if(dDebugLevel > 0)
4296  cout << "q0: Need " << locNumNeutralNeeded << ", Have " << locNumDetectedShowers << ", Max allowed: " << dMaxNumNeutrals << endl;
4297  if(locNumDetectedShowers < locNumNeutralNeeded)
4298  return false;
4299 
4300  for(const auto& locPIDPair : locPIDMap)
4301  {
4302  auto locNumParticlesForComboing = Get_ParticlesForComboing(locPIDPair.first, d_MixedStage).size();
4303  if(dDebugLevel > 0)
4304  cout << ParticleType(locPIDPair.first) << ": Need " << locPIDPair.second << ", Have " << locNumParticlesForComboing << endl;
4305  if(locNumParticlesForComboing < locPIDPair.second)
4306  return false;
4307  if(locPIDPair.first != Gamma)
4308  continue;
4309 
4310  //check if these photons can even at least agree on a beam bunch, regardless of vertex position
4311  size_t locMaxNumPhotonsSameBunch = 0;
4312  for(const auto& locZBinPair : dShowersByBeamBunchByZBin) //loop over z-bins
4313  {
4314  for(const auto& locBunchPair : locZBinPair.second) //loop over bunches
4315  {
4316  if(locBunchPair.first.empty())
4317  continue;
4318  if(locBunchPair.second.size() > locMaxNumPhotonsSameBunch)
4319  locMaxNumPhotonsSameBunch = locBunchPair.second.size();
4320  }
4321  }
4322  if(dDebugLevel > 0)
4323  cout << ParticleType(locPIDPair.first) << ": Need " << locPIDPair.second << ", Have at most " << locMaxNumPhotonsSameBunch << " that agree on any beam bunch." << endl;
4324  if(locMaxNumPhotonsSameBunch < locPIDPair.second)
4325  return false;
4326  }
4327  return true;
4328 }
4329 
4331 {
4332  cout << "Num combos by use (charged):" << endl;
4333  for(const auto& locCombosByUsePair : dSourceCombosByUse_Charged)
4334  {
4335  cout << locCombosByUsePair.second->size() << " of ";
4336  Print_SourceComboUse(locCombosByUsePair.first);
4337 
4338  //save
4339  auto locIterator = dNumMixedCombosMap_Charged.find(locCombosByUsePair.first);
4340  if(locIterator == dNumMixedCombosMap_Charged.end())
4341  dNumMixedCombosMap_Charged.emplace(locCombosByUsePair.first, locCombosByUsePair.second->size());
4342  else
4343  locIterator->second += locCombosByUsePair.second->size();
4344  }
4345 
4346  //get #mixed by use (must merge results for different charged combos)
4347  map<DSourceComboUse, size_t> locNumMixedCombosMap;
4348  for(const auto& locChargedComboPair : dMixedCombosByUseByChargedCombo)
4349  {
4350  const auto& locCombosByUseMap = locChargedComboPair.second;
4351  for(const auto& locCombosByUsePair : locCombosByUseMap)
4352  {
4353  auto locIterator = locNumMixedCombosMap.find(locCombosByUsePair.first);
4354  if(locIterator == locNumMixedCombosMap.end())
4355  locNumMixedCombosMap.emplace(locCombosByUsePair.first, locCombosByUsePair.second->size());
4356  else
4357  locIterator->second += locCombosByUsePair.second->size();
4358  }
4359  }
4360 
4361  cout << "Num combos by use (neutral/mixed):" << endl;
4362  for(const auto& locNumCombosByUsePair : locNumMixedCombosMap)
4363  {
4364  cout << locNumCombosByUsePair.second << " of ";
4365  Print_SourceComboUse(locNumCombosByUsePair.first);
4366 
4367  //save
4368  auto locIterator = dNumMixedCombosMap_Mixed.find(locNumCombosByUsePair.first);
4369  if(locIterator == dNumMixedCombosMap_Mixed.end())
4370  dNumMixedCombosMap_Mixed.emplace(locNumCombosByUsePair.first, locNumCombosByUsePair.second);
4371  else
4372  locIterator->second += locNumCombosByUsePair.second;
4373  }
4374 }
4375 
4376 void DSourceComboer::Check_ForDuplicates(const vector<const DSourceCombo*>& locCombos) const
4377 {
4378  //Check for dupes & reuses!
4379  if(std::any_of(locCombos.begin(), locCombos.end(), DSourceComboChecker_ReusedParticle()))
4380  {
4381  cout << "Re-used particles, event = " << dEventNumber << ". Aborting!" << endl;
4382  abort();
4383  }
4384  for(size_t loc_i = 0; loc_i < locCombos.size(); ++loc_i)
4385  {
4386  for(size_t loc_j = loc_i + 1; loc_j < locCombos.size(); ++loc_j)
4387  {
4388  if(!DAnalysis::Check_AreDuplicateCombos(locCombos[loc_i], locCombos[loc_j]))
4389  continue;
4390  cout << "Duplicate particles, event = " << dEventNumber << ". Aborting!" << endl;
4391  cout << "DUPE COMBO 1:" << endl;
4392  Print_SourceCombo(locCombos[loc_i]);
4393  cout << "DUPE COMBO 2:" << endl;
4394  Print_SourceCombo(locCombos[loc_j]);
4395  abort();
4396  }
4397  }
4398 }
4399 
4400 } //end DAnalysis namespace
void Combo_WithNeutralsAndBeam(const vector< const DReaction * > &locReactions, const DReactionVertexInfo *locReactionVertexInfo, const DSourceComboUse &locPrimaryComboUse, const DSourceCombo *locReactionChargedCombo, const vector< int > &locBeamBunches_Charged, DCombosByReaction &locOutputComboMap)
bool Cut_InvariantMass_MissingMassVertex(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionFullCombo, const DKinematicData *locBeamParticle, int locRFBunch)
pair< DSourceComboUse, size_t > Get_StepSourceComboUse(const DReaction *locReaction, size_t locDesiredStepIndex, DSourceComboUse locVertexPrimaryComboUse, size_t locVertexPrimaryStepIndex) const
int Get_ParticleIndex(const DReactionStep *locStep, Particle_t locInputPID, size_t locInstance)
bool Get_IsFirstStepBeam(const DReaction *locReaction)
Definition: DReaction.h:178
void Combo_Horizontally_All(const DSourceComboUse &locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, unsigned char locNumTabs)
double getEnergy() const
Definition: DFCALShower.h:156
bool Select_RFBunches_PhotonVertices(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionFullCombo, vector< int > &locValidRFBunches)
void Set_SourceComboTimeHandler(const DSourceComboTimeHandler *locSourceComboTimeHandler)
DSourceCombosByUse_Large dSourceCombosByUse_Charged
set< const DSourceComboInfo *, DCompare_SourceComboInfos > dSourceComboInfoSet
const vector< const DSourceCombo * > & Get_CombosByBeamBunch(const DSourceComboUse &locComboUse, DCombosByBeamBunch &locCombosByBunch, const vector< int > &locBeamBunches, ComboingStage_t locComboingStage)
map< tuple< const JObject *, Particle_t, vector< int >, signed char >, size_t > dResumeSearchAfterIndices_Particles
signed char Get_VertexZBin_NoBeam(bool locIsProductionVertex, const DSourceCombo *locPrimaryVertexCombo, bool locIsCombo2ndVertex) const
int Select_RFBunch_Full(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionFullCombo, const vector< int > &locRFBunches)
static char * ParticleName_ROOT(Particle_t p)
Definition: particleType.h:851
const DReaction * Get_Reaction(void) const
const DSourceComboInfo * GetOrMake_SourceComboInfo(const vector< pair< Particle_t, unsigned char >> &locNumParticles, const vector< pair< DSourceComboUse, unsigned char >> &locFurtherDecays, unsigned char locNumTabs)
vector< int > Get_CommonRFBunches(const vector< int > &locRFBunches1, const JObject *locObject, signed char locVertexZBin) const
DSourceComboP4Handler * dSourceComboP4Handler
static constexpr int Get_ParticleIndex_Initial(void)
map< vector< int >, vector< const DSourceCombo * >> DCombosByBeamBunch
Charge_t Get_ChargeContent_ThisVertex(const DSourceComboInfo *locSourceComboInfo)
Definition: DSourceCombo.h:503
void Setup(const vector< const DNeutralShower * > &locNeutralShowers, const DEventRFBunch *locInitialEventRFBunch, const DDetectorMatches *locDetectorMatches)
pair< int, int > Get_InitialParticleDecayFromIndices(const DReaction *locReaction, int locStepIndex)
Definition: DReaction.cc:93
bool Get_HasPhotons(const DSourceComboInfo *locComboInfo)
Definition: DSourceCombo.h:521
unordered_map< const DReaction *, map< size_t, DSourceComboUse > > dSourceComboUseReactionStepMap
bool Cut_InvariantMass_HasMassiveNeutral_OrPhotonVertex(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionFullCombo, vector< int > &locValidRFBunches)
map< vector< int >, vector< const JObject * >> DPhotonShowersByBeamBunch
void Create_Combo_OneParticle(const DSourceComboUse &locComboUseToCreate, ComboingStage_t locComboingStage, unsigned char locNumTabs)
DSourceComboUse Find_ZDependentUse_AtThisStep(const DSourceComboUse &locSourceComboUse, DSourceComboUse locUseToFind, size_t locDecayInstanceIndex) const
bool Check_AreDuplicateCombos(const DSourceCombo *lhs, const DSourceCombo *rhs)
Definition: DSourceCombo.h:559
map< Particle_t, map< DetectorSystem_t, vector< double > > > dEOverPCuts_TF1Params
map< Particle_t, map< DetectorSystem_t, TH2 * > > dHistMap_dEdx
Definition: GlueX.h:16
map< const DReaction *, TH1 * > dNumCombosSurvivedStageMap
size_t Get_ResumeAtIndex_Combos(const DSourceComboUse &locSourceComboUse, const DSourceCombo *locPreviousCombo, const vector< int > &locBeamBunches, ComboingStage_t locComboingStage) const
unordered_set< const DSourceComboInfo * > dComboInfosWithPhotons
const DSourceComboInfo * MakeOrGet_SourceComboInfo(const vector< pair< Particle_t, unsigned char >> &locNumParticles, const vector< pair< DSourceComboUse, unsigned char >> &locFurtherDecays, unsigned char locNumTabs)
vector< const DReaction * > Get_Reactions(JEventLoop *locEventLoop)
Definition: DReaction.cc:218
void Build_ComboIndices(const DSourceComboUse &locSourceComboUse, const vector< int > &locBeamBunches, const vector< const DSourceCombo * > &locCombos, ComboingStage_t locComboingStage)
const DReactionStepVertexInfo * Get_StepVertexInfo(size_t locStepIndex) const
char string[256]
void Print_SourceComboInfo(const DSourceComboInfo *locComboInfo, unsigned char locNumTabs=0)
Definition: DSourceCombo.h:316
bool Cut_dEdxAndEOverP(const DChargedTrackHypothesis *locHypo)
vector< DSourceCombo * > dCreatedCombos
unordered_set< const DSourceComboInfo * > dComboInfosWithMassiveNeutrals
DCombosByReaction Build_ParticleCombos(const DReactionVertexInfo *locReactionVertexInfo)
vector< pair< DSourceComboUse, vector< const DSourceCombo * >>> DSourceCombosByUse_Small
Definition: DSourceCombo.h:35
void Combo_Vertically_AllDecays(const DSourceComboUse &locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, unsigned char locNumTabs)
vector< int > Get_ValidRFBunches(const JObject *locObject, signed char locVertexZBin) const
map< pair< const DSourceCombo *, signed char >, vector< int > > dValidRFBunches_ByCombo
bool Check_Skims(const DReaction *locReaction) const
const DTrackTimeBased * Get_TrackTimeBased(void) const
const DSourceCombo * Get_StepSourceCombo(const DReaction *locReaction, size_t locDesiredStepIndex, const DSourceCombo *locVertexPrimaryCombo, size_t locVertexPrimaryStepIndex=0) const
void Combo_Horizontally_AddParticles(const DSourceComboUse &locComboUseToCreate, const DSourceComboUse &locComboUseAllBut1, const pair< Particle_t, unsigned char > &locParticlePairToAdd, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, bool locExpandAllBut1Flag, unsigned char locNumTabs)
bool Get_HasMassiveNeutrals(const DSourceComboInfo *locComboInfo)
Definition: DSourceCombo.h:508
map< pair< const DReactionVertexInfo *, vector< signed char > >, DSourceComboUse > dSourceComboUseVertexZMap
unordered_map< const DReaction *, size_t > dRFBunchCutsByReaction
static signed char Get_VertexZIndex_OutOfRange(void)
Definition: DSourceCombo.h:80
unordered_map< const DSourceCombo *, DSourceCombosByBeamBunchByUse > dSourceCombosByBeamBunchByUse
void Create_SourceComboInfos(const DReactionVertexInfo *locReactionVertexInfo)
vector< const DKinematicData * > Get_BeamParticlesByRFBunch(int locRFBunch, unsigned int locPlusMinusBunchRange) const
bool Cut_InvariantMass_NoMassiveNeutrals(const DSourceCombo *locVertexCombo, Particle_t locDecayPID, Particle_t locTargetPIDToSubtract, const DVector3 &locVertex, signed char locVertexZBin, bool locAccuratePhotonsFlag)
static unsigned short int IsDetachedVertex(Particle_t p)
Definition: particleType.h:817
const DSourceCombo * Find_Combo_AtThisStep(const DSourceCombo *locSourceCombo, DSourceComboUse locUseToFind, size_t locDecayInstanceIndex) const
static char * ParticleType(Particle_t p)
Definition: particleType.h:142
map< bool, vector< const JObject * > > dTracksByCharge
map< Particle_t, map< DetectorSystem_t, vector< pair< double, double > > > > dEOverPValueMap
DetectorSystem_t
Definition: GlueX.h:15
int Calc_RFBunchShift(double locTimeToStepTo) const
const vector< const DSourceCombo * > & Get_CombosForComboing(const DSourceComboUse &locComboUse, ComboingStage_t locComboingStage, const vector< int > &locBeamBunches, const DSourceCombo *locChargedCombo_PresidingPrevious)
static signed char Get_VertexZIndex_Unknown(void)
Definition: DSourceCombo.h:82
map< Particle_t, map< DetectorSystem_t, pair< TF1 *, TF1 * > > > ddEdxCutMap
unordered_map< const DReactionVertexInfo *, size_t > dMaxRFBunchCuts
vector< pair< Particle_t, unsigned char > > Get_NumParticles(bool locEntireChainFlag=false) const
Definition: DSourceCombo.h:247
static int ParticleCharge(Particle_t p)
map< pair< const DSourceCombo *, const DReactionStepVertexInfo * >, const DSourceCombo * > dVertexPrimaryComboMap
const DParticleCombo * Build_ParticleCombo(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locFullCombo, const DKinematicData *locBeamParticle, int locRFBunchShift, DKinFitType locKinFitType)
Charge_t
DSourceComboUse Build_NewZDependentUse(const DReaction *locReaction, size_t locStepIndex, signed char locVertexZBin, const DSourceComboUse &locOrigUse, const unordered_map< size_t, DSourceComboUse > &locCreatedUseMap)
DSourceComboTimeHandler * dSourceComboTimeHandler
DResourcePool< DSourceCombo > dResourcePool_SourceCombo
bool Get_IsComboingZIndependent(const JObject *locObject, Particle_t locPID) const
static signed char Get_VertexZIndex_ZIndependent(void)
Definition: DSourceCombo.h:81
bool Check_IfMissingDecayProduct(const DReaction *locReaction, size_t locStepIndex)
Definition: DReaction.h:259
Definition: GlueX.h:17
DPhotonKinematicsByZBin Get_PhotonKinematics(void) const
Definition: GlueX.h:19
shared_ptr< const DTOFHitMatchParams > Get_TOFHitMatchParams(void) const
unordered_map< const DSourceComboInfo *, Charge_t > dComboInfoChargeContent
JApplication * japp
tuple< Particle_t, signed char, const DSourceComboInfo *, bool, Particle_t > DSourceComboUse
Definition: DSourceCombo.h:34
bool Cut_InvariantMass_AccuratePhotonKinematics(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionFullCombo, const DKinematicData *locBeamParticle, int locRFBunch)
vector< Particle_t > Get_FinalPIDs(int locStepIndex=-1, bool locIncludeMissingFlag=true, bool locIncludeDecayingFlag=true, Charge_t locCharge=d_AllCharges, bool locIncludeDuplicatesFlag=true) const
Definition: DReaction.cc:17
map< pair< const DReactionStepVertexInfo *, DSourceComboUse >, size_t > dSourceComboInfoStepMap
void Create_SourceCombos_Unknown(const DSourceComboUse &locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, unsigned char locNumTabs)
vector< vector< const DSourceCombo * > * > dCreatedComboVectors
const DESSkimData * dESSkimData
DParticleComboCreator * dParticleComboCreator
DSourceComboUse Get_ZIndependentUse(const DSourceComboUse &locZDependentUse)
const DReaction * Get_Reaction(void) const
Charge_t Get_ChargeContent(const DSourceComboInfo *locSourceComboInfo) const
unordered_map< signed char, DPhotonShowersByBeamBunch > Get_ShowersByBeamBunchByZBin(void) const
pair< bool, map< DSourceComboUse, unsigned char > > Get_FinalStateDecayingComboUses(const DReaction *locReaction, size_t locStepIndex, const map< size_t, DSourceComboUse > &locStepComboUseMap) const
Charge_t Get_ChargeContent(const DSourceComboInfo *locSourceComboInfo)
Definition: DSourceCombo.h:467
void Set_BeamParticles(const vector< const DBeamPhoton * > &locBeamParticles)
const vector< const JObject * > & Get_ShowersByBeamBunch(const vector< int > &locBeamBunches, DPhotonShowersByBeamBunch &locShowersByBunch, signed char locVertexZBin)
void Combo_WithBeam(const vector< const DReaction * > &locReactions, const DReactionVertexInfo *locReactionVertexInfo, const DSourceComboUse &locReactionFullComboUse, const DSourceCombo *locReactionFullCombo, int locRFBunch, DCombosByReaction &locOutputComboMap)
void Vote_OldMethod(const DSourceCombo *locReactionFullCombo, vector< int > &locValidRFBunches)
Definition: GlueX.h:20
void Create_Combo_OneDecay(const DSourceComboUse &locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, unsigned char locNumTabs)
const DSourceCombo * Get_ChargedCombo_WithNow(const DSourceCombo *locChargedCombo_Presiding, const DSourceComboInfo *locToCreateComboInfo, ComboingStage_t locComboingStage) const
Definition: GlueX.h:18
map< const DReaction *, TH1 * > dNumEventsSurvivedStageMap
const DSourceCombo * Get_NextChargedCombo(const DSourceCombo *locChargedCombo_Presiding, const DSourceComboUse &locNextComboUse, ComboingStage_t locComboingStage, bool locGetPresidingFlag, size_t locInstance) const
void Combo_Vertically_NDecays(const DSourceComboUse &locComboUseToCreate, const DSourceComboUse &locNMinus1ComboUse, const DSourceComboUse &locSourceComboDecayUse, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, unsigned char locNumTabs)
bool Get_HasPhotons(const DSourceComboInfo *locSourceComboInfo) const
vector< const DSourceComboInfo * > dSourceComboInfos
signed char Get_VertexZBin(bool locIsProductionVertex, const DSourceCombo *locReactionCombo, const DSourceCombo *locPrimaryVertexCombo, const DKinematicData *locBeamParticle, bool locIsCombo2ndVertex) const
map< DSourceComboUse, size_t > dNumMixedCombosMap_Charged
DGeometry * GetDGeometry(unsigned int run_number)
string Get_ReactionName(void) const
Definition: DReaction.h:75
Definition: GlueX.h:22
DSourceComboUse Create_ZDependentSourceComboUses(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionChargedCombo)
DSourceCombosByBeamBunchByUse & Get_SourceCombosByBeamBunchByUse(Charge_t locChargeContent_SearchForUse, const DSourceCombo *locChargedCombo=nullptr)
bool Cut_Timing_MissingMassVertices(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionFullCombo, const DKinematicData *locBeamParticle, int locRFBunch)
bool Get_ExpandAllBut1Flag(ComboingStage_t locComboingStage, const DSourceComboUse &locAllBut1ComboUse, Charge_t locToAddChargeContent)
bool Cut_dEdx(Particle_t locPID, DetectorSystem_t locSystem, double locP, double locdEdx)
map< Particle_t, vector< const JObject * > > dTracksByPID
void Calc_VertexTimeOffsets_WithBeam(const DReactionVertexInfo *locReactionVertexInfo, const DSourceComboUse &locReactionFullComboUse, const DSourceCombo *locReactionFullCombo, const DKinematicData *locBeamParticle)
void Check_ForDuplicates(const vector< const DSourceCombo * > &locCombos) const
DSourceComboUse Make_ComboUse(Particle_t locInitPID, const map< Particle_t, unsigned char > &locNumParticles, const map< DSourceComboUse, unsigned char > &locFurtherDecays, bool locMissingDecayProductFlag, Particle_t locTargetToInclude)
unordered_map< const DSourceCombo *, DSourceCombosByUse_Large > dMixedCombosByUseByChargedCombo
void Create_SourceCombos(const DSourceComboUse &locComboUseToCreate, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, unsigned char locNumTabs)
h_means GetYaxis() -> SetTitleOffset(1.8)
map< Particle_t, map< DetectorSystem_t, vector< pair< double, double > > > > ddEdxValueMap
void Register_ValidRFBunches(const DSourceComboUse &locSourceComboUse, const DSourceCombo *locSourceCombo, const vector< int > &locRFBunches, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding)
TH2D * locHist
vector< signed char > Get_VertexZBins(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionCombo, const DKinematicData *locBeamParticle, bool locComboIsFullyCharged) const
h1eff_eff GetXaxis() -> SetTitleSize(0.05)
void Combo_Horizontally_AddCombo(const DSourceComboUse &locComboUseToCreate, const DSourceComboUse &locAllBut1ComboUse, const DSourceComboUse &locSourceComboUseToAdd, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, bool locExpandAllBut1Flag, unsigned char locNumTabs)
DSourceCombosByUse_Large & Get_CombosSoFar(ComboingStage_t locComboingStage, Charge_t locChargeContent_SearchForUse, const DSourceCombo *locChargedCombo=nullptr)
vector< const JObject * > Get_SourceParticles(const vector< pair< Particle_t, const JObject * >> &locSourceParticles, Particle_t locPID=Unknown)
Definition: DSourceCombo.h:369
map< Particle_t, unsigned char > Build_ParticleMap(const DReaction *locReaction, size_t locStepIndex, Charge_t locCharge) const
void Set_Vertex(const DVertex *locVertex)
bool Select_RFBunches_Charged(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locChargedCombo, vector< int > &locValidRFBunches)
map< DSourceComboUse, DSourceComboUse > dZDependentUseToIndependentMap
shared_ptr< const DSCHitMatchParams > Get_SCHitMatchParams(void) const
vector< pair< DSourceComboUse, unsigned char > > Get_FurtherDecays(void) const
Definition: DSourceCombo.h:76
const char * SystemName(DetectorSystem_t sys)
Definition: GlueX.h:38
bool Select_RFBunches_AllVerticesUnknown(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionFullCombo, Charge_t locCharge, vector< int > &locValidRFBunches)
unordered_map< const DReactionStepVertexInfo *, DSourceComboUse > dSourceComboUseReactionMap
size_t Get_ParticleInstanceIndex(const DReactionStep *locStep, size_t locParticleIndex)
Definition: DReaction.cc:118
static double ParticleMass(Particle_t p)
vector< const DReactionStepVertexInfo * > Get_StepVertexInfos(void) const
vector< const DSourceCombo * > * Get_SourceComboVectorResource(void)
bool Check_Reactions(vector< const DReaction * > &locReactions)
shared_ptr< const DFCALShowerMatchParams > Get_FCALShowerMatchParams(void) const
const DSourceCombo * Get_VertexPrimaryCombo(const DSourceCombo *locReactionCombo, const DReactionStepVertexInfo *locStepVertexInfo) const
map< vector< const JObject * >, const DSourceCombo * > dNPhotonsToComboMap
map< Particle_t, map< DetectorSystem_t, TH2 * > > dHistMap_EOverP
void Print_SourceComboUse(const DSourceComboUse &locComboUse, unsigned char locNumTabs=0, bool locIgnoreTabs=false)
Definition: DSourceCombo.h:337
void Reset_NewEvent(JEventLoop *locEventLoop)
void Set_SourceComboVertexer(const DSourceComboVertexer *locSourceComboVertexer)
DResourcePool< vector< const DSourceCombo * > > dResourcePool_SourceComboVector
void Build_ComboResumeIndices(const DSourceComboUse &locSourceComboUse, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding)
vector< const DReactionStepVertexInfo * > Get_StepVertexInfos_ReverseOrderByStep(const DReactionVertexInfo *locReactionVertexInfo)
void Copy_ZIndependentMixedResults(const DSourceComboUse &locComboUseToCreate, const DSourceCombo *locChargedCombo_Presiding)
map< DSourceComboUse, size_t > dNumMixedCombosMap_Mixed
const DVector3 & momentum(void) const
unordered_map< signed char, DPhotonShowersByBeamBunch > dShowersByBeamBunchByZBin
void Combo_Vertically_NParticles(const DSourceComboUse &locComboUseToCreate, const DSourceComboUse &locNMinus1ComboUse, ComboingStage_t locComboingStage, unsigned char locNumTabs)
bool Cut_EOverP(Particle_t locPID, DetectorSystem_t locSystem, double locP, double locEOverP)
size_t Get_ResumeAtIndex_Particles(Particle_t locPID, const JObject *locPreviousObject, vector< int > locBeamBunches, signed char locVertexZBin) const
size_t Get_NumObjectsAllThreads(void) const
bool Check_NumParticles(const DReaction *locReaction)
map< const DReaction *, map< DConstructionStage, size_t > > dNumCombosSurvivedStageTracker
const DReactionStep * Get_ReactionStep(size_t locStepIndex) const
Definition: DReaction.h:84
void Set_ControlParams(size_t locGetBatchSize, size_t locNumToAllocateAtOnce, size_t locMaxLocalPoolSize)
int Get_DecayStepIndex(const DReaction *locReaction, size_t locStepIndex, size_t locParticleIndex)
Definition: DReaction.cc:135
map< Particle_t, map< DetectorSystem_t, string > > dEOverPCuts_TF1FunctionStrings
shared_ptr< const DBCALShowerMatchParams > Get_BCALShowerMatchParams(void) const
map< Particle_t, map< DetectorSystem_t, pair< string, string > > > ddEdxCuts_TF1FunctionStrings
bool Get_IsVertexKnown_NoBeam(bool locIsProductionVertex, const DSourceCombo *locVertexCombo, bool locIsCombo2ndVertex) const
DSourceCombosByUse_Small Get_FurtherDecayCombos(void) const
Definition: DSourceCombo.h:176
void Combo_Horizontally_AddDecay(const DSourceComboUse &locComboUseToCreate, const DSourceComboUse &locComboUseAllBut1, const DSourceComboUse &locComboUseToAdd, ComboingStage_t locComboingStage, const DSourceCombo *locChargedCombo_Presiding, bool locExpandAllBut1Flag, unsigned char locNumTabs)
void Set_SourceComboTimeHandler(const DSourceComboTimeHandler *locSourceComboTimeHandler)
vector< const DReactionStep * > Get_ReactionSteps(void) const
Definition: DReaction.h:85
void Calc_VertexTimeOffsets_WithCharged(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionChargedCombo)
void Set_PhotonKinematics(const DPhotonKinematicsByZBin &locPhotonKinematics)
vector< size_t > Get_StepIndices(void) const
DetectorSystem_t t1_detector(void) const
bool Get_PromoteFlag(ComboingStage_t locComboingStage, Particle_t locDecayPID_UseToCheck, const DSourceComboInfo *locComboInfo_UseToCreate, const DSourceComboInfo *locComboInfo_UseToCheck, DSourceComboUse &locNonNeutralUse) const
void Calc_VertexTimeOffsets_WithPhotons(const DReactionVertexInfo *locReactionVertexInfo, const DSourceCombo *locReactionChargedCombo, const DSourceCombo *locReactionFullCombo)