TTruthVerticesModule.cxx

Go to the documentation of this file.

#include "TTruthVerticesModule.hxx"

#include <sstream>
// Root
#include <TClonesArray.h>
// oaEvent
#include <TG4PrimaryParticle.hxx>
#include <TG4PrimaryVertex.hxx>
#include <TGeometryId.hxx>
// oaGeomInfo
#include <TGeomInfo.hxx>
#include <TP0DGeom.hxx>

ClassImp(ND::TTruthVerticesModule::TTruthVertex);

#define CVSTAG \
  "\
  $Name: 7.0$"
#define CVSID \
  "\
  $Id: eventAnalysis TTruthVerticesModule.cxx,2024/03/20:09:46:11,Alexander_J_Finch,lapw.lancs.ac.uk $"

ND::TTruthVerticesModule::TTruthVerticesModule(const char* name,
                                               const char* title)
    : fMaxNVertices(500),
      fNVtx(0),
      fNVtxFGD1(0),
      fNVtxFGD2(0),
      fNVtxP0D(0),
      fNVtxDsECal(0),
      fNVtxBrlECal(0),
      fNVtxP0DECal(0),
      fNVtxTPC1(0),
      fNVtxTPC2(0),
      fNVtxTPC3(0),
      fNVtxSMRD(0),
      fNVtxRestOfOffAxis(0),
      fNVtxINGRID(0),
      fNVtxOther(0),
      fVertices(new TClonesArray("ND::TTruthVerticesModule::TTruthVertex",
                                 fMaxNVertices)) {
  SetNameTitle(name, title);
  fIsEnabled = kTRUE;
  fDescription = "Module storing the truth information for primary vertices";
  fCVSTagName = CVSTAG;
  fCVSID = CVSID;
}

ND::TTruthVerticesModule::~TTruthVerticesModule() { delete fVertices; }

void ND::TTruthVerticesModule::InitializeBranches() {
  fOutputTree->Branch("NVtx", &fNVtx, "NVtx/I", fBufferSize)->SetTitle("The total number of vertices recorded");
  fOutputTree->Branch("NVtxFGD1", &fNVtxFGD1, "NVtxFGD1/I", fBufferSize)->SetTitle("The number of vertices recorded in FGD 1");
  fOutputTree->Branch("NVtxFGD2", &fNVtxFGD2, "NVtxFGD2/I", fBufferSize)->SetTitle("The number of vertices recorded in FGD 2");
  fOutputTree->Branch("NVtxSFG", &fNVtxSFG, "NVtxSFG/I", fBufferSize)->SetTitle("The number of vertices recorded in SFG");
  fOutputTree->Branch("NVtxP0D", &fNVtxP0D, "NVtxP0D/I", fBufferSize)->SetTitle("The number of vertices recorded in the P0D");
  fOutputTree->Branch("NVtxDsECal", &fNVtxDsECal, "NVtxDsECal/I", fBufferSize)->SetTitle("The number of vertices recorded in the Downstream ECal");
  fOutputTree->Branch("NVtxBrlECal", &fNVtxBrlECal, "NVtxBrlECal/I",
                      fBufferSize)->SetTitle("The number of vertices recorded in the Barrel ECal");
  fOutputTree->Branch("NVtxP0DECal", &fNVtxP0DECal, "NVtxP0DECal/I",
                      fBufferSize)->SetTitle("The number of vertices recorded in the P0D ECal");
  fOutputTree->Branch("NVtxTPC1", &fNVtxTPC1, "NVtxTPC1/I", fBufferSize)->SetTitle("The number of vertices recorded in TPC 1");
  fOutputTree->Branch("NVtxTPC2", &fNVtxTPC2, "NVtxTPC2/I", fBufferSize)->SetTitle("The number of vertices recorded in TPC 2");
  fOutputTree->Branch("NVtxTPC3", &fNVtxTPC3, "NVtxTPC3/I", fBufferSize)->SetTitle(" The number of vertices recorded in TPC 3");
  fOutputTree->Branch("NVtxSMRD", &fNVtxSMRD, "NVtxSMRD/I", fBufferSize)->SetTitle("The number of vertices recorded in the SMRD");
  fOutputTree->Branch("NVtxRestOfOffAxis", &fNVtxRestOfOffAxis,
                      "NVtxRestOfOffAxis/I", fBufferSize)->SetTitle(" The number of vertices recorded in the rest of the off-axis detector");
  fOutputTree->Branch("NVtxINGRID", &fNVtxINGRID, "NVtxINGRID/I", fBufferSize)->SetTitle("The number of vertices recorded in INGRID");
  fOutputTree->Branch("NVtxOther", &fNVtxOther, "NVtxOther/I", fBufferSize)->SetTitle("The number of vertices recorded anywhere which does not fall  into the other available categories");

  fOutputTree->Branch("Vertices", "TClonesArray", &fVertices, fBufferSize,
                      fSplitLevel)->SetTitle("The TClonesArray storing the TTruthVertex objects holding the information relating to each vertex.");
}

Bool_t ND::TTruthVerticesModule::ProcessFirstEvent(ND::TND280Event& event) {
  bool IsDataFile = event.GetContext().IsDetector();
  if (IsDataFile) {
    throw ND::EDataFile();
  }
  return true;
}

bool ND::TTruthVerticesModule::FillTree(ND::TND280Event& event) {
  if (!event.GetContext().IsMC()) {
    ND280Error(
        "Event Context reports event is non-MC. Assuming entire file is "
        "non-MC.");
    return false;  // disable module
  }

  fNVtx = 0;
  fNVtxFGD1 = 0;
  fNVtxFGD2 = 0;
  fNVtxSFG  = 0;
  fNVtxP0D = 0;
  fNVtxDsECal = 0;
  fNVtxBrlECal = 0;
  fNVtxP0DECal = 0;
  fNVtxTPC1 = 0;
  fNVtxTPC2 = 0;
  fNVtxTPC3 = 0;
  fNVtxSMRD = 0;
  fNVtxRestOfOffAxis = 0;
  fNVtxINGRID = 0;
  fNVtxOther = 0;
  fVertices->Clear();

  if (!event.Has<ND::TG4PrimaryVertexContainer>("truth/G4PrimVertex00")) {
    ND280Warn("No vertices container found, skipping event.");
    return true;  // this happens occasionally, so don't disable the module
  }

  const ND::THandle<ND::TG4PrimaryVertexContainer>& vertices =
      event.Get<ND::TG4PrimaryVertexContainer>("truth/G4PrimVertex00");

  std::vector<ND::TG4PrimaryVertex>::const_iterator vertices_it =
      vertices->begin();
  std::vector<ND::TG4PrimaryVertex>::const_iterator vertices_end;

  if (vertices->size() <= fMaxNVertices) {
    vertices_end = vertices->end();
  } else {
    ND280Warn("Event contains " << vertices->size() << " vertices."
                                << " Only the first " << fMaxNVertices
                                << " will be saved.");
    vertices_end = vertices->begin();
    vertices_end += fMaxNVertices;
  }

  for (; vertices_it != vertices_end; ++vertices_it) {
    // Firstly, check if the vertex is the geantino vertex, if it is, don't
    // bother with it!
    if (ND::TTruthVerticesModule::IsGeantinoVtx((*vertices_it))) continue;

    TTruthVertex* vtxToFill = new ((*fVertices)[fNVtx]) TTruthVertex();
    ++fNVtx;

    vtxToFill->ID = vertices_it->GetInteractionNumber();
    vtxToFill->Position = vertices_it->GetPosition();
    vtxToFill->Generator = vertices_it->GetGeneratorName();
    vtxToFill->ReactionCode = vertices_it->GetReaction();

    TLorentzVector pos = vertices_it->GetPosition();

    gGeoManager->FindNode(pos.X(), pos.Y(), pos.Z());
    std::string path = gGeoManager->GetPath();
    vtxToFill->Subdetector = eventAnalysisEnums::PathToSubdetector(path);

    // Fill vertex's fiducial distance and increment the relevant
    // vertex counter
    switch (vtxToFill->Subdetector) {
      case eventAnalysisEnums::kFGD1:
        ++fNVtxFGD1;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kFGD2:
        ++fNVtxFGD2;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kSFG:
        ++fNVtxSFG;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kP0D:
        ++fNVtxP0D;
        vtxToFill->Fiducial =
            TGeomInfo::P0D().WaterFiducialDistance(pos.X(), pos.Y(), pos.Z()) *
            unit::mm;
        break;
      case eventAnalysisEnums::kDsECal:
        ++fNVtxDsECal;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kBrlECalBottom:
      case eventAnalysisEnums::kBrlECalLeft:
      case eventAnalysisEnums::kBrlECalTop:
      case eventAnalysisEnums::kBrlECalRight:
        ++fNVtxBrlECal;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kP0DECalBottom:
      case eventAnalysisEnums::kP0DECalLeft:
      case eventAnalysisEnums::kP0DECalTop:
      case eventAnalysisEnums::kP0DECalRight:
        ++fNVtxP0DECal;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kTPC1:
        ++fNVtxTPC1;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kTPC2:
        ++fNVtxTPC2;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kTPC3:
        ++fNVtxTPC3;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kSMRD:
        ++fNVtxSMRD;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kOffAxis:
        ++fNVtxRestOfOffAxis;
        vtxToFill->Fiducial = -1.0;
        break;
      case eventAnalysisEnums::kINGRID:
        ++fNVtxINGRID;
        vtxToFill->Fiducial = -1.0;
        break;
      default:
        ++fNVtxOther;
        vtxToFill->Fiducial = -1.0;
        break;
    }
    // END SWITCH ON SUBDETECTOR

    ///
    /// Fill primary particle trajectory information.
    /// These are particles leaving the vertex.

    std::vector<ND::TG4PrimaryParticle>::const_iterator primaries_it =
        vertices_it->GetPrimaryParticles().begin();
    std::vector<ND::TG4PrimaryParticle>::const_iterator primaries_end =
        vertices_it->GetPrimaryParticles().end();
    vtxToFill->NPrimaryTraj = 0;

    for (; primaries_it != primaries_end; ++primaries_it) {
      vtxToFill->PrimaryTrajIDs.push_back(primaries_it->GetTrackId());
      ++(vtxToFill->NPrimaryTraj);
    }

    // For the info particles, it is a bit hard to specify the
    // behaviour without knowing exactly what the generator is
    // giving us.
    // For now, just look at particles with Id == -1 (-2 means
    // stuff going on inside the nucleus). Look for neutrinos and
    // non-neutrinos, assuming the latter are targets.
    // Loop through all info vertices and their particles until both
    // neutrino and target have been seen, then stop.

    bool found_neutrino = false;
    bool found_target = false;

    std::vector<ND::TG4PrimaryVertex>::const_iterator info_it =
        vertices_it->GetInfoVertex().begin();
    std::vector<ND::TG4PrimaryVertex>::const_iterator info_end =
        vertices_it->GetInfoVertex().end();

    for (; info_it != info_end; ++info_it) {
      std::vector<ND::TG4PrimaryParticle>::const_iterator particles_it =
          info_it->GetPrimaryParticles().begin();
      std::vector<ND::TG4PrimaryParticle>::const_iterator particles_end =
          info_it->GetPrimaryParticles().end();

      Int_t temp_pdg;

      for (; particles_it != particles_end; ++particles_it) {
        if (particles_it->GetTrackId() ==
            -1)  // is an incoming particle (including target)
        {
          temp_pdg = particles_it->GetPDGCode();
          switch (temp_pdg) {
            case 12:   // fall through
            case 14:   // fall through
            case 16:   // fall through
            case -12:  // fall through
            case -14:  // fall through
            case -16:
              vtxToFill->NeutrinoPDG = temp_pdg;
              vtxToFill->NeutrinoMomentum = particles_it->GetMomentum();
              found_neutrino = true;
              break;
            default:
              vtxToFill->TargetPDG = temp_pdg;
              vtxToFill->TargetMomentum = particles_it->GetMomentum();
              found_target = true;
              break;
          }
          if (found_neutrino && found_target) break;
        }
      }
      if (found_neutrino && found_target) break;
    }
    // END LOOP OVER INFO VERTICES
  }
  // END LOOP OVER VERTICES

  fVertices->Sort();

  return true;
}

//
// Truth Vertex
// =====================================================================
//

ND::TTruthVerticesModule::TTruthVertex::TTruthVertex()
    : ID(-1),
      Position(-999999.9, -999999.9, -999999.9, -999999.9),
      Generator("Not Set"),
      ReactionCode("Not Set"),
      Subdetector(eventAnalysisEnums::kNSubdetectors),
      Fiducial(-999999.9),
      NPrimaryTraj(-1),
      NeutrinoPDG(-1),
      NeutrinoMomentum(-999999.9, -999999.9, -999999.9, -999999.9),
      TargetPDG(-1),
      TargetMomentum(-999999.9, -999999.9, -999999.9, -999999.9) {}

Int_t ND::TTruthVerticesModule::TTruthVertex::Compare(
    const TObject* obj) const {
  // Sorry about the casting, but I don't want to slow down sorting
  // by checking that everything compared with is of the correct
  // type. In any case, you get what you ask for if you try and
  // use this with anything other than another vertex.

  if (this->ID <
      static_cast<const ND::TTruthVerticesModule::TTruthVertex*>(obj)->ID) {
    return -1;
  } else if (this->ID >
             static_cast<const ND::TTruthVerticesModule::TTruthVertex*>(obj)
                 ->ID) {
    return 1;
  } else  // they must be equal
  {
    return 0;
  }
}

Bool_t ND::TTruthVerticesModule::IsGeantinoVtx(ND::TG4PrimaryVertex vtx) {
  if (vtx.GetPrimaryParticles().size() == 1 &&
      vtx.GetInteractionNumber() == 0 && vtx.GetProbability() == 0 &&
      vtx.GetCrossSection() == 0 && vtx.GetDiffCrossSection() == 0 &&
      vtx.GetPrimaryParticles().begin()->GetPDGCode() == 0)
    return true;

  return false;
}