Merge branch '21.6-sp-seup-fix' into '21.6'

Fix for  Pythia8 UserHook to enable decay to SUEP

See merge request atlas/athena!35491

Generators/Pythia8_i/src/UserHooks/DecayToSUEP.cxx

@@ -128,7 +128,29 @@ namespace Pythia8{
 
 	// Generate the shower, output are 4 vectors in the rest frame of the shower
 	higgs4mom=process[ii].p();
-	suep_shower4momenta=suep_shower.generate_shower();
+	int nShowerAttempts=0;
+	do {
+	  try {
+	    nShowerAttempts++;
+	    suep_shower4momenta=suep_shower.generate_shower();
+	    if( suep_shower4momenta.size()<3){
+	      //Failed to balance energy or less than 3 particles in the shower
+	      //Try again until nShowerAttempts >= 3
+	    } else {
+	      //All ok!
+	      nShowerAttempts = -1; //exit condition
+	    }
+	  } catch (std::exception &e) {
+	    //Failed to generate the shower!
+	    //Can happen in some rare circumstances, 
+	    //Try again until nShowerAttempts >= 3
+	  }
+	} while ((nShowerAttempts > 0) && (nShowerAttempts < 3));
+	if (nShowerAttempts >= 3) {
+	  //Something is seriously wrong then, print warning and skip to next event
+	  std::cout << "[SUEP] WARNING: Something went wrong in generating the shower. Skipping the event." << std::endl;
+	  return true; //veto the event!
+	}
 
 	// Loop over hidden sector mesons and append to the event	
 	for (unsigned j = 0; j < suep_shower4momenta.size(); ++j){

Generators/Pythia8_i/src/UserHooks/suep_shower.cxx

@@ -146,11 +146,30 @@ vector< Vec4 > Suep_shower::generate_shower(){
       event[n][i] += correction;
     } 
   }
+
+  //Shield against an exception in the calculation of "p_scale" further down. If it fails, abort the event.
+  if(Suep_shower::reballance_func(2.0,event)<0.0){
+    // failed to balance energy. 
+    event.clear();
+    return event;
+  } 
   
   // finally, ballance the total energy, without destroying momentum conservation
   tolerance tol = 0.00001;
   double p_scale;
-  p_scale = (bisect(boost::bind(&Suep_shower::reballance_func, this, _1, event),0.0,2.0, tol)).first;
+  try {
+    p_scale = (bisect(boost::bind(&Suep_shower::reballance_func, this, _1, event),0.0,2.0, tol)).first;
+  } catch (std::exception &e) {
+    //in some rare circumstances, this balancing might fail
+    std::cout << "[SUEP_SHOWER] WARNING: Failed to rebalance the following event; Printing for debug and throw exception" << std::endl;
+    std::cout << e.what() << std::endl;
+    std::cout << "N. Particle, px, py, pz, E" << std::endl;
+    for (size_t jj=0; jj < event.size(); jj++) {
+      //print out event
+      std::cout << jj << ": " << event[jj].px() << ", " << event[jj].py() << ", " << event[jj].pz() << ", " << event[jj].e() << std::endl;
+    }
+    throw e;
+  }
   
   for(int n=0;n<len;n++){
     event[n].px(p_scale*event[n].px());