Changes in classes/DelphesHepMC3Reader.cc [b4786d3:a26a130] in git

File:

• classes/DelphesHepMC3Reader.cc (modified) (14 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);
+    AnalyzeParticle(factory);
   }
 
   if(EventReady())
   {
-    FinalizeParticles(allParticleOutputArray);
+    FinalizeParticles(allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
   }
 
@@ -363 +349 @@
   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;
@@ -387 +373 @@
 {
   Weight *element;
-  vector<double>::const_iterator itWeight;
-
-  for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
+  vector<double>::const_iterator itWeights;
+
+  for(itWeights = fWeights.begin(); itWeights != fWeights.end(); ++itWeights)
   {
     element = static_cast<Weight *>(branch->NewEntry());
 
-    element->Weight = *itWeight;
-  }
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory,
-  TObjArray *allParticleOutputArray,
+    element->Weight = *itWeights;
+  }
+}
+
+//---------------------------------------------------------------------------
+
+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);
+    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;
+    }
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory)
+{
+  Candidate *candidate;
+
+  candidate = factory->NewCandidate();
+
+  candidate->PID = fPID;
+
+  candidate->Status = fParticleStatus;
+
+  candidate->Mass = fMass;
+
+  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
+
+  candidate->D1 = fParticleCode;
+
+  AnalyzeVertex(factory, fOutVertexCode, candidate);
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray,
   TObjArray *stableParticleOutputArray,
   TObjArray *partonOutputArray)
 {
+  TLorentzVector *position;
+  TObjArray *array;
   Candidate *candidate;
   TParticlePDG *pdgParticle;
   int pdgCode;
-
-  candidate = factory->NewCandidate();
-
-  candidate->PID = fPID;
-  pdgCode = TMath::Abs(candidate->PID);
-
-  candidate->Status = fParticleStatus;
-
-  pdgParticle = fPDG->GetParticle(fPID);
-  candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
-  candidate->Mass = fMass;
-
-  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
-  if(fMomentumCoefficient != 1.0)
-  {
-    candidate->Momentum *= fMomentumCoefficient;
-  }
-
-  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);
-
-  if(!pdgParticle) return;
-
-  if(fParticleStatus == 1)
-  {
-    stableParticleOutputArray->Add(candidate);
-  }
-  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
-  {
-    partonOutputArray->Add(candidate);
-  }
-}
-
-//---------------------------------------------------------------------------
-
-void DelphesHepMC3Reader::FinalizeParticles(TObjArray *allParticleOutputArray)
-{
-  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;
+      }
+
+      if(fMomentumCoefficient != 1.0)
+      {
+        candidate->Momentum *= fMomentumCoefficient;
+      }
+
+      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;
+
+      pdgParticle = fPDG->GetParticle(candidate->PID);
+
+      candidate->Charge = pdgParticle ? int(pdgParticle->Charge() / 3.0) : -999;
+
+      if(!pdgParticle) continue;
+
+      pdgCode = TMath::Abs(candidate->PID);
+
+      if(candidate->Status == 1)
+      {
+        stableParticleOutputArray->Add(candidate);
+      }
+      else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
+      {
+        partonOutputArray->Add(candidate);
+      }
+    }
+  }
 
   for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
@@ -472 +549 @@
     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 +568 @@
     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;
+      }
+    }
+  }
+}
+
+//---------------------------------------------------------------------------