TGRooTrackerVtxModule.cxx

Go to the documentation of this file.

#include <cassert>
#include <cmath>

#include "TG4PrimaryParticle.hxx"
#include "TG4PrimaryVertex.hxx"
#include "TGRooTrackerVtxModule.hxx"
#include "TGeomInfo.hxx"
#include "TJNuBeamFlux.hxx"
#include "TND280Log.hxx"
#include "TOADatabase.hxx"
#include "TP0DGeom.hxx"

ND::TGRooTrackerVtxModule::TGRooTrackerVtxModule(const char *name,
                                                 const char *title) {
  fVtx = new TClonesArray("ND::GRooTrackerVtx", 100);
  SetNameTitle(name, title);
  ///< Enable this module by default:
  fIsEnabled = kTRUE;
  fDescription = "Genie specific module which fills mc-truth vertex info";
  fCurrGenieVtx = new GRooTrackerVtx();
  // ND::TND280Log::Configure();
}

ND::TGRooTrackerVtxModule::~TGRooTrackerVtxModule() {}

Bool_t ND::TGRooTrackerVtxModule::IsKnownGenerator(const char *name) {
  std::string name2compare(name);
  bool is_genie1 = (name2compare.find("genie") != name2compare.npos);
  bool is_genie2 = (name2compare.find("gRooTracker") != name2compare.npos);
  bool is_genie = (is_genie1 || is_genie2);
  ND280Debug("Returning is generator = " << is_genie
                                         << " for generator name = " << name);
  return is_genie;
}

Bool_t ND::TGRooTrackerVtxModule::CheckMismatch(ND::TG4PrimaryVertex vtx) {
  if (fCurrGenieVtx == NULL) {
    ND280Debug("Could not compare pass-through to G4 primary vertex as"
               << " there is no passed through vertex to compare to!");
    return false;
  }

  ///< Compare the event codes.
  bool evtcode =
      (fCurrGenieVtx->EvtCode->GetString().CompareTo(vtx.GetReaction()) == 0);
  ND280Debug("Comparing " << fCurrGenieVtx->EvtCode->GetString().Data()
                          << " and " << vtx.GetReaction());
  ///< And the vertex position.
  double smallnum = 0.000001;
  bool position =
      (std::abs(fCurrGenieVtx->EvtVtx[0] - vtx.GetPosition().X()) < smallnum) &&
      (std::abs(fCurrGenieVtx->EvtVtx[1] - vtx.GetPosition().Y()) < smallnum) &&
      (std::abs(fCurrGenieVtx->EvtVtx[2] - vtx.GetPosition().Z()) < smallnum);
  ND280Debug("Pass through vertex position check bool passed = " << position);
  ///< In future can add more but for now this should suffice along
  ///< with the validation scripts.
  if (position && evtcode) {
    ND280Debug("Passed through vertex passed mismatch test!");
    return true;
  }
  if (ND::TGRooTrackerVtxModule::IsGeantinoVtx(vtx)) return true;
  ND280Error("SEVERE ERROR! Mismatch between G4Vertex and pass-through info!");
  return false;
}

Bool_t ND::TGRooTrackerVtxModule::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;
}

void ND::TGRooTrackerVtxModule::FillVtx(int entry) {
  ND280Verbose(
      "Filling fCurrGenieVtx with OrigEvtNum = " << fOrigTreeEntryNumber);

  fCurrGenieVtx->Reset();
  fRooTrackerTree->GetEntry(entry);

  ///< first fill up all the book-keeping info
  fCurrGenieVtx->OrigFileName->SetString(fOrigInputFileName.c_str());
  fCurrGenieVtx->OrigTreeName->SetString(fInputTreeName.c_str());
  fCurrGenieVtx->OrigEvtNum = fOrigTreeEntryNumber;
  fCurrGenieVtx->OrigTreePOT = fOrigInputTreePOT;
  fCurrGenieVtx->OrigTreeEntries = fOrigInputTreeEntries;
  fCurrGenieVtx->TimeInSpill = fTimeInSpill;
  fCurrGenieVtx->TruthVertexID = fTruthVertexID;

  ///< Need to convert to the nd280 length units
  fCurrGenieVtx->EvtVtx[0] *= 1000.0;
  fCurrGenieVtx->EvtVtx[1] *= 1000.0;
  fCurrGenieVtx->EvtVtx[2] *= 1000.0;

  ///< Fill in the geometry path for the current vertex.
  gGeoManager->FindNode(fCurrGenieVtx->EvtVtx[0], fCurrGenieVtx->EvtVtx[1],
                        fCurrGenieVtx->EvtVtx[2]);

  ///< Get the name of the generator that generated this event.
  fCurrGenieVtx->GeneratorName->SetString(fGeneratorName.c_str());

  ///< Get the geometry path for this event.
  fCurrGenieVtx->GeomPath->SetString(gGeoManager->GetPath());

  ///< Now create another fRooTracker vtx to add to the vertex
  ///< container.
  GRooTrackerVtx *fGRooTrackerVtx;
  fGRooTrackerVtx = new ((*fVtx)[fNVtx++]) GRooTrackerVtx();

  ///< Copy the current genie event. Make sure this is done last so we copy a
  ///< complete event.
  fGRooTrackerVtx->Copy(fCurrGenieVtx);
}

void ND::TGRooTrackerVtxModule::SetGeneratorTreeAddresses() {
  // Set the flux branch addresses common to both generators
  ND::JNuBeamFlux *flux = (ND::JNuBeamFlux *)fCurrGenieVtx;
  ND::TRooTrackerVtxModuleBase::SetFluxTreeAddresses(flux);

  // Now set the generator specific branch addresses
  fRooTrackerTree->SetBranchAddress("G2NeutEvtCode",
                                    &(fCurrGenieVtx->G2NeutEvtCode));
  fRooTrackerTree->SetBranchAddress("EvtCode", &(fCurrGenieVtx->EvtCode));
  fRooTrackerTree->SetBranchAddress("EvtNum", &(fCurrGenieVtx->EvtNum));
  fRooTrackerTree->SetBranchAddress("EvtXSec", &(fCurrGenieVtx->EvtXSec));
  fRooTrackerTree->SetBranchAddress("EvtDXSec", &(fCurrGenieVtx->EvtDXSec));
  fRooTrackerTree->SetBranchAddress("EvtWght", &(fCurrGenieVtx->EvtWght));
  fRooTrackerTree->SetBranchAddress("EvtProb", &(fCurrGenieVtx->EvtProb));
  fRooTrackerTree->SetBranchAddress("EvtVtx", (fCurrGenieVtx->EvtVtx));
  fRooTrackerTree->SetBranchAddress("StdHepN", &(fCurrGenieVtx->StdHepN));
  fRooTrackerTree->SetBranchAddress("StdHepPdg",
                                    &(fCurrGenieVtx->StdHepPdgTemp));
  fRooTrackerTree->SetBranchAddress("StdHepRescat",
                                    &(fCurrGenieVtx->StdHepRescatTemp));
  fRooTrackerTree->SetBranchAddress("StdHepStatus",
                                    &(fCurrGenieVtx->StdHepStatusTemp));
  fRooTrackerTree->SetBranchAddress("StdHepX4", (fCurrGenieVtx->StdHepX4));
  fRooTrackerTree->SetBranchAddress("StdHepP4", (fCurrGenieVtx->StdHepP4));
  fRooTrackerTree->SetBranchAddress("StdHepPolz", (fCurrGenieVtx->StdHepPolz));
  fRooTrackerTree->SetBranchAddress("StdHepFd", (fCurrGenieVtx->StdHepFdTemp));
  fRooTrackerTree->SetBranchAddress("StdHepLd", (fCurrGenieVtx->StdHepLdTemp));
  fRooTrackerTree->SetBranchAddress("StdHepFm", (fCurrGenieVtx->StdHepFmTemp));
  fRooTrackerTree->SetBranchAddress("StdHepLm", (fCurrGenieVtx->StdHepLmTemp));
}


void ND::TGRooTrackerVtxModule::FillConfigTree(TTree* configTree){

  if(fFilledConfigTree){
    ND280NamedWarn("TGRooTrackerVtxModule", "Module: " << this->GetName()
      << " has already written the config tree this run.");
  }

  int kGStdHepNPmax_hldr = kGStdHepNPmax;
  int kGStdHepIdxPx_hldr = kGStdHepIdxPx;
  int kGStdHepIdxPy_hldr = kGStdHepIdxPy;
  int kGStdHepIdxPz_hldr = kGStdHepIdxPz;
  int kGStdHepIdxE_hldr = kGStdHepIdxE;
  int kGStdHepIdxX_hldr = kGStdHepIdxX;
  int kGStdHepIdxY_hldr = kGStdHepIdxY;
  int kGStdHepIdxZ_hldr = kGStdHepIdxZ;
  int kGStdHepIdxT_hldr = kGStdHepIdxT;

  configTree->Branch("kGStdHepNPmax", &kGStdHepNPmax_hldr, "kGStdHepNPmax/I");
  configTree->Branch("kGStdHepIdxPx", &kGStdHepIdxPx_hldr, "kGStdHepIdxPx/I");
  configTree->Branch("kGStdHepIdxPy", &kGStdHepIdxPy_hldr, "kGStdHepIdxPy/I");
  configTree->Branch("kGStdHepIdxPz", &kGStdHepIdxPz_hldr, "kGStdHepIdxPz/I");
  configTree->Branch("kGStdHepIdxE", &kGStdHepIdxE_hldr, "kGStdHepIdxE/I");
  configTree->Branch("kGStdHepIdxX", &kGStdHepIdxX_hldr, "kGStdHepIdxX/I");
  configTree->Branch("kGStdHepIdxY", &kGStdHepIdxY_hldr, "kGStdHepIdxY/I");
  configTree->Branch("kGStdHepIdxZ", &kGStdHepIdxZ_hldr, "kGStdHepIdxZ/I");
  configTree->Branch("kGStdHepIdxT", &kGStdHepIdxT_hldr, "kGStdHepIdxT/I");

  TRooTrackerVtxModuleBase::FillConfigTree(configTree);
}