Bugfixes in SLDCorrection

Analysis/PIDTools/include/ComprehensivePIDProcessor.h

@@ -151,6 +151,7 @@ class ComprehensivePIDProcessor : public Processor {
   bool _modeExtract = false;
   bool _modeTrain = false;
   bool _modeInfer = false;
+  bool _addMCPID = false;
 
   std::string _TTreeFileName{};
   std::vector<std::string> _inputAlgoSpecs{};

Analysis/PIDTools/src/ComprehensivePIDProcessor.cc

@@ -51,6 +51,11 @@ ComprehensivePIDProcessor::ComprehensivePIDProcessor() : Processor("Comprehensiv
                              "you need to provide the reference and weight files; default: false.",
                              _modeInfer, false);
 
+  registerProcessorParameter("addMCPID",
+                             "Set true to add a PID output called 'MCPID' to the PFOs which is directly derived from "
+                             "MC information (no inference), not compatible with training mode; default: false.",
+                             _addMCPID, false);
+
   registerProcessorParameter(
       "TTreeFileName",
       "Name of the root file in which the TTree with all observables is stored; in case of extraction it is an "
@@ -146,6 +151,10 @@ void ComprehensivePIDProcessor::init() {
     sloE << "I cannot train and infer in the same process!" << std::endl;
     throw std::runtime_error("mode error");
   }
+  if (_modeTrain && _addMCPID) {
+    sloE << "I cannot train and add MCPID in the same process!" << std::endl;
+    throw std::runtime_error("mode error");
+  }
 
   if (_TTreeFileName == "") {
     sloM << "TTreeFileName is empty, no root file is created" << std::endl;
@@ -410,7 +419,7 @@ void ComprehensivePIDProcessor::processEvent(LCEvent* evt) {
   if (_nAlgos == 0)
     return;
 
-  if (_modeExtract || _modeInfer) {
+  if (_modeExtract || _modeInfer || _addMCPID) {
     LCCollection *col_pfo{}, *col_pfo2mc{};
 
     try {
@@ -427,7 +436,7 @@ void ComprehensivePIDProcessor::processEvent(LCEvent* evt) {
     PIDHandler pidh(col_pfo);
     std::vector<int> algoID{};
     std::vector<int> allPDGs{};
-    if (_modeInfer) {
+    if (_modeInfer || _addMCPID) {
       std::vector<std::string> PDGness{};
       for (int pdg : _signalPDGs) {
         std::stringstream pn;
@@ -441,8 +450,12 @@ void ComprehensivePIDProcessor::processEvent(LCEvent* evt) {
         PDGness.push_back(pn.str());
         allPDGs.push_back(pdg);
       }
-      for (int m = 0; m < _nModels; ++m)
-        algoID.push_back(pidh.addAlgorithm(_trainModelNames[m], PDGness));
+      if (_modeInfer)
+        for (int m = 0; m < _nModels; ++m)
+          algoID.push_back(pidh.addAlgorithm(_trainModelNames[m], PDGness));
+
+      if (_addMCPID)
+        algoID.push_back(pidh.addAlgorithm("MCPID", PDGness));
     }
 
     for (int i = 0; i < n_pfo; ++i) {
@@ -567,6 +580,20 @@ void ComprehensivePIDProcessor::processEvent(LCEvent* evt) {
             break;
           }
       }
+
+      if (_addMCPID) {
+        sloD << "start adding MCPID" << std::endl;
+        std::vector<float> mcPID{};
+
+        for (int ipdg = 0; ipdg < (int)allPDGs.size(); ++ipdg) {
+          if (allPDGs[ipdg] == abs(MCPDG))
+            mcPID.push_back(1);
+          else
+            mcPID.push_back(0);
+        }
+
+        pidh.setParticleID(pfo, 0, abs(MCPDG), 0, pidh.getAlgorithmID("MCPID"), mcPID);
+      }
     } // for pfos
   }
 
@@ -633,8 +660,9 @@ void ComprehensivePIDProcessor::end() {
 
   for (int i = 0; i < _nAlgos; ++i)
     delete _inputAlgorithms[i];
-  for (int i = 0; i < _nModels; ++i)
-    delete _trainingModels[i];
+  if (_modeTrain || _modeInfer)
+    for (int i = 0; i < _nModels; ++i)
+      delete _trainingModels[i];
 
   sloM << "Numer of Events: " << _nEvt << "    Numer of PFOs: " << _nPFO << std::endl;
   sloM << "-------------------------------------------------" << std::endl;

Analysis/SLDCorrection/README.md

@@ -71,4 +71,18 @@ Green histograms are used to estimate the overall uncertainty due to neutrino co
 
 sigmaAlphaNu = 0.100 rad , sigmaENu = 4.0 GeV
 
-In the end, neutrinos are assumed to be massless and sigmaAlphaNu and sigmaENu are transformed to a covariance matrix in (p,E) space.
+In the end, neutrinos are assumed to be massless and sigmaAlphaNu and sigmaENu are transformed to a covariance matrix in (p,E) space.
+
+### SLDStatus Reference
+| SLDStatus | Interpretation                                  |
+|-----------|-------------------------------------------------|
+|         0 | PrimaryVertex collection found, but empty       |
+|         1 | Reconstructed Lepton is not found               |
+|         2 | Reconstructed Lepton doesn't belong to any jet  |
+|         3 | Infer SLD from reconstructed lepton in primary vertex    |
+|         4 | Infer SLD from reconstructed lepton in a secondary (BuildUp) Vertex.    |
+|         5 | Infer SLD frmo BuildUp Vertex in a jet. Intersection point of Lepton and other BuildUp Vertices in jet is used as vertex of semi-leptonic decay.   |
+|         6 | Inferred from lepton track alone.               |
+|         7 | other                                           |
+
+In the h_SLDStatus histogram, these values are inserted with a shift of -0.5

Analysis/SLDCorrection/include/SLDCorrection.h

@@ -14,20 +14,19 @@
 #include <EVENT/Track.h>
 #include <IMPL/ReconstructedParticleImpl.h>
 #include <IMPL/VertexImpl.h>
+#include <TF1.h>
+#include <TFile.h>
+#include <TH1F.h>
+#include <TH1I.h>
+#include <TH2F.h>
+#include <TH2I.h>
+#include <TTree.h>
 #include <marlin/Global.h>
 #include <marlin/Processor.h>
 #include <math.h>
 #include <set>
 #include <string>
 #include <vector>
-class TFile;
-class TDirectory;
-class TH1F;
-class TH1I;
-class TH2I;
-class TH2F;
-class TTree;
-class TF1;
 
 using namespace lcio;
 using namespace marlin;
@@ -59,24 +58,24 @@ class SLDCorrection : public Processor {
   bool hasUpStreamSLDecay(const MCP& parentHadron);
 
   //	checkBHadronSLDecay checks the flavour of the Hadron decays semi-leptonicall, if B-Hadron: true, if not: false
-  bool checkBHadronSLDecay(const MCP& SLDLepton);
+  bool checkBHadronSLDecay(const MCP& parentHadron);
 
   //	checkCHadronSLDecay checks the flavour of the Hadron decays semi-leptonicall, if C-Hadron: true, if not: false
-  bool checkCHadronSLDecay(const MCP& SLDLepton);
+  bool checkCHadronSLDecay(const MCP& parentHadron);
 
   //	checkTauLeptonSLDecay checks the flavour of the Hadron decays semi-leptonicall, if tau-lepton: true, if not:
   // false
-  bool checkTauLeptonSLDecay(const MCP& SLDLepton);
+  bool checkTauLeptonSLDecay(const MCP& parentHadron);
 
   //	doSLDCorrection prepares all inputs for nu-correction and makes output LCIO elements (SLDVertex, associated RP
   // of SLDVertex, solutions for neutrino as Reconstructed Particle, etc)
-  virtual void doSLDCorrection(EVENT::LCEvent* pLCEvent, const MCP& SLDLepton, VertexVector& semiLeptonicVertices,
-                               PFOVector& semiLeptonicVertexRecoParticles, PFOVector& jetsOfSemiLeptonicDecays,
-                               PFOVectorVector& neutrinos, IntVector& sldStatus, IntVector& pvaStatus,
-                               IntVector& solutionSigns, MCPVector& trueNeutrinos);
+  virtual void doSLDCorrection(EVENT::LCEvent* pLCEvent, const MCP& SLDLepton, size_t parentHadronIDx,
+                               VertexVector& semiLeptonicVertices, PFOVector& semiLeptonicVertexRecoParticles,
+                               PFOVector& jetsOfSemiLeptonicDecays, PFOVectorVector& neutrinos, IntVector& sldStatus,
+                               IntVector& pvaStatus, IntVector& solutionSigns, MCPVector& trueNeutrinos);
 
   //	showTrueParameters prints true input for nu-correction
-  void showTrueParameters(const MCP& SLDLepton);
+  void showTrueParameters(const MCP& SLDLepton, size_t parentHadronIDx);
 
   //	getNeutrinoFourMomentum corrects neutrino energy using rapidity-based approach
   TLorentzVector getNeutrinoFourMomentum(const TVector3& flightDirection, const TLorentzVector& visibleFourMomentum,
@@ -94,7 +93,7 @@ class SLDCorrection : public Processor {
                                                        const double& parentHadronMass, const float& solutionSign);
 
   //	getTrueNeutrinogets true four-momentum of neutrino
-  MCP getTrueNeutrino(const MCP& SLDLepton);
+  MCP getTrueNeutrino(const MCP& SLDLepton, size_t parent_idx);
 
   //	fillTrueRecoFourMomentummakes performance evaluation variables stored in root tree
   void fillTrueRecoFourMomentum(
@@ -582,6 +581,6 @@ class SLDCorrection : public Processor {
   TH1F* h_FlightDirectionError{};
   TH1F* h_distRecoLeptonToDownStreamVertex{};
   TFile* m_pTFile{};
-  TTree* m_pTTree1{};
+  TTree* m_pTTree1 = new TTree("SLDCorrection", "SLDCorrection");
 };
 #endif

Analysis/SLDCorrection/src/FlightDirection.cc

@@ -48,6 +48,12 @@ int getRecoFlightDirection(const RecoParticle& linkedRecoLepton, TVector3& recoF
   daughterHadronFlightDistance = 0.0;
   daughterHadronFlightDirection = TVector3(0.0, 0.0, 0.0);
   sldVertexPosition.clear();
+  if (linkedRecoLepton->getTracks().size() == 0) {
+    streamlog_out(DEBUG1) << "	(" << SLDStatus << ") No track for linkedRecoLepton. SLDCorrection aborts."
+                          << std::endl;
+    return SLDStatus;
+  }
+
   if (recoLeptonIsInVertex) {
     SLDStatus = 4;
     streamlog_out(DEBUG1) << "	(" << SLDStatus
@@ -151,8 +157,6 @@ int getRecoFlightDirection(const RecoParticle& linkedRecoLepton, TVector3& recoF
       //	thirdVertex = testVertex;
       // }
     }
-    if (thirdVertex != NULL) {
-    }
   } else {
     streamlog_out(DEBUG1) << "	There is NO BuildUp Vertex in jet" << std::endl;
     streamlog_out(DEBUG1) << "	Other scenarios are investigated, but Nu-correction will not be performed" << std::endl;