Changeset 4006893 in git

Timestamp:

Jun 1, 2021, 8:10:27 PM (3 years ago)

Author:

Pavel Demin <pavel.demin@…>

Branches:

master

Children:

fc6bce2

Parents:

06630ef

Message:

fix order of decay products in DelphesHepMC3Reader

Location:

classes

Files:

• DelphesHepMC3Reader.cc (modified) (16 diffs)
• DelphesHepMC3Reader.h (modified) (4 diffs)

classes/DelphesHepMC3Reader.cc

@@ -56 +56 @@
 DelphesHepMC3Reader::DelphesHepMC3Reader() :
   fInputFile(0), fBuffer(0), fPDG(0),
-  fVertexCounter(-1), fParticleCounter(-1)
+  fVertexCounter(-2), fParticleCounter(-1)
 {
   fBuffer = new char[kBufferSize];
@@ -81 +81 @@
 void DelphesHepMC3Reader::Clear()
 {
-  fWeight.clear();
+  fWeights.clear();
   fMomentumCoefficient = 1.0;
   fPositionCoefficient = 1.0;
-  fVertexCounter = -1;
+  fVertexCounter = -2;
   fParticleCounter = -1;
-  fVertexMap.clear();
+  fVertices.clear();
+  fParticles.clear();
+  fInVertexMap.clear();
+  fOutVertexMap.clear();
+  fMotherMap.clear();
   fDaughterMap.clear();
 }
@@ -94 +98 @@
 bool DelphesHepMC3Reader::EventReady()
 {
-  return (fParticleCounter == 0);
+  return (fVertexCounter == -1) && (fParticleCounter == 0);
 }
 
@@ -119 +123 @@
     Clear();
 
-    fX = 0.0;
-    fY = 0.0;
-    fZ = 0.0;
-    fT = 0.0;
-
     rc = bufferStream.ReadInt(fEventNumber)
       && bufferStream.ReadInt(fVertexCounter)
@@ -168 +167 @@
     while(bufferStream.ReadDbl(weight))
     {
-      fWeight.push_back(weight);
+      fWeights.push_back(weight);
     }
   }
@@ -257 +256 @@
   else if(key == 'V')
   {
+    fParticles.clear();
+
     fX = 0.0;
     fY = 0.0;
@@ -262 +263 @@
     fT = 0.0;
 
-    fM1 = 0;
-    fM2 = 0;
-
     rc = bufferStream.ReadInt(fVertexCode)
       && bufferStream.ReadInt(fVertexStatus);
@@ -286 +284 @@
     while(bufferStream.ReadInt(code))
     {
-      if(code < fM1 || fM1 == 0) fM1 = code;
-      if(code > fM2) fM2 = code;
-      fVertexMap[code] = fVertexCode;
-    }
-
-    if(fM1 == fM2) fM2 = 0;
+      fParticles.push_back(code);
+      bufferStream.FindChr(',');
+    }
 
     if(bufferStream.FindChr('@'))
@@ -307 +302 @@
       }
     }
+
+    AnalyzeVertex(factory, fVertexCode);
   }
   else if(key == 'P' && fParticleCounter > 0)
@@ -329 +326 @@
     }
 
-    itDaughterMap = fDaughterMap.find(fOutVertexCode);
-    if(itDaughterMap == fDaughterMap.end())
-    {
-      fDaughterMap[fOutVertexCode] = make_pair(fParticleCode, fParticleCode);
-    }
-    else
-    {
-      itDaughterMap->second.second = fParticleCode;
-    }
-
     AnalyzeParticle(factory, allParticleOutputArray,
       stableParticleOutputArray, partonOutputArray);
@@ -363 +350 @@
   element->ProcessID = fProcessID;
   element->MPI = fMPI;
-  element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
+  element->Weight = fWeights.size() > 0 ? fWeights[0] : 1.0;
   element->CrossSection = fCrossSection;
   element->CrossSectionError = fCrossSectionError;
@@ -389 +376 @@
   vector<double>::const_iterator itWeight;
 
-  for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
+  for(itWeight = fWeights.begin(); itWeight != fWeights.end(); ++itWeight)
   {
     element = static_cast<Weight *>(branch->NewEntry());
 
     element->Weight = *itWeight;
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate)
+{
+  int index;
+  TLorentzVector *position;
+  TObjArray *array;
+  vector<int>::iterator itParticle;
+  map<int, int>::iterator itVertexMap;
+
+  itVertexMap = fOutVertexMap.find(code);
+  if(itVertexMap == fOutVertexMap.end())
+  {
+    --fVertexCounter;
+
+    index = fVertices.size();
+    fOutVertexMap[code] = index;
+    if(candidate && code > 0) fInVertexMap[code] = index;
+
+    position = factory->New<TLorentzVector>();
+    array = factory->NewArray();
+    position->SetXYZT(0.0, 0.0, 0.0, 0.0);
+    array->Clear();
+    fVertices.push_back(make_pair(position, array));
+  }
+  else
+  {
+    index = itVertexMap->second;
+    position = fVertices[index].first;
+    array = fVertices[index].second;
+  }
+
+  if(candidate)
+  {
+    array->Add(candidate);
+  }
+  else
+  {
+    position->SetXYZT(fX, fY, fZ, fT);
+    for(itParticle = fParticles.begin(); itParticle != fParticles.end(); ++itParticle)
+    {
+      fInVertexMap[*itParticle] = index;
+    }
   }
 }
@@ -425 +458 @@
   }
 
-  candidate->Position.SetXYZT(fX, fY, fZ, fT);
-  if(fPositionCoefficient != 1.0)
-  {
-    candidate->Position *= fPositionCoefficient;
-  }
-
-  candidate->D1 = -1;
-  candidate->D2 = -1;
-
-  if(fOutVertexCode < 0)
-  {
-    candidate->M1 = fM1 - 1;
-    candidate->M2 = fM2 - 1;
-  }
-  else
-  {
-    candidate->M1 = fOutVertexCode - 1;
-    candidate->M2 = -1;
-  }
-
-  allParticleOutputArray->Add(candidate);
+  candidate->D1 = fParticleCode;
+
+  AnalyzeVertex(factory, fOutVertexCode, candidate);
 
   if(!pdgParticle) return;
@@ -463 +478 @@
 void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray)
 {
+  TLorentzVector *position;
+  TObjArray *array;
   Candidate *candidate;
   map<int, int >::iterator itVertexMap;
+  map<int, pair<int, int> >::iterator itMotherMap;
   map<int, pair<int, int> >::iterator itDaughterMap;
-  int i, index;
+  int i, j, code, counter;
+
+  counter = 0;
+  for(i = 0; i < fVertices.size(); ++i)
+  {
+    position = fVertices[i].first;
+    array = fVertices[i].second;
+
+    for(j = 0; j < array->GetEntriesFast(); ++j)
+    {
+      candidate = static_cast<Candidate *>(array->At(j));
+
+      candidate->Position = *position;
+      if(fPositionCoefficient != 1.0)
+      {
+        candidate->Position *= fPositionCoefficient;
+      }
+
+      candidate->M1 = i;
+
+      itDaughterMap = fDaughterMap.find(i);
+      if(itDaughterMap == fDaughterMap.end())
+      {
+        fDaughterMap[i] = make_pair(counter, counter);
+      }
+      else
+      {
+        itDaughterMap->second.second = counter;
+      }
+
+      code = candidate->D1;
+
+      itVertexMap = fInVertexMap.find(code);
+      if(itVertexMap == fInVertexMap.end())
+      {
+        candidate->D1 = -1;
+      }
+      else
+      {
+        code = itVertexMap->second;
+
+        candidate->D1 = code;
+
+        itMotherMap = fMotherMap.find(code);
+        if(itMotherMap == fMotherMap.end())
+        {
+          fMotherMap[code] = make_pair(counter, -1);
+        }
+        else
+        {
+          itMotherMap->second.second = counter;
+        }
+      }
+
+      allParticleOutputArray->Add(candidate);
+
+      ++counter;
+    }
+  }
 
   for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
@@ -472 +548 @@
     candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
 
-    index = i + 1;
-
-    itVertexMap = fVertexMap.find(index);
-    if(itVertexMap != fVertexMap.end())
-    {
-      index = itVertexMap->second;
-    }
-
-    itDaughterMap = fDaughterMap.find(index);
-    if(itDaughterMap == fDaughterMap.end())
+    itMotherMap = fMotherMap.find(candidate->M1);
+    if(itMotherMap == fMotherMap.end())
+    {
+      candidate->M1 = -1;
+      candidate->M2 = -1;
+    }
+    else
+    {
+      candidate->M1 = itMotherMap->second.first;
+      candidate->M2 = itMotherMap->second.second;
+    }
+
+    if(candidate->D1 < 0)
     {
       candidate->D1 = -1;
@@ -488 +567 @@
     else
     {
-      candidate->D1 = itDaughterMap->second.first - 1;
-      candidate->D2 = itDaughterMap->second.second - 1;
-    }
-  }
-}
-
-//---------------------------------------------------------------------------
+      itDaughterMap = fDaughterMap.find(candidate->D1);
+      if(itDaughterMap == fDaughterMap.end())
+      {
+        candidate->D1 = -1;
+        candidate->D2 = -1;
+      }
+      else
+      {
+        candidate->D1 = itDaughterMap->second.first;
+        candidate->D2 = itDaughterMap->second.second;
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------

classes/DelphesHepMC3Reader.h

@@ -36 +36 @@
 class TStopwatch;
 class TDatabasePDG;
+class TLorentzVector;
 class ExRootTreeBranch;
 class DelphesFactory;
+class Candidate;
 
 class DelphesHepMC3Reader
@@ -61 +63 @@
 
 private:
+  void AnalyzeVertex(DelphesFactory *factory, int code, Candidate *candidate = 0);
+
   void AnalyzeParticle(DelphesFactory *factory,
     TObjArray *allParticleOutputArray,
@@ -79 +83 @@
   double fMomentumCoefficient, fPositionCoefficient;
 
-  std::vector<double> fWeight;
+  std::vector<double> fWeights;
 
   double fCrossSection, fCrossSectionError;
@@ -92 +96 @@
   double fPx, fPy, fPz, fE, fMass;
 
-  int fM1, fM2;
+  std::vector<std::pair<TLorentzVector *, TObjArray *> > fVertices;
+  std::vector<int> fParticles;
 
-  std::map<int, int> fVertexMap;
+  std::map<int, int> fInVertexMap;
+  std::map<int, int> fOutVertexMap;
+
+  std::map<int, std::pair<int, int> > fMotherMap;
   std::map<int, std::pair<int, int> > fDaughterMap;
 };