fixed processes with multiple non-interfering cross sections

madgraph5 · Nov 6, 2024 · 63bbbfc · 63bbbfc
1 parent 910bfce
commit 63bbbfc
Show file tree

Hide file tree

Showing 3 changed files with 29 additions and 8 deletions.
diff --git a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/trex.py b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/trex.py
@@ -656,8 +656,9 @@ def do_quit(self, line):
             for line in self.banner['init'].split('\n'):
                 split = line.split()
                 if len(split) == 4:
-                    cross, error = float(split[0]), float(split[1])
-
+                    cross += float(split[0])
+                    error += float(split[1])**2
+            error = error**0.5
         if not self.multicore == 'create':
             # No print of results for the multicore mode for the one printed on screen
             if 'orig' not in self.all_cross_section:

diff --git a/tools/REX/rwgt_runner.cc b/tools/REX/rwgt_runner.cc
@@ -21,6 +21,16 @@
 namespace %(process_namespace)s{
 //namespace dummy{
 
+    std::vector<std::vector<std::string_view>> getInitPrts(){
+        static std::vector<std::vector<std::string_view>> initPrts = {%(init_prt_ids)s};
+        return initPrts;
+    }
+
+    std::vector<std::vector<std::string_view>> getFinPrts(){
+        static std::vector<std::vector<std::string_view>> finPrts = {%(fin_prt_ids)s};
+        return finPrts;
+    }
+
     std::shared_ptr<std::vector<FORTRANFPTYPE>> amp( int& nEvt, int& nPar, int& nMom, std::vector<FORTRANFPTYPE>& momenta, std::vector<FORTRANFPTYPE>& alphaS, std::vector<FORTRANFPTYPE>& rndHel, std::vector<FORTRANFPTYPE>& rndCol, std::vector<int>& selHel, std::vector<int>& selCol, unsigned int& chanId, bool& goodHel ){
         CppObjectInFortran *bridgeInst;
         auto evalScatAmps = std::make_shared<std::vector<FORTRANFPTYPE>>( nEvt );
@@ -45,8 +55,8 @@ namespace %(process_namespace)s{
     }
 
     std::shared_ptr<std::vector<size_t>> procSort( std::string_view status, std::vector<std::string_view> arguments, size_t index ){
-        std::vector<std::vector<std::string_view>> initPrts = {%(init_prt_ids)s};
-        std::vector<std::vector<std::string_view>> finPrts = {%(fin_prt_ids)s};
+        std::vector<std::vector<std::string_view>> initPrts = getInitPrts();
+        std::vector<std::vector<std::string_view>> finPrts = getFinPrts();
         std::shared_ptr<std::vector<size_t>> refOrder;
     if( index == REX::npos ){
 	if( status == "-1" ){
@@ -81,11 +91,13 @@ namespace %(process_namespace)s{
 
     bool checkProc( REX::event& process, std::vector<std::string>& relStats ){
         size_t no_evts = %(no_events)s;
+        auto finPrts = getFinPrts();
         for( size_t k = 0 ; k < no_evts ; ++k ){
             REX::statSort locSort = [ind = k](std::string_view status, std::vector<std::string_view> arguments){
                 return procSort( status, arguments, ind );
             };
             auto order = process.getProcOrder( locSort );
+            if( order.at("1").size() != finPrts[k].size() ){ continue; }
             for( size_t j = 0 ; j < relStats.size() ; ++j ){
                 auto currPts = order.at( relStats[j] );
                 if( std::find(currPts.begin(), currPts.end(), REX::npos) != currPts.end() ){ break; }

diff --git a/tools/REX/teawREX.cc b/tools/REX/teawREX.cc
@@ -507,11 +507,14 @@ namespace REX::teaw
             if( this->ampNorm != 0.0 ){ return; }
             auto xSecLines = this->lheFile->getInit()->getLines();
             if( xSecLines.size() > 1 ){
-                std::cout << "\n\033[1;33mWarning: Multiple cross-section lines found in LHE file.\nteawREX only supports single (inclusive) process reweighting.\nWill proceed assuming all events belong to first process type.\033[0m\n";
+                std::cout << "\n\033[1;33mWarning: Multiple cross-section lines found in LHE file.\nAssuming total cross section given by sum of all cross sections.\033[0m\n";
             }
             if( xSecLines.size() == 0 )
                 throw std::runtime_error( "No cross-section information found in LHE file." );
-            auto xSec = std::stod(std::string(xSecLines[0]->xsecup));
+            double xSec = 0.0;
+            for( size_t k = 0 ; k < xSecLines.size() ; ++k ){
+                xSec += std::stod(std::string(xSecLines[k]->xsecup));
+            }
             double div = 0.0;
             bool sameWeight = true;
             for( size_t k = 1 ; k < this->flatWgts->size() - 1 ; k += size_t(flatWgts->size()/21) ){
@@ -658,8 +661,13 @@ namespace REX::teaw
             double invN = 1. / double(reWgts->at(0)->size());
             double sqrtInvN = std::sqrt( invN );
             auto xSecLines = this->lheFile->getInit()->getLines();
-            double xSec = std::stod(std::string(xSecLines[0]->xsecup));
-            double xErr = std::stod(std::string(xSecLines[0]->xerrup));
+            double xSec = 0.0;
+            double xErr = 0.0;
+            for( size_t k = 0 ; k < xSecLines.size() ; ++k ){
+                xSec += std::stod(std::string(xSecLines[k]->xsecup));
+                xErr += std::pow(std::stod(std::string(xSecLines[k]->xerrup)),2);
+            }
+            xErr = std::sqrt( xErr );
             for( size_t k = 0 ; k < reWgts->size() ; ++k ){
                 double xSecCurr = normXSecs->at(k);
                 auto locWgts = reWgts->at(k);