LArStackingAction.cxx

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//

#include "larsim/LegacyLArG4/LArStackingAction.h"
#include "larcore/Geometry/Geometry.h"

#include "CLHEP/Units/SystemOfUnits.h"

#include "Geant4/G4MuonMinus.hh"
#include "Geant4/G4MuonPlus.hh"
#include "Geant4/G4ParticleDefinition.hh"
#include "Geant4/G4StackManager.hh"
#include "Geant4/G4String.hh"
#include "Geant4/G4ThreeVector.hh"
#include "Geant4/G4Track.hh"
#include "Geant4/G4TrackStatus.hh"
#include "Geant4/G4VProcess.hh"

// ROOT includes
#include "RtypesCore.h"
#include "TString.h"
#include "TVector3.h"

// Framework includes
#include "art/Framework/Services/Registry/ServiceHandle.h"

LArStackingAction::LArStackingAction(G4int dum)
  : fstage(0), freqMuon(2), freqIsoMuon(0), freqIso(10), fangRoI(30. * CLHEP::deg)
{
  fStack = dum;
  // Positive values effect action in this routine. Negative values effect action in
  // G4BadIdeaAction.
}

G4ClassificationOfNewTrack LArStackingAction::ClassifyNewTrack(const G4Track* aTrack)
{
  G4ClassificationOfNewTrack classification = fWaiting;
  art::ServiceHandle<geo::Geometry const> geom;
  TString volName(InsideTPC(aTrack));
  Double_t buffer = 500; // Keep muNucl neutrals within 5m (for now) of larTPC.

  // These 3 for now for investigation in gdb.
  int ppdg = aTrack->GetParentID();
  TString process("NA");
  if (ppdg) process = (TString)aTrack->GetCreatorProcess()->GetProcessName();

  switch (fstage) {
  case 0: // Fstage 0 : Primary muons only
    if (aTrack->GetParentID() == 0) {
      G4ParticleDefinition* particleType = aTrack->GetDefinition();
      if (((particleType == G4MuonPlus::MuonPlusDefinition()) ||
           (particleType == G4MuonMinus::MuonMinusDefinition())) &&
          !volName.Contains("unknown")) {
        classification = fUrgent;
      }
    }
    if (volName.Contains("unknown")) classification = fKill;
    break;

  case 1: // Stage 1 : K0,Lambda,n's made urgent here.
          //           Suspended tracks will be sent to the waiting stack
    if (aTrack->GetTrackStatus() == fSuspend) { break; }

    if ((aTrack->GetDefinition()->GetPDGEncoding() == 2112 ||
         aTrack->GetDefinition()->GetPDGEncoding() == 130 ||
         aTrack->GetDefinition()->GetPDGEncoding() == 310 ||
         aTrack->GetDefinition()->GetPDGEncoding() == 311 ||
         aTrack->GetDefinition()->GetPDGEncoding() == 3122) &&
        (aTrack->GetParentID() == 1) && !volName.Contains("unknown")) {

      const G4ThreeVector tr4Pos = aTrack->GetPosition();
      // G4 returns positions in mm, have to convert to cm for LArSoft coordinate systems
      const TVector3 trPos(tr4Pos.x() / CLHEP::cm, tr4Pos.y() / CLHEP::cm, tr4Pos.z() / CLHEP::cm);
      classification = fUrgent;

      auto const& tpc = geom->TPC({0, 0});
      if (trPos.X() < (tpc.HalfWidth() * 2.0 + buffer) && trPos.X() > (-buffer) &&
          trPos.Y() < (tpc.HalfHeight() * 2.0 + buffer) &&
          trPos.Y() > (-tpc.HalfHeight() * 2.0 - buffer) && trPos.Z() < (tpc.Length() + buffer) &&
          trPos.Z() > (-buffer)) {
        classification = fUrgent;
        break;
      }
      // These tracks need to be "scored" cuz every now and then they might get to the
      // LAr.
      else {
        classification = fKill;
        break;
      }
    }

    if (volName.Contains("unknown")) classification = fKill;
    break;

  default:
    // Track all other Primaries. Accept all secondaries in TPC.  Kill muon ionization
    // electrons outside TPC ignore primaries since they have no creator process
    if (aTrack->GetParentID() == 0 && !volName.Contains("unknown")) {
      classification = fUrgent;
      break;
    }

    auto const& tpc = geom->TPC({0, 0});
    if (volName.Contains(tpc.ActiveVolume()->GetName()) && aTrack->GetParentID() != 0) {
      classification = fUrgent;
      if (fStack & 0x4 && aTrack->GetCreatorProcess()->GetProcessName().contains("muIoni")) {
        classification = fKill;
      }
      break;
    }
    else if (volName.Contains("unknown")) {
      classification = fKill;
      break;
    }
    // Leave this here, even though I claim we've Killed these in stage 2.
    if (aTrack->GetDefinition()->GetPDGEncoding() == 11 &&
        aTrack->GetCreatorProcess()->GetProcessName().contains("muIoni")) {
      classification = fKill;
      break;
    }
    // For now, kill every other thing, no matter where it is.
    classification = fKill;

  } // end switch

  return classification;
}

std::string LArStackingAction::InsideTPC(const G4Track* aTrack)
{
  art::ServiceHandle<geo::Geometry const> geom;
  const G4ThreeVector tr4Pos = aTrack->GetPosition();

  // G4 returns positions in mm, have to convert to cm for LArSoft coordinate systems
  const geo::Point_t trPos(tr4Pos.x() / CLHEP::cm, tr4Pos.y() / CLHEP::cm, tr4Pos.z() / CLHEP::cm);

  return geom->VolumeName(trPos);
}

void LArStackingAction::NewStage()
{
  // Here when Urgent stack is empty. Waiting stack about to be made Urgent, upon saying
  // ReClassify().
  fstage++;

  // I yanked the ExN04's use here of stackManager->clear(), which clears stack and
  // prepares to end the event. Think I may wanna do something like this if muon has been
  // tracked, its doca is large, there are no hit voxels in the TPC, and I'm in further
  // need of optimization.

  if (fstage == 1) {
    // Stage 0->1 : check if at least "reqMuon" hits on muon chamber otherwise abort
    //              current event

    stackManager->ReClassify();
    return;
  }

  else if (fstage == 2) {
    // Stage 1->2 : check the isolation of muon tracks at least "reqIsoMuon" isolated
    //              muons otherwise abort current event.  Isolation requires "reqIso" or
    //              less hits (including own hits) in the RoI region in the tracker
    //              layers.
    stackManager->ReClassify();
    return;
  }

  else {
    // Other stage change : just re-classify
    stackManager->ReClassify();
  }
}

void LArStackingAction::PrepareNewEvent()
{
  fstage = 0;
}