TNRooTrackerVtxModule.cxx

Go to the documentation of this file.

#include "TNRooTrackerVtxModule.hxx"
#include "TG4PrimaryParticle.hxx"
#include "TG4PrimaryVertex.hxx"
#include "TGeomInfo.hxx"
#include "TND280Log.hxx"
#include "TOADatabase.hxx"
#include "TP0DGeom.hxx"

ND::TNRooTrackerVtxModule::TNRooTrackerVtxModule(const char *name,
                                                 const char *title) {
  fVtx = new TClonesArray("ND::NRooTrackerVtx", 100);
  SetNameTitle(name, title);
  // Enable this module by default:
  fIsEnabled = kTRUE;
  fDescription = "Neut-specific module which fills mc-truth vertex info";
  fCurrNeutVtx = new NRooTrackerVtx();
}

ND::TNRooTrackerVtxModule::~TNRooTrackerVtxModule() {}

Bool_t ND::TNRooTrackerVtxModule::IsKnownGenerator(const char *name) {
  std::string name2compare(name);
  bool is_neut1 = (name2compare.find("neut") != name2compare.npos);
  bool is_neut2 = (name2compare.find("nRooTracker") != name2compare.npos);
  bool is_neut = (is_neut1 || is_neut2);
  ND280Debug("Returning is generator = " << is_neut
                                         << " for generator name = " << name);
  return is_neut;
}

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

Bool_t ND::TNRooTrackerVtxModule::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::TNRooTrackerVtxModule::FillVtx(int entry) {
  fCurrNeutVtx->Reset();
  fRooTrackerTree->GetEntry(entry);

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

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

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

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

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

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

  /// Copy the current neut event. Make sure this is done last so we copy a
  /// complete event.
  fNRooTrackerVtx->Copy(fCurrNeutVtx);
}

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

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

  // NEUT > v5.0.7 && MCP > 1 (>10a)
  fRooTrackerTree->SetBranchAddress("NEnvc", &(fCurrNeutVtx->NEnvc));
  fRooTrackerTree->SetBranchAddress("NEipvc", (fCurrNeutVtx->NEipvcTemp));
  fRooTrackerTree->SetBranchAddress("NEpvc", (fCurrNeutVtx->NEpvc));
  fRooTrackerTree->SetBranchAddress("NEiorgvc", (fCurrNeutVtx->NEiorgvcTemp));
  fRooTrackerTree->SetBranchAddress("NEiflgvc", (fCurrNeutVtx->NEiflgvcTemp));
  fRooTrackerTree->SetBranchAddress("NEicrnvc", (fCurrNeutVtx->NEicrnvcTemp));

  fRooTrackerTree->SetBranchAddress("NEnvert", &(fCurrNeutVtx->NEnvert));
  fRooTrackerTree->SetBranchAddress("NEposvert", (fCurrNeutVtx->NEposvert));
  fRooTrackerTree->SetBranchAddress("NEiflgvert",
                                    (fCurrNeutVtx->NEiflgvertTemp));
  fRooTrackerTree->SetBranchAddress("NEnvcvert", &(fCurrNeutVtx->NEnvcvert));
  fRooTrackerTree->SetBranchAddress("NEdirvert", (fCurrNeutVtx->NEdirvert));
  fRooTrackerTree->SetBranchAddress("NEabspvert",
                                    (fCurrNeutVtx->NEabspvertTemp));
  fRooTrackerTree->SetBranchAddress("NEabstpvert",
                                    (fCurrNeutVtx->NEabstpvertTemp));
  fRooTrackerTree->SetBranchAddress("NEipvert", (fCurrNeutVtx->NEipvertTemp));
  fRooTrackerTree->SetBranchAddress("NEiverti", (fCurrNeutVtx->NEivertiTemp));
  fRooTrackerTree->SetBranchAddress("NEivertf", (fCurrNeutVtx->NEivertfTemp));
  // end NEUT > v5.0.7 && MCP > 1 (>10a)

  // NEUT > v5.1.2 && MCP >= 5
  fRooTrackerTree->SetBranchAddress("NEcrsx", &(fCurrNeutVtx->NEcrsx));
  fRooTrackerTree->SetBranchAddress("NEcrsy", &(fCurrNeutVtx->NEcrsy));
  fRooTrackerTree->SetBranchAddress("NEcrsz", &(fCurrNeutVtx->NEcrsz));
  fRooTrackerTree->SetBranchAddress("NEcrsphi", &(fCurrNeutVtx->NEcrsphi));

  // Nucleon FSI information
  fHaveNFBranches = !fRooTrackerTree->SetBranchAddress("NFnvert", &(fCurrNeutVtx->NFnvert));
  if(fHaveNFBranches){
    fRooTrackerTree->SetBranchAddress("NFiflag", (fCurrNeutVtx->NFiflagTEMP));
    fRooTrackerTree->SetBranchAddress("NFx", (fCurrNeutVtx->NFxTEMP));
    fRooTrackerTree->SetBranchAddress("NFy", (fCurrNeutVtx->NFyTEMP));
    fRooTrackerTree->SetBranchAddress("NFz", (fCurrNeutVtx->NFzTEMP));
    fRooTrackerTree->SetBranchAddress("NFpx", (fCurrNeutVtx->NFpxTEMP));
    fRooTrackerTree->SetBranchAddress("NFpy", (fCurrNeutVtx->NFpyTEMP));
    fRooTrackerTree->SetBranchAddress("NFpz", (fCurrNeutVtx->NFpzTEMP));
    fRooTrackerTree->SetBranchAddress("NFe", (fCurrNeutVtx->NFeTEMP));
    fRooTrackerTree->SetBranchAddress("PCascProb", &(fCurrNeutVtx->PCascProb));

    fRooTrackerTree->SetBranchAddress("NFfirststep",
                                      (fCurrNeutVtx->NFfirststepTEMP));
    fRooTrackerTree->SetBranchAddress("NFnstep", &(fCurrNeutVtx->NFnstep));
    fRooTrackerTree->SetBranchAddress("NFecms2", (fCurrNeutVtx->NFecms2TEMP));
    fRooTrackerTree->SetBranchAddress("Prob", (fCurrNeutVtx->ProbTEMP));
    fRooTrackerTree->SetBranchAddress("VertFlagStep", (fCurrNeutVtx->VertFlagStepTEMP));
    fRooTrackerTree->SetBranchAddress("VertFsiRhon", (fCurrNeutVtx->VertFsiRhonTEMP));
    fRooTrackerTree->SetBranchAddress("StepPel", (fCurrNeutVtx->StepPelTEMP));
    fRooTrackerTree->SetBranchAddress("StepPsp", (fCurrNeutVtx->StepPspTEMP));
    fRooTrackerTree->SetBranchAddress("StepPdp", (fCurrNeutVtx->StepPdpTEMP));
  } else {
    ND280NamedLog("TNRooTrackerVtxModule", "Passthrough tree doesn't contain"
      " nucleon FSI tracking information, these branches will remain off.");
  }
  fRooTrackerTree->SetBranchAddress("SPIDelta", &(fCurrNeutVtx->SPIDelta));
  fRooTrackerTree->SetBranchAddress("IRadCorrPht", &(fCurrNeutVtx->IRadCorrPht));

}

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

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

  int kNStdHepNPmax_hldr = kNStdHepNPmax;
  int kNStdHepIdxPx_hldr = kNStdHepIdxPx;
  int kNStdHepIdxPy_hldr = kNStdHepIdxPy;
  int kNStdHepIdxPz_hldr = kNStdHepIdxPz;
  int kNStdHepIdxE_hldr = kNStdHepIdxE;
  int kNStdHepIdxX_hldr = kNStdHepIdxX;
  int kNStdHepIdxY_hldr = kNStdHepIdxY;
  int kNStdHepIdxZ_hldr = kNStdHepIdxZ;
  int kNStdHepIdxT_hldr = kNStdHepIdxT;

  int kNEmaxvc_hldr = kNEmaxvc;
  int kNEmaxvert_hldr = kNEmaxvert;
  int kNEmaxvertp_hldr = kNEmaxvertp;

  configTree->Branch("kNStdHepNPmax", &kNStdHepNPmax_hldr, "kNStdHepNPmax/I");
  configTree->Branch("kNStdHepIdxPx", &kNStdHepIdxPx_hldr, "kNStdHepIdxPx/I");
  configTree->Branch("kNStdHepIdxPy", &kNStdHepIdxPy_hldr, "kNStdHepIdxPy/I");
  configTree->Branch("kNStdHepIdxPz", &kNStdHepIdxPz_hldr, "kNStdHepIdxPz/I");
  configTree->Branch("kNStdHepIdxE", &kNStdHepIdxE_hldr, "kNStdHepIdxE/I");
  configTree->Branch("kNStdHepIdxX", &kNStdHepIdxX_hldr, "kNStdHepIdxX/I");
  configTree->Branch("kNStdHepIdxY", &kNStdHepIdxY_hldr, "kNStdHepIdxY/I");
  configTree->Branch("kNStdHepIdxZ", &kNStdHepIdxZ_hldr, "kNStdHepIdxZ/I");
  configTree->Branch("kNStdHepIdxT", &kNStdHepIdxT_hldr, "kNStdHepIdxT/I");
  configTree->Branch("kNEmaxvc", &kNEmaxvc_hldr, "kNEmaxvc/I");
  configTree->Branch("kNEmaxvert", &kNEmaxvert_hldr, "kNEmaxvert/I");
  configTree->Branch("kNEmaxvertp", &kNEmaxvertp_hldr, "kNEmaxvertp/I");

  TRooTrackerVtxModuleBase::FillConfigTree(configTree);
}