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

File:

• classes/DelphesHepMC3Reader.cc (modified) (14 diffs)

classes/DelphesHepMC3Reader.cc

@@ -56 +56 @@
 DelphesHepMC3Reader::DelphesHepMC3Reader() :
   fInputFile(0), fBuffer(0), fPDG(0),
-  fVertexCounter(-2), fParticleCounter(-1)
+  fVertexCounter(-1), fParticleCounter(-1)
 {
   fBuffer = new char[kBufferSize];
@@ -81 +81 @@
 void DelphesHepMC3Reader::Clear()
 {
-  fWeights.clear();
+  fWeight.clear();
   fMomentumCoefficient = 1.0;
   fPositionCoefficient = 1.0;
-  fVertexCounter = -2;
+  fVertexCounter = -1;
   fParticleCounter = -1;
-  fVertices.clear();
-  fParticles.clear();
-  fInVertexMap.clear();
-  fOutVertexMap.clear();
-  fMotherMap.clear();
+  fVertexMap.clear();
   fDaughterMap.clear();
 }
@@ -98 +94 @@
 bool DelphesHepMC3Reader::EventReady()
 {
-  return (fVertexCounter == -1) && (fParticleCounter == 0);
+  return (fParticleCounter == 0);
 }
 
@@ -123 +119 @@
     Clear();
 
+    fX = 0.0;
+    fY = 0.0;
+    fZ = 0.0;
+    fT = 0.0;
+
     rc = bufferStream.ReadInt(fEventNumber)
       && bufferStream.ReadInt(fVertexCounter)
@@ -167 +168 @@
     while(bufferStream.ReadDbl(weight))
     {
-      fWeights.push_back(weight);
+      fWeight.push_back(weight);
     }
   }
@@ -256 +257 @@
   else if(key == 'V')
   {
-    fParticles.clear();
-
     fX = 0.0;
     fY = 0.0;
@@ -263 +262 @@
     fT = 0.0;
 
+    fM1 = 0;
+    fM2 = 0;
+
     rc = bufferStream.ReadInt(fVertexCode)
       && bufferStream.ReadInt(fVertexStatus);
@@ -284 +286 @@
     while(bufferStream.ReadInt(code))
     {
-      fParticles.push_back(code);
-      bufferStream.FindChr(',');
-    }
+      if(code < fM1 || fM1 == 0) fM1 = code;
+      if(code > fM2) fM2 = code;
+      fVertexMap[code] = fVertexCode;
+    }
+
+    if(fM1 == fM2) fM2 = 0;
 
     if(bufferStream.FindChr('@'))
@@ -302 +307 @@
       }
     }
-
-    AnalyzeVertex(factory, fVertexCode);
   }
   else if(key == 'P' && fParticleCounter > 0)
@@ -326 +329 @@
     }
 
-    AnalyzeParticle(factory);
+    itDaughterMap = fDaughterMap.find(fOutVertexCode);
+    if(itDaughterMap == fDaughterMap.end())
+    {
+      fDaughterMap[fOutVertexCode] = make_pair(fParticleCode, fParticleCode);
+    }
+    else
+    {
+      itDaughterMap->second.second = fParticleCode;
+    }
+
+    AnalyzeParticle(factory, allParticleOutputArray,
+      stableParticleOutputArray, partonOutputArray);
   }
 
   if(EventReady())
   {
-    FinalizeParticles(allParticleOutputArray, stableParticleOutputArray, partonOutputArray);
+    FinalizeParticles(allParticleOutputArray);
   }
 
@@ -349 +363 @@
   element->ProcessID = fProcessID;
   element->MPI = fMPI;
-  element->Weight = fWeights.size() > 0 ? fWeights[0] : 1.0;
+  element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
   element->CrossSection = fCrossSection;
   element->CrossSectionError = fCrossSectionError;
@@ -373 +387 @@
 {
   Weight *element;
-  vector<double>::const_iterator itWeights;
-
-  for(itWeights = fWeights.begin(); itWeights != fWeights.end(); ++itWeights)
+  vector<double>::const_iterator itWeight;
+
+  for(itWeight = fWeight.begin(); itWeight != fWeight.end(); ++itWeight)
   {
     element = static_cast<Weight *>(branch->NewEntry());
 
-    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,
+    element->Weight = *itWeight;
+  }
+}
+
+//---------------------------------------------------------------------------
+
+void DelphesHepMC3Reader::AnalyzeParticle(DelphesFactory *factory,
+  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, 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);
-      }
-    }
-  }
+  int i, index;
 
   for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
@@ -549 +472 @@
     candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
 
-    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)
+    index = i + 1;
+
+    itVertexMap = fVertexMap.find(index);
+    if(itVertexMap != fVertexMap.end())
+    {
+      index = itVertexMap->second;
+    }
+
+    itDaughterMap = fDaughterMap.find(index);
+    if(itDaughterMap == fDaughterMap.end())
     {
       candidate->D1 = -1;
@@ -568 +488 @@
     else
     {
-      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;
-      }
-    }
-  }
-}
-
-//---------------------------------------------------------------------------
+      candidate->D1 = itDaughterMap->second.first - 1;
+      candidate->D2 = itDaughterMap->second.second - 1;
+    }
+  }
+}
+
+//---------------------------------------------------------------------------