Merge pull request #19 from vlvovch/devel

Fixed bug in CE Event generator where electric charge was mixed up with charm, resulting in small violation of exact charge conservation
vlvovch · Apr 12, 2019 · 1e91614 · 1e91614
2 parents 155fa46 + 9f992ac
commit 1e91614
Show file tree

Hide file tree

Showing 3 changed files with 62 additions and 2 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -8,7 +8,7 @@ project (ThermalFIST)
 # The version number.
 set (ThermalFIST_VERSION_MAJOR 1)
 set (ThermalFIST_VERSION_MINOR 1)
-set (ThermalFIST_VERSION_DEVEL 3)
+set (ThermalFIST_VERSION_DEVEL 4)
 
 # configure a header file to pass some of the CMake settings
 # to the source code

diff --git a/include/HRGEventGenerator/SimpleParticle.h b/include/HRGEventGenerator/SimpleParticle.h
@@ -51,6 +51,12 @@ namespace thermalfist {
     double GetY() const {
       return 0.5*log((p0 + pz) / (p0 - pz));
     }
+
+    /// The longitudinal pseudorapidity
+    double GetEta() const {
+      //return atanh(pz/GetP());
+      return 0.5*log((GetP() + pz) / (GetP() - pz));
+    }
   };
 
 } // namespace thermalfist

diff --git a/src/library/HRGEventGenerator/EventGeneratorBase.cpp b/src/library/HRGEventGenerator/EventGeneratorBase.cpp
@@ -621,6 +621,18 @@ namespace thermalfist {
         if (tind >= static_cast<int>(fAntiStrangeMesonsc.size())) tind = fAntiStrangeMesonsc.size() - 1;
         totals[fAntiStrangeMesonsc[tind].second]++;
       }
+
+      // Cross-check that all resulting strangeness is zero
+      int finS = 0;
+      for (size_t i = 0; i < totals.size(); ++i) {
+        finS += totals[i] * m_THM->TPS()->Particles()[i].Strangeness();
+      }
+
+      if (finS != 0) {
+        printf("**ERROR** EventGeneratorBase::GenerateTotalsSCESubVolume(): Generated strangeness is non-zero!");
+        exit(1);
+      }
+
       return totals;
     }
     return totals;
@@ -799,6 +811,18 @@ namespace thermalfist {
       if (tind >= static_cast<int>(fAntiCharmAllc.size())) tind = fAntiCharmAllc.size() - 1;
       totals[fAntiCharmAllc[tind].second]++;
     }
+
+    // Cross-check that all resulting strangeness is zero
+    int finC = 0;
+    for (size_t i = 0; i < totals.size(); ++i) {
+      finC += totals[i] * m_THM->TPS()->Particles()[i].Charm();
+    }
+
+    if (finC != m_THM->Parameters().C) {
+      printf("**ERROR** EventGeneratorBase::GenerateTotalsCCESubVolume(): Generated charm is non-zero!");
+      exit(1);
+    }
+
     return totals;
   }
 
@@ -887,7 +911,8 @@ namespace thermalfist {
         if (tind >= static_cast<int>(fAntiBaryonsc.size())) tind = fAntiBaryonsc.size() - 1;
         totals[fAntiBaryonsc[tind].second]++;
         netS += m_THM->TPS()->Particles()[fAntiBaryonsc[tind].second].Strangeness();
-        netQ += m_THM->TPS()->Particles()[fAntiBaryonsc[tind].second].Charm();
+        netQ += m_THM->TPS()->Particles()[fAntiBaryonsc[tind].second].ElectricCharge();
+        netC += m_THM->TPS()->Particles()[fAntiBaryonsc[tind].second].Charm();
       }
 
       // Total numbers of (anti)strange mesons
@@ -974,6 +999,35 @@ namespace thermalfist {
         }
       }
 
+      // Cross-check that all resulting charges are OK
+      int finB = 0, finQ = 0, finS = 0, finC = 0;
+      for (size_t i = 0; i < totals.size(); ++i) {
+        finB += totals[i] * m_THM->TPS()->Particles()[i].BaryonCharge();
+        finQ += totals[i] * m_THM->TPS()->Particles()[i].ElectricCharge();
+        finS += totals[i] * m_THM->TPS()->Particles()[i].Strangeness();
+        finC += totals[i] * m_THM->TPS()->Particles()[i].Charm();
+      }
+
+      if (finB != m_THM->Parameters().B) {
+        printf("**ERROR** EventGeneratorBase::GenerateTotalsCE(): Generated baryon number does not match the input!");
+        exit(1);
+      }
+
+      if (finQ != m_THM->Parameters().Q) {
+        printf("**ERROR** EventGeneratorBase::GenerateTotalsCE(): Generated electric charge does not match the input!");
+        exit(1);
+      }
+
+      if (finS != m_THM->Parameters().S) {
+        printf("**ERROR** EventGeneratorBase::GenerateTotalsCE(): Generated strangeness does not match the input!");
+        exit(1);
+      }
+
+      if (finC != m_THM->Parameters().C) {
+        printf("**ERROR** EventGeneratorBase::GenerateTotalsCE(): Generated charm does not match the input!");
+        exit(1);
+      }
+
       return totals;
     }
     return totals;