diff --git a/FPFSim.cc b/FPFSim.cc index 1910f7f..69ae1eb 100644 --- a/FPFSim.cc +++ b/FPFSim.cc @@ -17,10 +17,11 @@ using namespace std; * - run macros * - start interactive UI mode (no arguments) */ -int main(int argc, char** argv) { - std::cout<<"Application starting..."<ApplyCommand("/control/execute macros/vis.mac"); ui->SessionStart(); @@ -61,7 +62,7 @@ int main(int argc, char** argv) { ui->SessionStart(); delete ui; } else { - std::cout<<"Please specify the second argument as vis to visualize the event"<> temp.mac /gen/genie/genieInput nue_kling_ar40_e5000.ghep.root /gen/genie/genieIStart ${istart} -/histo/addDiffusion false -/histo/saveHit false -/histo/save3DEvd true -/histo/save2DEvd false -/histo/circleFit true -/histo/fileName ${outputfile} +/out/flare/addDiffusion false + +/out/flare/save3DEvd true +/out/flare/save2DEvd false + +/out/fileName ${outputfile} /run/beamOn ${num_evt_per_file} EOF diff --git a/GridUtils/dune-gpvms/numu/README.md b/GridUtils/dune-gpvms/numu/README.md index 2de2a71..2cad557 100644 --- a/GridUtils/dune-gpvms/numu/README.md +++ b/GridUtils/dune-gpvms/numu/README.md @@ -1,27 +1,32 @@ # Description This directory contains scripts for submitting batch jobs to simulate numu interactions in FLArE. -* `batch_script.sh` is the script that is actually being submitted and executed on the worker nodes. - Worker nodes have a limited scope, so the script takes care of: - 1. copying over from `pnfs` all inputs (including the executable); - 2. running the executable; - 3. copying back to `pnfs` output and log files. -* `submit.sh` is the script that handles the job submission by calling `jobsub_submit`. - It sets the required computing resources and passes the input filepaths to `batch_script.sh`. +- `batch_script.sh` is the script that is actually being submitted and executed on the worker nodes. + Worker nodes have a limited scope, so the script takes care of: -* `numu_job.sh` is the main script: - 1. it sets production parameters (number of files, number of events) and paths; - 2. it generates several `.mac` files starting from a template, keeping track of them in a list; - 3. it passes the lists of configuration files and other paths to `submit.sh`. + 1. copying over from `pnfs` all inputs (including the executable); + 2. running the executable; + 3. copying back to `pnfs` output and log files. -* `glob_output.sh` is the script that retrieves all the job output files in `/pnfs/dune/scratch`. - The files are copied over to `dune\data` for storage. +- `submit.sh` is the script that handles the job submission by calling `jobsub_submit`. + It sets the required computing resources and passes the input filepaths to `batch_script.sh`. -* `merge_outputs.sh` is the script that merges all the output files into a single one. +- `numu_job.sh` is the main script: + + 1. it sets production parameters (number of files, number of events) and paths; + 2. it generates several `.mac` files starting from a template, keeping track of them in a list; + 3. it passes the lists of configuration files and other paths to `submit.sh`. + +- `glob_output.sh` is the script that retrieves all the job output files in `/pnfs/dune/scratch`. + The files are copied over to `dune\data` for storage. + +- `merge_outputs.sh` is the script that merges all the output files into a single one. ## Instructions: -* The `.mac` template is + +- The `.mac` template is + ``` /control/execute ${geomacro} @@ -32,19 +37,21 @@ This directory contains scripts for submitting batch jobs to simulate numu inter /genie/genieInput numu_kling_ar40_e5000.ghep.root /genie/genieIStart ${istart} -/histo/saveEvd false -/histo/saveHit false -/histo/circleFit true -/histo/fileName ${outputfile} +/out/saveEvd false + + +/out/fileName ${outputfile} /run/beamOn ${num_evt_per_file} -``` +``` + The script relies on GENIE-simulated numu neutrino vertices stored in `numu_kling_ar40_e5000.ghep.root`. In order for this file to be visibile on the worker nodes, `batch_script.sh` copies it over from an hardcoded `pnfs` location: + ``` ifdh cp /pnfs/dune/persistent/users/mvicenzi/numu/numu_kling_ar40_e5000.ghep.root numu_kling_ar40_e5000.ghep.root ``` -* The current implementation requires the user to provide the name of the geometry macro from `macro\geometry_options` by writing it in `numu_job.sh`. The name of this file is used to create the output directory as well as for the filenames themselves. It is the only argument that is needed by `glob_output.sh` and `merge_outputs.sh`. -* `${istart}` is set according to the number of files requested and the number of events per file chosen, in order to avoid using the same events. -* Most of the paths are currently hardcoded to point to `/pnfs/dune/scratch/users/mvicenzi` to store (temporary) submission files as well as the output files. They can be changed as necessary, but they need to be on `pnfs` to be accessible by the worker nodes with `ifdh cp`. - + +- The current implementation requires the user to provide the name of the geometry macro from `macro\geometry_options` by writing it in `numu_job.sh`. The name of this file is used to create the output directory as well as for the filenames themselves. It is the only argument that is needed by `glob_output.sh` and `merge_outputs.sh`. +- `${istart}` is set according to the number of files requested and the number of events per file chosen, in order to avoid using the same events. +- Most of the paths are currently hardcoded to point to `/pnfs/dune/scratch/users/mvicenzi` to store (temporary) submission files as well as the output files. They can be changed as necessary, but they need to be on `pnfs` to be accessible by the worker nodes with `ifdh cp`. diff --git a/GridUtils/dune-gpvms/numu/numu_job.sh b/GridUtils/dune-gpvms/numu/numu_job.sh index 6ffedf4..add4042 100644 --- a/GridUtils/dune-gpvms/numu/numu_job.sh +++ b/GridUtils/dune-gpvms/numu/numu_job.sh @@ -33,9 +33,9 @@ cat << EOF >> temp.mac /gen/genie/genieInput numu_kling_ar40_e5000.ghep.root /gen/genie/genieIStart ${istart} -/histo/saveHit false -/histo/circleFit true -/histo/fileName ${outputfile} + + +/out/fileName ${outputfile} /run/beamOn ${num_evt_per_file} EOF diff --git a/GridUtils/dune-gpvms/nutau/nutau_job.sh b/GridUtils/dune-gpvms/nutau/nutau_job.sh index 0972e3c..08276e1 100644 --- a/GridUtils/dune-gpvms/nutau/nutau_job.sh +++ b/GridUtils/dune-gpvms/nutau/nutau_job.sh @@ -33,12 +33,12 @@ cat << EOF >> temp.mac /gen/genie/genieInput nutau_bai_ar40_e2000.ghep.root /gen/genie/genieIStart ${istart} -/histo/addDiffusion false -/histo/saveHit false -/histo/save3DEvd true -/histo/save2DEvd false -/histo/circleFit true -/histo/fileName ${outputfile} +/out/flare/addDiffusion false + +/out/flare/save3DEvd true +/out/flare/save2DEvd false + +/out/fileName ${outputfile} /run/beamOn ${num_evt_per_file} EOF diff --git a/GridUtils/dune-gpvms/single_particle/gen_jobfile.sh b/GridUtils/dune-gpvms/single_particle/gen_jobfile.sh index 1a8ac19..89428ee 100644 --- a/GridUtils/dune-gpvms/single_particle/gen_jobfile.sh +++ b/GridUtils/dune-gpvms/single_particle/gen_jobfile.sh @@ -30,8 +30,8 @@ cat << EOF >> ${filename} /gps/particle ${particle} /gps/ene/mono ${particle_kin} GeV -/histo/saveEvd false -/histo/fileName ${outputfile} +/out/flare/save2DEvd false +/out/fileName ${outputfile} /run/beamOn ${num_evt_per_file} EOF diff --git a/GridUtils/dune-gpvms/single_particle/gen_muons.sh b/GridUtils/dune-gpvms/single_particle/gen_muons.sh index 896d818..4f056e1 100644 --- a/GridUtils/dune-gpvms/single_particle/gen_muons.sh +++ b/GridUtils/dune-gpvms/single_particle/gen_muons.sh @@ -35,10 +35,10 @@ cat << EOF >> ${filename} /gps/particle ${particle} /gps/ene/mono ${particle_kin} GeV -/histo/saveEvd false -/histo/saveHit false -/histo/circleFit true -/histo/fileName ${outputfile} +/out/flare/save2DEvd false + + +/out/fileName ${outputfile} /run/beamOn ${num_evt_per_file} EOF diff --git a/GridUtils/lxplus/nujobs/nu_job.sh b/GridUtils/lxplus/nujobs/nu_job.sh index 6235dcb..e624631 100755 --- a/GridUtils/lxplus/nujobs/nu_job.sh +++ b/GridUtils/lxplus/nujobs/nu_job.sh @@ -81,12 +81,12 @@ function generate_macros { /gen/genie/genieInput ${gst} /gen/genie/genieIStart ${istart} -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/addDiffusion false -/histo/circleFit true -/histo/fileName ${outputfile} +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/flare/addDiffusion false + +/out/fileName ${outputfile} /run/beamOn ${n_evt_per_job} EOF diff --git a/GridUtils/sdcc/bkgjobs/batch_bkg_script.sh b/GridUtils/sdcc/bkgjobs/batch_bkg_script.sh new file mode 100755 index 0000000..61ae870 --- /dev/null +++ b/GridUtils/sdcc/bkgjobs/batch_bkg_script.sh @@ -0,0 +1,26 @@ +#!/bin/bash + +cluster=$1 +process=$2 +fpfsim=$3 +macrolist=$4 +setup=$5 + +echo "Executing JOBID ${cluster}.${process}" + +# source the environment +source $setup + +# select the macro file +let num=$((${process}+1)) +echo "Selecting macro from line ${num} in list:" + +macropath=$(tail -n+${num} ${macrolist} | head -n1) +macro=`basename "$macropath"` +echo "$macro" + +# running !! +echo "Running ${fpfsim} ${macropath}" +${fpfsim} ${macropath} + +echo "Completed JOBID ${cluster}.${process}" diff --git a/GridUtils/sdcc/bkgjobs/bkg_job.sh b/GridUtils/sdcc/bkgjobs/bkg_job.sh new file mode 100755 index 0000000..0a79176 --- /dev/null +++ b/GridUtils/sdcc/bkgjobs/bkg_job.sh @@ -0,0 +1,176 @@ +#!/bin/bash + +###----------------------------------------- +### DEFINITION OF PRODUCTION PARAMETERS +###---------------------------------------- + +# Production name for output directories +# This is used to place output logs and files +export prodname="test_campaign" + +# Define how many jobs, how many files +# Jobs will be placed in the same cluster +export n_jobs=10 +export n_evt_per_job=200 + +# Paths to FPFsim build directory +export builddir="/direct/lbne+u/${USER}/FPFSim/build" +export fpfsim="${builddir}/FPFSim" +export setup="${builddir}/../sdcc_setup.sh" +export libdict="${builddir}/libFPFClasses_rdict.pcm" + +# Generator options +export twindow="200 us" +export inputbkg="${builddir}/../backgrounds/background_input.root" + +# Path to geometry macro +export geometry="${builddir}/macros/geometry_options/FLArE_only_BabyMIND.mac" +export rock="false" + +# Path to log/output directory +export outdir="/gpfs01/lbne/users/fpf/$USER/CONDOR_OUTPUT" +export logdir="/gpfs01/lbne/users/fpf/$USER/CONDOR_LOGS" + +###------------------------------------------------------------------ +###------------------------------------------------------------------ + +function generate_macros { + + list=$1 + + # for each job + for i in $(seq 0 $((${n_jobs}-1))) + do + + # setup different seeds per jobs + seed1=$((42+i)) + seed2=$((47+i)) + + # select the event numbering offset for this job + istart=$((n_evt_per_job*i)) + + # select name for output file + outputfile="${prodname}_${i}.root" + + # final macro path + macfile="${logdir}/${prodname}/mac/${prodname}_${i}.mac" + if [ -f ${macfile} ]; then + echo "${macfile} exists, delete file first" + rm ${macfile} + fi + + # add macro path to list + echo $macfile >> $list + rm temp.mac + cat << EOF >> temp.mac +/control/execute ${geometry} +/det/enableRock ${rock} +/det/flare/useNativeG4Scorer true + +/random/setSeeds ${seed1} ${seed2} +/run/initialize + +/gen/select background +/gen/bkg/backgroundInput ${inputbkg} +/gen/bkg/backgroundWindow ${twindow} +/gen/bkg/eventOffset ${istart} + +/out/flare/save3DEvd false +/out/flare/save2DEvd false +/out/flare/addDiffusion false +/out/fileName ${outputfile} + +/run/beamOn ${n_evt_per_job} +EOF + + # move the macro to final destination + cp temp.mac ${macfile} + done + +} + +###------------------------------------------------------------------ +###------------------------------------------------------------------ + +function generate_submission_file { + + sub=$1; shift + listpath=$1; shift + + cat << EOF >> ${sub} +universe = vanilla +notification = never +executable = batch_bkg_script.sh +arguments = \$(ClusterId) \$(ProcId) ${fpfsim} ${listpath} ${setup} +initialdir = ${outdir}/${prodname} +output = ${logdir}/${prodname}/out/\$(ClusterId).\$(ProcId).out +error = ${logdir}/${prodname}/err/\$(ClusterId).\$(ProcId).err +log = ${logdir}/${prodname}/log/\$(ClusterId).\$(ProcId).log +getenv = True +request_memory = 8000 +queue ${n_jobs} +EOF + +} + +###------------------------------------------------------------------ +###------------------------------------------------------------------ + +# create log directories +echo "Creating log directories for: '$prodname'" +mkdir -p ${logdir}/${prodname}/out +mkdir -p ${logdir}/${prodname}/err +mkdir -p ${logdir}/${prodname}/log +mkdir -p ${logdir}/${prodname}/mac + +# check if production with the same name already exists +if [ -d "${outdir}/${prodname}" ] && [ "$(ls -A ${outdir}/${prodname})" ]; then + echo "WARNING: ${outdir}/${prodname} is not empty!" + echo "Outputs for production ${prodname} already exist, delete them first" + echo "or change the name for this production!" + exit 1 +else + mkdir -p ${outdir}/${prodname} + if [ -d "${logdir}/${prodname}/out/" ] && [ "$(ls -A ${logdir}/${prodname}/out/)" ]; + then rm ${logdir}/${prodname}/out/*; fi + if [ -d "${logdir}/${prodname}/err/" ] && [ "$(ls -A ${logdir}/${prodname}/err/)" ]; + then rm ${logdir}/${prodname}/err/*; fi + if [ -d "${logdir}/${prodname}/log/" ] && [ "$(ls -A ${logdir}/${prodname}/log/)" ]; + then rm ${logdir}/${prodname}/log/*; fi + if [ -d "${logdir}/${prodname}/mac/" ] && [ "$(ls -A ${logdir}/${prodname}/mac/)" ]; + then rm ${logdir}/${prodname}/mac/*; fi +fi + + +# generate all macros for each job +# as well as a list of their paths +echo "Generating list of .mac files..." + +macrolist="macrolist.txt" +if [ -f ${macrolist} ]; then + rm ${macrolist} +fi +touch $macrolist + +generate_macros $macrolist + +macropath="${logdir}/${prodname}/${macrolist}" +cp $macrolist ${macropath} + +echo "Generating .sub file..." + +subfile=${prodname}.sub +if [ -f ${subfile} ]; then + rm ${subfile} +fi +touch ${subfile} + +generate_submission_file $subfile $macropath + +echo "Bookeeping submission files..." + +cp $subfile ${logdir}/${prodname} + +condor_submit $subfile + +rm ${macrolist} ${subfile} temp.mac diff --git a/GridUtils/sdcc/bkgjobs/merge_outputs.sh b/GridUtils/sdcc/bkgjobs/merge_outputs.sh new file mode 100644 index 0000000..8ade22a --- /dev/null +++ b/GridUtils/sdcc/bkgjobs/merge_outputs.sh @@ -0,0 +1,23 @@ +#!/bin/bash + +export prodname=$1 + +export dir="/gpfs01/lbne/users/fpf/${USER}/CONDOR_OUTPUT/${prodname}/${prodname}_*.root" +export out="/gpfs01/lbne/users/fpf/${USER}/CONDOR_OUTPUT/${prodname}/${prodname}.root" +list="" + +for f in $(ls ${dir}); +do + #echo $f + list+=" ${f}" +done + +#echo $list +echo "Merging into ${out}" + +if test -f "$out"; then + echo "$out exists. Removing old file" + rm $out +fi + +hadd ${out} ${list} diff --git a/README.md b/README.md index 2e1b973..1321443 100644 --- a/README.md +++ b/README.md @@ -181,18 +181,18 @@ Older versions of FORESEE output events in the HepMC2 format. To run over HepMC2 |/gen/hepmc/vtxOffset | if `hepmc` is selected, give an x, y, z offset to each vertex. Useful if you need to align events from another generator (i.e. FORESEE) with this G4 simulation.| |/gen/hepmc/placeInDecayVolume | if `hepmc` is selected, will automatically translate vertex into the FASER2 decay volume. Assumes that the vertices start from (0,0,0). If this is not the case then `/gen/hepmc/vtxOffset` must also be set.| -### Analysis +### Output |Command |Description | |:--|:--| -|/histo/fileName | option for AnalysisManagerMessenger, set name of the file saving all analysis variables| -|/histo/saveHit | if `true` save info for all hits, `false` in default to save space| -|/histo/saveTrack | if `true` save info for all tracks, `false` in default, requires `\tracking\storeTrajectory 1`| -|/histo/save3DEvd | if `true` save 3D spatial distribution of energy deposition, `false` in default| -|/histo/save2DEvd | if `true` save 2D spatial distribution of energy deposition, `false` in default| -|/histo/circleFit | if `true` run circle fitting and save information in output, `false` in default to save space| -|/histo/addDiffusion | if `toy` diffuse energy, if `single` diffuse single electron, `false` in default without diffusion| -|/histo/actsHits | if `true` save the truth particle and hit information in a format that can be read by the FASER2 Acts tracking software, `true` by default if FASER2 is enabled| +|/out/fileName | option for AnalysisManagerMessenger, set name of the file saving all analysis variables| +|/out/saveTrack | if `true` save all tracks, `false` by default, requires `\tracking\storeTrajectory 1`| +|/out/flare/enableOutput | if `true` save FLArE output, `true` by default | +|/out/flare/save3DEvd | if `true` save 3D spatial distribution of energy deposition, `false` by default | +|/out/flare/save2DEvd | if `true` save 2D spatial distribution of energy deposition, `false` by default | +|/out/flare/addDiffusion | if `toy` diffuse energy, if `single` diffuse single electron, `false` by default | +|/out/flare/pseudoReco | if `true` save pseudo-reco ouput tree, `false` by default | +|/out/faser/actsHits | if `true` save the truth particle and hit information in a format that can be read by the FASER2 Acts tracking software, `true` by default if FASER2 is enabled| ## How to save an event display with high resolution using the DAWN Event Display diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index dd42ae9..cbe7829 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -2,11 +2,15 @@ #define ANALYSISMANAGER_HH #include +#include +#include + #include "globals.hh" -#include -#include -#include -#include +#include "G4Event.hh" +#include "TFile.h" +#include "TTree.h" +#include "TH2F.h" + #include "AnalysisManagerMessenger.hh" #include "PixelMap3D.hh" #include "FPFParticle.hh" @@ -16,179 +20,203 @@ class AnalysisManager { public: + AnalysisManager(); ~AnalysisManager(); static AnalysisManager* GetInstance(); - void bookEvtTree(); - void bookTrkTree(); - void BeginOfRun(); + + //------------------------------------------------ + // Functions executed at specific times + void BeginOfRun(); void EndOfRun(); void BeginOfEvent(); void EndOfEvent(const G4Event* event); - TFile* GetOutputFile() { return thefile; } - void SetTrackPTPair(G4int PID, G4int TID) { allTracksPTPair.insert(std::make_pair(PID, TID)); } - void AddOnePrimaryTrack() { nTestNPrimaryTrack++; } - public: - // function for controlling from the configuration file - void setFileName(std::string val) { m_filename = val; } - void saveHit(G4bool val) { m_saveHit = val; } - void saveTrack(G4bool val) { m_saveTrack = val; } - void save3DEvd(G4bool val) { m_save3DEvd = val; } - void save2DEvd(G4bool val) { m_save2DEvd = val; } - void circleFit(G4bool val) { m_circularFit = val; } - void saveActs(G4bool val) { m_saveActs = val; } - void addDiffusion(G4String val) { m_addDiffusion = val; } + //------------------------------------------------ + // functions for controlling from the configuration file + void setFileName(std::string val) { fFilename = val; } + void saveTrack(G4bool val) { fSaveTrack = val; } + void saveActs(G4bool val) { fSaveActs = val; } + void savePseudoReco(G4bool val) { fSavePseudoReco = val; } + void addDiffusion(G4String val) { fAddDiffusion = val; } + void save3DEvd(G4bool val) { fSave3DEvd = val; } + void save2DEvd(G4bool val) { fSave2DEvd = val; } + void enableFLArE(G4bool val) { fEnableFLArE = val;} + + // build TID to primary ancestor association + // filled progressively from StackingAction + void SetTrackPrimaryAncestor(G4int trackID, G4int ancestorID) { trackToPrimaryAncestor[trackID] = ancestorID; } + G4int GetTrackPrimaryAncestor(G4int trackID) { return trackToPrimaryAncestor.at(trackID); } + + // TODO: needed??? + void AddOnePrimaryTrack() { nTestNPrimaryTrack++; } private: - static AnalysisManager* instance; - AnalysisManagerMessenger* messenger; - - TFile* thefile; - std::string m_filename; - TTree* evt; - TTree* trk; - TTree* acts_hits_tree; - TTree* acts_particles_tree; - std::string fH5Filename; - hep_hpc::hdf5::File fH5file; - G4int evtID; - FPFNeutrino neutrino; + //------------------------------------------------ + // Book ROOT output TTrees + // common + detector specific + void bookEvtTree(); + void bookTrkTree(); + void bookPrimTree(); + void bookFLArETrees(); + void bookFASER2Trees(); + + void FillPrimariesTree(const G4Event* event); + void FillTrajectoriesTree(const G4Event* event); + + void FillFLArEOutput(); + void FillFLArEPseudoReco(); + + void FillFASER2Output(); + + float_t GetTotalEnergy(float_t px, float_t py, float_t pz, float_t m); + + static AnalysisManager* fInstance; + AnalysisManagerMessenger* fMessenger; + + G4bool fSaveTrack; + G4bool fSave3DEvd; + G4bool fSave2DEvd; + G4bool fSavePseudoReco; + TString fAddDiffusion; + G4bool fSaveActs; + G4bool fEnableFLArE; + + std::map fSDNamelist; + std::vector fFlareSDs; + std::vector fFaser2SDs; + + G4HCofThisEvent* fHCofEvent; + + G4int nPrimaryVertex; std::vector primaries; std::vector primaryIDs; - // Truth information from genie - G4int nuIdx; ///<- neutrino index (for genie neutrino interaction) - G4int nuPDG; ///<- neutrino PDG code (for genie neutrino interaction) - G4double nuE; ///<- neutrino energy - G4double nuX; ///<- neutrino vertex X - G4double nuY; ///<- neutrino vertex Y - G4double nuZ; ///<- neutrino vertex Z - G4int nuIntType; ///<- interaction type: CC, NC, et.al. - G4int nuScatteringType; ///<- scattering type: QE, DIS, RES, MEC, et. al. - G4double nuW; ///<- invariant hadronic mass - G4int nuFSLPDG; ///<- Final state lepton PDG code (for genie neutrino interaction) - G4double nuFSLPx; ///<- Final state lepton Px - G4double nuFSLPy; ///<- Final state lepton Py - G4double nuFSLPz; ///<- Final state lepton Pz - G4double nuFSLE; ///<- Final state lepton total energy (GeV) - - G4int nTestNPrimaryTrack; - G4int countPrimaryParticle; - G4int nPrimaryVertex; - G4int nPrimaryParticle; ///<- number of primary particle - ///<- (in case of genie neutrino interaction, number of stable particle in the final state) - ///<- (in case of the FSL decay, decay products counted as primary particle) - ///<- (in case of the final state pizero, decay products counted as primary particle) - //// Geant4 truth - G4int primaryParentPDG[1000]; ///<- parent PDG of primary particles - G4double primaryTrackLength[1000]; ///<- track length of primary particles - G4double primaryTrackLengthInTPC[1000]; ///<- track length of primary particles in TPC region - - // pseudo-reco - G4double ProngEInDetector[1000]; - G4double ProngEInLAr[1000]; - G4double ProngEInHadCal[1000]; - G4double ProngEInMuonFinder[1000]; - G4double ProngEInMuonFinderLayer1X[1000]; - G4double ProngEInMuonFinderLayer1Y[1000]; - G4double ProngEInMuonFinderLayer2X[1000]; - G4double ProngEInMuonFinderLayer2Y[1000]; - G4double ProngAngleToBeamDir[1000]; - G4double ShowerLength[1000]; - G4double ShowerLengthInLAr[1000]; - G4double ShowerWidth[1000]; - G4double ShowerWidthInLAr[1000]; - G4double ProngAvgdEdx[1000]; - G4double ProngAvgdEdxInLAr[1000]; - G4double ProngdEdxAlongTrack[1000][100]; - G4int ProngdEdxTrackLength[1000][100]; - G4double TotalDedxLongitudinal[3000]; - G4double TrueTotalDedxLongitudinal[3000]; - // reco - // direction - G4double dir_pol_x[1000]; - G4double dir_pol_y[1000]; - G4double dir_pol_z[1000]; - G4double dir_coc_x[1000]; - G4double dir_coc_y[1000]; - G4double dir_coc_z[1000]; - - G4double edepInLAr; - G4double edepInHadCalX; - G4double edepInHadCalY; - G4double edepInMuonFinderX; - G4double edepInMuonFinderY; - G4double edepInHadAborb; - G4double edepInMuonFinderAbsorb; - G4double missCountedEnergy; - - G4int nFromFSLParticles; - G4int nFromFSPizeroParticles; - G4int nFromFSLDecayPizeroParticles; - G4int fromFSLParticlePDG[2000000]; - - G4int nHits; - G4int HitTID[40000000]; - G4int HitPID[40000000]; - G4int HitPDG[40000000]; - G4int HitTrackStatus[40000000]; - G4double HitPrePositionX[40000000]; - G4double HitPrePositionY[40000000]; - G4double HitPrePositionZ[40000000]; - G4double HitPosPositionX[40000000]; - G4double HitPosPositionY[40000000]; - G4double HitPosPositionZ[40000000]; - G4double HitEdep[40000000]; - - G4bool m_saveHit; - G4bool m_saveTrack; - G4bool m_save3DEvd; - G4bool m_save2DEvd; - G4bool m_circularFit; - G4bool m_saveActs; - TString m_addDiffusion; - PixelMap3D* pm3D; - G4double sparseFractionMem; - G4double sparseFractionBins; - - // Circular fit in HadCat + MF - G4int circNhits; //number of hits - G4double xc, zc, rc; //circle fit - G4double p0, p1, cosDip; //dip angle fit - std::vector hitXFSL; // MC truth - std::vector hitZFSL; - std::vector hitYFSL; - std::vector hitPFSL; - - // Circular fit in FASER magnet(s) - G4int Nmagnets; // size of output vectors - std::vector magzpos; - G4int trkNhits; // hits in tracking stations - std::vector trkxc; //circle fit results - std::vector trkzc; - std::vector trkrc; - std::vector trkmIn; // entering track fit - std::vector trkqIn; - std::vector trkmOut; // exiting track fit - std::vector trkqOut; - G4double trkp0, trkp1, trkcosDip; //dip angle fit - std::vector trkXFSL; // MC truth - std::vector trkZFSL; - std::vector trkYFSL; - std::vector trkPFSL; - - // track information - int trackTID; - int trackPID; + //------------------------------------------------ + // output files and trees + std::string fFilename; + std::string fH5Filename; + hep_hpc::hdf5::File fH5file; + TFile* fFile; + TTree* fEvt; + TTree* fTrk; + TTree* fPrim; + + TDirectory* fFLArEDir; + TTree* fFLArEHits; + TTree* fFLArEHCALHits; + TTree* fFLArEPseudoReco; + + TDirectory* fFASER2Dir; + TTree* fActsHitsTree; + TTree* fActsParticlesTree; + + // track to primary ancestor + std::map trackToPrimaryAncestor; + + // TODO: no longer needed? + G4int nTestNPrimaryTrack; + + //--------------------------------------------------- + // OUTPUT VARIABLES FOR COMMON TREES + // TODO: review carefully + // TODO: need to make evt tree less FLARE-centric + // TODO: remove pseudo-reco or add it as FLARE-only tree + // TODO: turn arrays in std::vector! + + G4int evtID; + FPFNeutrino neutrino; //TODO: remove?? + + //--------------------------------------------------- + // Output variables for TRAJECTORIES tree + int trackTID; + int trackPID; int trackPDG; double trackKinE; - int trackNPoints; - std::vector trackPointX; - std::vector trackPointY; + int trackNPoints; + std::vector trackPointX; + std::vector trackPointY; std::vector trackPointZ; + //--------------------------------------------------- + // Output variables for PRIMARIES tree + UInt_t primVtxID; + UInt_t primParticleID; + UInt_t primTrackID; + UInt_t primPDG; // why unsigned? + float_t primM; + float_t primQ; + float_t primEta; + float_t primPhi; + float_t primPt; + float_t primP; + float_t primVx; + float_t primVy; + float_t primVz; + float_t primVt; + float_t primPx; + float_t primPy; + float_t primPz; + float_t primE; + float_t primKE; + + //--------------------------------------------------- + // OUTPUT VARIABLES FOR FLArE TREES + // TODO: merge hit variables? no need to use different names? + // TODO: somehow need to add back here info from evt tree + + PixelMap3D* pm3D; + + UInt_t flareTrackID; + UInt_t flareParticleID; + UInt_t flareParentID; + UInt_t flarePDG; + UInt_t flareCopyNum; + float_t flareT; + float_t flareX; + float_t flareY; + float_t flareZ; + float_t flarePx; + float_t flarePy; + float_t flarePz; + float_t flareDeltaPx; + float_t flareDeltaPy; + float_t flareDeltaPz; + float_t flareEdep; + bool flareIsZX; + + // flare pseudo-reco + float_t edepInLAr; + float_t edepInHCALX; + float_t edepInHCALY; + float_t sparseFractionMem; + float_t sparseFractionBins; + std::vector TotalDedxLongitudinal; + std::vector TrueTotalDedxLongitudinal; + G4int nprimaries; + std::vector primaryPDG; + std::vector primaryTrackLength; + std::vector primaryTrackLengthInTPC; + std::vector ProngEInLAr; + std::vector ProngEInHadCal; + std::vector ProngAngleToBeamDir; + std::vector ShowerLength; + std::vector ShowerLengthInLAr; + std::vector ShowerWidth; + std::vector ShowerWidthInLAr; + std::vector ProngAvgdEdx; + std::vector ProngAvgdEdxInLAr; + std::vector dir_pol_x; + std::vector dir_pol_y; + std::vector dir_pol_z; + std::vector dir_coc_x; + std::vector dir_coc_y; + std::vector dir_coc_z; + + //--------------------------------------------------- + // OUTPUT VARIABLES FOR FASER2 TREES + // Acts Hit Information - the types are set to match the types expected by Acts::RootSimHitReader UInt_t ActsHitsEventID; ULong64_t ActsHitsGeometryID; @@ -232,7 +260,7 @@ class AnalysisManager { std::vector ActsParticlesVertexPrimary; std::vector ActsParticlesVertexSecondary; std::vector ActsParticlesParticle; - + std::vector ActsParticlesGeneration; std::vector ActsParticlesSubParticle; std::vector ActsParticlesELoss; @@ -240,25 +268,7 @@ class AnalysisManager { std::vector ActsParticlesPathInL0; std::vector ActsParticlesNumberOfHits; std::vector ActsParticlesOutcome; - - private: - void FillPrimaryTruthTree(G4int sdId, std::string sdName); - void FillTrueEdep(G4int sdId, std::string sdName); - double GetTotalEnergy(double px, double py, double pz, double m); - void FillPseudoRecoVar(); - - G4int NumberOfSDs; - std::set > SDNamelist; - G4HCofThisEvent* hcofEvent; - - std::set > allTracksPTPair; - std::vector > trackClusters; - std::set tracksFromFSL; - std::set tracksFromFSLSecondary; - std::set tracksFromFSPizeroSecondary; - std::set tracksFromFSLDecayPizeroSecondary; - int fPrimIdxFSL; }; #endif diff --git a/include/AnalysisManagerMessenger.hh b/include/AnalysisManagerMessenger.hh index 8a89ee7..ba305cb 100644 --- a/include/AnalysisManagerMessenger.hh +++ b/include/AnalysisManagerMessenger.hh @@ -52,18 +52,22 @@ class AnalysisManagerMessenger: public G4UImessenger private: - AnalysisManager* anamanager; - - G4UIdirectory* histoDir; - G4UIcmdWithAString* factoryCmd; - G4UIcmdWithABool* saveHitCmd; - G4UIcmdWithABool* saveTrackCmd; - G4UIcmdWithABool* save3DEvdCmd; - G4UIcmdWithABool* save2DEvdCmd; - G4UIcmdWithABool* circleFitCmd; - G4UIcmdWithABool* saveActsCmd; - G4UIcmdWithAnInteger* histoCmd; - G4UIcmdWithAString* addDiffusionCmd; + AnalysisManager* fAnalysisManager; + + G4UIdirectory* fOutDir; + G4UIcmdWithAString* fFileCmd; + G4UIcmdWithABool* fSaveTrackCmd; + + G4UIdirectory* fFLArEDir; + G4UIcmdWithABool* fEnableFLArEOutCmd; + G4UIcmdWithABool* fSave3DEvdCmd; + G4UIcmdWithABool* fSave2DEvdCmd; + G4UIcmdWithABool* fPseudoRecoCmd; + G4UIcmdWithAString* fAddDiffusionCmd; + + G4UIdirectory* fFASER2Dir; + G4UIcmdWithABool* fSaveActsCmd; + }; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/include/LArBoxHit.hh b/include/LArBoxHit.hh index 67424a2..00d70f7 100644 --- a/include/LArBoxHit.hh +++ b/include/LArBoxHit.hh @@ -23,6 +23,12 @@ class LArBoxHit : public G4VHit { void SetPID(G4int pid) { fPID = pid; } void SetTID(G4int tid) { fTID = tid; } void SetStepNo(G4int stepno) { fStepno = stepno; } + //fran added + void SetCopyNum(G4int boxNum) { fBoxNum = boxNum; } + void SetPDG(G4int PDG) { fPDG = PDG; } + void SetDeltaMom(G4ThreeVector& TrackDeltaMom) { fDeltaMom = TrackDeltaMom; } + void SetTime(G4double time) { fTime = time; } + //end fran added void SetPreStepPosition(G4ThreeVector& PreStepPosition) { fPreStepPosition = PreStepPosition; } void SetPostStepPosition(G4ThreeVector& PostStepPosition) { fPostStepPosition = PostStepPosition; } void SetInitMomentum(G4ThreeVector& InitMomentum) { fInitMomentum = InitMomentum; } @@ -41,20 +47,24 @@ class LArBoxHit : public G4VHit { //getter G4int GetTrackStatus() const { return fTrackStatus; } G4ThreeVector GetTrackVertex() const { return fTrackVertex; } - G4double GetTrackLength() const { return fTrackLength; } + G4float GetTrackLength() const { return fTrackLength; } G4int GetParticle() const { return fParticle; } - G4double GetParticleMass() const { return fParticleMass; } + G4float GetParticleMass() const { return fParticleMass; } G4int GetPID() const { return fPID; } G4int GetTID() const { return fTID; } G4int GetStepNo() const { return fStepno; } + G4int GetCopyNum() const { return fBoxNum; } + G4int GetPDG() const { return fPDG; } + G4ThreeVector GetDeltaMom() const { return fDeltaMom; } + G4float GetTime() const { return fTime; } G4ThreeVector GetPreStepPosition() const { return fPreStepPosition; } G4ThreeVector GetPostStepPosition() const { return fPostStepPosition; } G4ThreeVector GetInitMomentum() const { return fInitMomentum; } - G4double GetInitKinEnergy() const { return fInitKinEne; } + G4float GetInitKinEnergy() const { return fInitKinEne; } G4String GetCreatorProcess() const { return fCreatorProcess; } G4String GetProcessName() const { return fProcessName; } - G4double GetStepLength() const { return fStepLength; } - G4double GetEdep() const { return fEdep; } + G4float GetStepLength() const { return fStepLength; } + G4float GetEdep() const { return fEdep; } G4ThreeVector GetEdepPosition() const { return fEdepPosition; } G4String GetVolume() const { return fVolumeName; } G4int GetStepStatus() const { return fStepStatus; } @@ -65,16 +75,20 @@ class LArBoxHit : public G4VHit { private: G4int fTrackStatus; G4ThreeVector fTrackVertex; - G4double fTrackLength; + G4float fTrackLength; G4int fParticle; - G4double fParticleMass; + G4float fParticleMass; G4int fPID; G4int fTID; G4int fStepno; + G4int fBoxNum; + G4int fPDG; + G4ThreeVector fDeltaMom; + G4float fTime; G4ThreeVector fPreStepPosition; G4ThreeVector fPostStepPosition; G4ThreeVector fInitMomentum; - G4double fInitKinEne; + G4float fInitKinEne; G4String fCreatorProcess; G4String fProcessName; G4double fStepLength; diff --git a/include/PixelMap3D.hh b/include/PixelMap3D.hh index 0092117..9c92578 100644 --- a/include/PixelMap3D.hh +++ b/include/PixelMap3D.hh @@ -3,11 +3,12 @@ #include -#include -#include -#include +#include "THnSparse.h" +#include "TFile.h" +#include "TDirectory.h" +#include "TH2.h" -#include +#include "G4ThreeVector.hh" #include "FPFNeutrino.hh" @@ -25,7 +26,7 @@ class PixelMap3D { void FillEntry(const Double_t* pos_xyz, const Double_t* vtx_xyz, const Double_t edep, const Int_t idxPrim); void FillEntryWithToyElectronTransportation(const Double_t* pos_xyz, const Double_t* vtx_xyz, Double_t edep, const Int_t idxPrim); void FillEntryWithToySingleElectronTransportation(const Double_t* pos_xyz, const Double_t* vtx_xyz, Double_t edep, const Int_t idxPrim); - void Write2DPMToFile(TFile* thefile); + void Write2DPMToFile(TFile* thefile, TDirectory *thedir); void Process3DPM(hep_hpc::hdf5::File &h5file, FPFNeutrino neutrino, G4bool save3D); // this should go to a Geometry Service class diff --git a/include/PrimaryParticleInformation.hh b/include/PrimaryParticleInformation.hh index 920a066..767dc05 100644 --- a/include/PrimaryParticleInformation.hh +++ b/include/PrimaryParticleInformation.hh @@ -38,12 +38,13 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// Get the particle PDG code inline G4int GetPDG() const { return fPDG; }; - /// Get the particle mass in MeV - inline G4double GetMass() const { return fMass; }; /// Get the particle charge inline G4double GetCharge() const { return fCharge; }; + /// Get the particle mass in MeV + inline G4double GetMass() const { return fMass; }; + /// Get the particle initial momentum. inline G4ThreeVector GetMomentumMC() const { return fMomentumMC; }; @@ -63,7 +64,7 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// Prints the information about the particle. virtual void Print() const; - + private: /// A particle unique ID. @@ -71,13 +72,10 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// A particle type (PDG code). G4int fPDG; - + /// A particle mass G4double fMass; - /// A particle charge - G4double fCharge; - /// A particle initial momentum (from particle generator) G4ThreeVector fMomentumMC; @@ -101,7 +99,7 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// Scattering type of the neutrino interaction if it is from GENIE G4int fScatteringType; - + /// invariant hadronic mass G4double fW; @@ -113,6 +111,11 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// fsl X4 of the neutrino interaction if it is from GENIE TLorentzVector fFSLX4; + + + /// particle's charge + G4double fCharge; + }; #endif diff --git a/include/geometry/GeometricalParameters.hh b/include/geometry/GeometricalParameters.hh index 1c8eb78..50b0718 100644 --- a/include/geometry/GeometricalParameters.hh +++ b/include/geometry/GeometricalParameters.hh @@ -1,7 +1,7 @@ #ifndef GeometricalParameters_hh #define GeometricalParameters_hh -#include +#include #include "G4String.hh" #include "G4ThreeVector.hh" @@ -219,16 +219,7 @@ class GeometricalParameters { // Sensitive detectors void AddSD2List(int idx, std::string val) { fSDNamelist.insert(std::make_pair(idx, val)); } - std::set > GetSDNamelist() { return fSDNamelist; } - - void SetAddFLArE(G4bool val) { fAddFLArE = val; } - bool GetAddFLArE() { return fAddFLArE; } - void SetAddFORMOSA(G4bool val) { fAddFORMOSA = val; } - bool GetAddFORMOSA() { return fAddFORMOSA; } - void SetAddFASERnu2(G4bool val) { fAddFASERnu2 = val; } - bool GetAddFASERnu2() { return fAddFASERnu2; } - void SetAddFASER2(G4bool val) { fAddFASER2 = val; } - bool GetAddFASER2() { return fAddFASER2; } + std::map GetSDNamelist() { return fSDNamelist; } private: //the singleton @@ -238,11 +229,6 @@ class GeometricalParameters { G4double fHallHeadDistance; ///<- distance between the entrance wall and the first detector G4double fHallOffsetX; // x offset of hall center from the LOS G4double fHallOffsetY; // x offset of hall center from the LOS - - bool fAddFLArE; - bool fAddFORMOSA; - bool fAddFASERnu2; - bool fAddFASER2; // rock envelope G4bool fEnableRockEnvelope; @@ -324,7 +310,6 @@ class GeometricalParameters { G4bool fFASER2FillCaloAndWall; // if true then fill volumes with material, otherwise use air to save CPU - // FASERnu2 Emulsion detector G4double fFASERnu2TotalSizeZ; // Emulsion/Tungsten @@ -353,7 +338,7 @@ class GeometricalParameters { G4ThreeVector fFORMOSAPos; // Sensitive detectors - std::set > fSDNamelist; + std::map fSDNamelist; }; #endif diff --git a/include/reco/Barcode.hh b/include/reco/Barcode.hh index 0402c37..cb35729 100644 --- a/include/reco/Barcode.hh +++ b/include/reco/Barcode.hh @@ -8,6 +8,7 @@ #pragma once +//#include "Acts/Utilities/MultiIndex.hpp" #include "reco/MultiIndex.hh" #include @@ -178,4 +179,4 @@ struct hash { return std::hash()(barcode.value()); } }; -} // namespace std \ No newline at end of file +} // namespace std diff --git a/include/reco/CircleFit.hh b/include/reco/CircleFit.hh deleted file mode 100644 index 5cc8e7e..0000000 --- a/include/reco/CircleFit.hh +++ /dev/null @@ -1,120 +0,0 @@ -//////////////////////////////////////////////////////////////////////// -// \file CircularFit.hh -// \brief Class headers to perform circular fits using analytical methods -// for momentum estimation in magnetic fields -// \author M. Vicenzi (mvicenzi@bnl.gov) -//////////////////////////////////////////////////////////////////////// - -#ifndef CIRCULARFIT_HH -#define CIRCULARFIT_HH - -#include - -namespace circularfitter { - struct line { - double m; - double q; - }; - - class CircleFit { - public: - CircleFit(const std::vector x, const std::vector z); - - double GetXc() { return fXc; }; - double GetZc() { return fZc; }; - double GetR() { return fR; }; - double GetChi2() { return fChi2; }; - int GetStatus() { return fStatus; }; - - ~CircleFit(); - - private: - double fXc; - double fZc; - double fR; - double fChi2; - int fStatus; - }; - - class LineFit { - public: - LineFit(const std::vector z, const std::vector y); - - double GetP0() { return fP0; }; - double GetP1() { return fP1; }; - double GetCosDip() { return fCosDip; }; - double GetChi2() { return fChi2; }; - int GetStatus() { return fStatus; }; - line GetLine(); - - ~LineFit(); - - private: - double fP0; - double fP1; - double fCosDip; - double fChi2; - int fStatus; - }; - - struct hit { - double x; - double y; - double z; - }; - - class CircleExtractor { - public: - CircleExtractor(const std::vector x, const std::vector y, const std::vector z); - - // splits hits between pre and post magnet(s) according to list of magnet z position(s). - void splitHits(std::vector magnetZs, - std::vector x, std::vector y, std::vector z, - std::vector> &hits); - - // for a given magnet, compute circle from pre and post hit collections + magnet z limits - void computeCircle(std::vector zpre, std::vector xpre, std::vector zpost, - std::vector xpost, double Zin, double Zout); - - // support functions - double extrapolateLine(double z, line l); - std::pair getPerpLineIntersection(double z1, double x1, double z2, double x2, line l1, line l2); - double getR(double za, double xa, double zb, double xb); - - // return methods - std::vector GetMagnetZs() { return fzpos; }; - std::vector GetXc() { return fXc; }; - std::vector GetZc() { return fZc; }; - std::vector GetR(); - std::vector GetPreLine() { return fpre; }; - std::vector GetPostLine() { return fpost; }; - - ~CircleExtractor(); - - private: - std::vector fpre; - std::vector fpost; - std::vector fzpos; - std::vector fZc; - std::vector fXc; - std::vector fR1; - std::vector fR2; - }; - - class ParabolicFit { - public: - ParabolicFit(const std::vector z, const std::vector x, double r0); - - double GetA() { return fA; }; - double GetB() { return fB; }; - double GetR() { return fR; }; - - ~ParabolicFit(); - - private: - double fA; - double fB; - double fR; - }; -} -#endif diff --git a/include/reco/MultiIndex.hh b/include/reco/MultiIndex.hh index 22124e6..23f6a75 100644 --- a/include/reco/MultiIndex.hh +++ b/include/reco/MultiIndex.hh @@ -163,4 +163,4 @@ struct hash> { return std::hash()(idx.value()); } }; -} // namespace std \ No newline at end of file +} // namespace std diff --git a/macros/backgrounds.mac b/macros/backgrounds.mac index b1704d2..a65515a 100644 --- a/macros/backgrounds.mac +++ b/macros/backgrounds.mac @@ -10,15 +10,15 @@ /gen/bkg/backgroundWindow 187.5 us # define output options -/histo/save3DEvd false -/histo/save2DEvd true -/histo/saveHit false -/histo/addDiffusion false -/histo/fileName test_backgrounds.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/addDiffusion false +/out/fileName test_backgrounds.root # store and save trajectories for EVD /tracking/storeTrajectory 1 -/histo/saveTrack true +/out/saveTrack true # shoot background equivalent to 1 "spill" /run/beamOn 1 diff --git a/macros/basic.mac b/macros/basic.mac index e166aa2..ef5d05b 100644 --- a/macros/basic.mac +++ b/macros/basic.mac @@ -19,10 +19,10 @@ /gps/ene/mono 1000 GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName 1000GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName 1000GeV.root # shoot 5000 particles /run/beamOn 5000 diff --git a/macros/dark_photon_hepmc.mac b/macros/dark_photon_hepmc.mac index a2690c0..999c2f7 100644 --- a/macros/dark_photon_hepmc.mac +++ b/macros/dark_photon_hepmc.mac @@ -10,12 +10,12 @@ /gen/hepmc/placeInDecayVolume true # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit true -/histo/addDiffusion false -/histo/circleFit true -/histo/fileName test_dark_photon_hepmc.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/flare/addDiffusion false + +/out/fileName test_dark_photon_hepmc.root # shoot 5 particles /run/beamOn 5 diff --git a/macros/dark_photon_hepmc2.mac b/macros/dark_photon_hepmc2.mac index 42b5ce3..c951ddc 100644 --- a/macros/dark_photon_hepmc2.mac +++ b/macros/dark_photon_hepmc2.mac @@ -11,12 +11,12 @@ /gen/placeInDecayVolume true # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit true -/histo/addDiffusion false -/histo/circleFit true -/histo/fileName test_dark_photon_hepmc.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/flare/addDiffusion false + +/out/fileName test_dark_photon_hepmc.root # shoot 5 particles /run/beamOn 5 diff --git a/macros/legacy/LAr_e-.mac b/macros/legacy/LAr_e-.mac index 9b536fa..d33af7a 100644 --- a/macros/legacy/LAr_e-.mac +++ b/macros/legacy/LAr_e-.mac @@ -24,6 +24,6 @@ ## 200MeV, 500MeV, 1GeV, 5GeV, 10GeV, 50GeV, 100GeV ## 200GeV, 300GeV, 500GeV, 800GeV, 1000GeV /gps/ene/mono 5 GeV -/histo/fileName 5GeV.root +/out/fileName 5GeV.root /run/beamOn 1 #/run/beamOn 10000 diff --git a/macros/legacy/LAr_mu-.mac b/macros/legacy/LAr_mu-.mac index 228e936..210ca64 100644 --- a/macros/legacy/LAr_mu-.mac +++ b/macros/legacy/LAr_mu-.mac @@ -13,6 +13,6 @@ /gps/ene/mono 5 GeV -/histo/fileName 5GeV.root +/out/fileName 5GeV.root /run/beamOn 1 diff --git a/macros/legacy/LAr_tau-.mac b/macros/legacy/LAr_tau-.mac index 887169c..28103e0 100644 --- a/macros/legacy/LAr_tau-.mac +++ b/macros/legacy/LAr_tau-.mac @@ -22,5 +22,5 @@ ## 200MeV, 500MeV, 1GeV, 5GeV, 10GeV, 50GeV, 100GeV ## 200GeV, 300GeV, 500GeV, 800GeV, 1000GeV /gps/ene/mono 5 GeV -/histo/fileName 5GeV.root +/out/fileName 5GeV.root /run/beamOn 10000 diff --git a/macros/legacy/LKr_e-.mac b/macros/legacy/LKr_e-.mac index a1bede5..486462a 100644 --- a/macros/legacy/LKr_e-.mac +++ b/macros/legacy/LKr_e-.mac @@ -21,5 +21,5 @@ ## 200MeV, 500MeV, 1GeV, 5GeV, 10GeV, 50GeV, 100GeV ## 200GeV, 300GeV, 500GeV, 800GeV, 1000GeV /gps/ene/mono 5 GeV -/histo/fileName 5GeV.root +/out/fileName 5GeV.root /run/beamOn 10000 diff --git a/macros/legacy/LKr_mu-.mac b/macros/legacy/LKr_mu-.mac index eb0983b..26d70dd 100644 --- a/macros/legacy/LKr_mu-.mac +++ b/macros/legacy/LKr_mu-.mac @@ -21,5 +21,5 @@ ## 200MeV, 500MeV, 1GeV, 5GeV, 10GeV, 50GeV, 100GeV ## 200GeV, 300GeV, 500GeV, 800GeV, 1000GeV /gps/ene/mono 5 GeV -/histo/fileName 5GeV.root +/out/fileName 5GeV.root /run/beamOn 10000 diff --git a/macros/legacy/LKr_tau-.mac b/macros/legacy/LKr_tau-.mac index e84271b..f826776 100644 --- a/macros/legacy/LKr_tau-.mac +++ b/macros/legacy/LKr_tau-.mac @@ -21,5 +21,5 @@ ## 200MeV, 500MeV, 1GeV, 5GeV, 10GeV, 50GeV, 100GeV ## 200GeV, 300GeV, 500GeV, 800GeV, 1000GeV /gps/ene/mono 5 GeV -/histo/fileName 5GeV.root +/out/fileName 5GeV.root /run/beamOn 10000 diff --git a/macros/nue_genie.mac b/macros/nue_genie.mac index 7f4b98f..ee60d11 100644 --- a/macros/nue_genie.mac +++ b/macros/nue_genie.mac @@ -10,11 +10,11 @@ /gen/genie/genieIStart 5 # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/addDiffusion false -/histo/fileName test_genie_nue.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/flare/addDiffusion false +/out/fileName test_genie_nue.root # shoot 5 particles /run/beamOn 5 diff --git a/macros/numu_genie.mac b/macros/numu_genie.mac index 1474dd8..ebe04dd 100644 --- a/macros/numu_genie.mac +++ b/macros/numu_genie.mac @@ -10,11 +10,12 @@ /gen/genie/genieIStart 5 # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/addDiffusion false -/histo/fileName test_genie_numu.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false +/out/flare/pseudoReco false + +/out/flare/addDiffusion false +/out/fileName test_genie_numu.root # shoot 5 particles /run/beamOn 5 diff --git a/macros/nutau_genie.mac b/macros/nutau_genie.mac index 3a48d45..458ccb3 100644 --- a/macros/nutau_genie.mac +++ b/macros/nutau_genie.mac @@ -10,11 +10,11 @@ /gen/genie/genieIStart 5 # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/addDiffusion false -/histo/fileName test_genie_nutau.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/flare/addDiffusion false +/out/fileName test_genie_nutau.root # shoot 5 particles /run/beamOn 5 diff --git a/macros/single_particle/LAr_Kminus_mono.mac b/macros/single_particle/LAr_Kminus_mono.mac index 81955e5..c434d78 100644 --- a/macros/single_particle/LAr_Kminus_mono.mac +++ b/macros/single_particle/LAr_Kminus_mono.mac @@ -14,10 +14,10 @@ /gps/ene/mono 49.50876 GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName Kminus_momentum_50GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName Kminus_momentum_50GeV.root # shoot 1 particles /run/beamOn 1 diff --git a/macros/single_particle/LAr_e-_mono.mac b/macros/single_particle/LAr_e-_mono.mac index f137e44..28ea895 100644 --- a/macros/single_particle/LAr_e-_mono.mac +++ b/macros/single_particle/LAr_e-_mono.mac @@ -14,10 +14,10 @@ /gps/ene/mono 19.999489 GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName e-_momentum_20GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName e-_momentum_20GeV.root # shoot 1 particles /run/beamOn 1 diff --git a/macros/single_particle/LAr_gamma_mono.mac b/macros/single_particle/LAr_gamma_mono.mac index ce46aca..1b0f38f 100644 --- a/macros/single_particle/LAr_gamma_mono.mac +++ b/macros/single_particle/LAr_gamma_mono.mac @@ -14,10 +14,10 @@ /gps/ene/mono 20. GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName gamma_momentum_20GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName gamma_momentum_20GeV.root # shoot 1 particles /run/beamOn 1 diff --git a/macros/single_particle/LAr_mu-_mono.mac b/macros/single_particle/LAr_mu-_mono.mac index 7b03b05..d60eb57 100644 --- a/macros/single_particle/LAr_mu-_mono.mac +++ b/macros/single_particle/LAr_mu-_mono.mac @@ -1,5 +1,6 @@ # define detector material before run initialization -/control/execute macros/geometry_options/FPF_hall_Reference.mac +/control/execute macros/geometry_options/FLArE_only_newBabyMIND.mac +/det/saveGdml false /run/initialize /random/setSeeds 324199 420475 @@ -8,16 +9,12 @@ /gen/select gun /gps/verbose 0 /gps/pos/type Point -/gps/pos/centre 0 0 -2.5 m +/gps/pos/centre 0 0 7.8 m /gps/direction 0 0 1 /gps/particle mu- -/gps/ene/mono 19.894621 GeV +/gps/ene/mono 5.0 GeV -# define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName mu-_momentum_20GeV.root +/out/fileName mu-_5GeV_newBabyMIND.root # shoot 1 particles -/run/beamOn 1 +/run/beamOn 50 diff --git a/macros/single_particle/LAr_neutron_mono.mac b/macros/single_particle/LAr_neutron_mono.mac index 2307594..462bff1 100644 --- a/macros/single_particle/LAr_neutron_mono.mac +++ b/macros/single_particle/LAr_neutron_mono.mac @@ -14,10 +14,10 @@ /gps/ene/mono 4.147947 GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName neutron_momentum_5GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName neutron_momentum_5GeV.root # shoot 1 particles /run/beamOn 1 diff --git a/macros/single_particle/LAr_pi0_mono.mac b/macros/single_particle/LAr_pi0_mono.mac index 1a94a97..b9b0baa 100644 --- a/macros/single_particle/LAr_pi0_mono.mac +++ b/macros/single_particle/LAr_pi0_mono.mac @@ -14,10 +14,10 @@ /gps/ene/mono 24.865388 GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName pi0_momentum_25GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName pi0_momentum_25GeV.root # shoot 1 particles /run/beamOn 1 diff --git a/macros/single_particle/LAr_piminus_mono.mac b/macros/single_particle/LAr_piminus_mono.mac index 15c0e67..483b4dc 100644 --- a/macros/single_particle/LAr_piminus_mono.mac +++ b/macros/single_particle/LAr_piminus_mono.mac @@ -14,10 +14,10 @@ /gps/ene/mono 24.860819 GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName piminus_momentum_25GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName piminus_momentum_25GeV.root # shoot 1 particles /run/beamOn 1 diff --git a/macros/single_particle/LAr_proton_mono.mac b/macros/single_particle/LAr_proton_mono.mac index 23e7f3a..99a63b5 100644 --- a/macros/single_particle/LAr_proton_mono.mac +++ b/macros/single_particle/LAr_proton_mono.mac @@ -14,10 +14,10 @@ /gps/ene/mono 7.116562 GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName proton_momentum_8GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName proton_momentum_8GeV.root # shoot 1 particles /run/beamOn 1 diff --git a/macros/single_particle/LAr_tau-_mono.mac b/macros/single_particle/LAr_tau-_mono.mac index 3aa366e..32e7a61 100644 --- a/macros/single_particle/LAr_tau-_mono.mac +++ b/macros/single_particle/LAr_tau-_mono.mac @@ -14,10 +14,10 @@ /gps/ene/mono 18.301916 GeV # define output options -/histo/save3DEvd false -/histo/save2DEvd false -/histo/saveHit false -/histo/fileName tau-_momentum_20GeV.root +/out/flare/save3DEvd false +/out/flare/save2DEvd false + +/out/fileName tau-_momentum_20GeV.root # shoot 1 particles /run/beamOn 1 diff --git a/sdcc_setup.sh b/sdcc_setup.sh new file mode 100644 index 0000000..8e744af --- /dev/null +++ b/sdcc_setup.sh @@ -0,0 +1,13 @@ +## get software stack from LCG in CVMFS +## this setups both ROOT and Geant4 +## docs: https://lcginfo.cern.ch/ +echo "Setting up LGC_105 software stack..." +source /cvmfs/sft.cern.ch/lcg/views/LCG_105/x86_64-centos7-gcc11-opt/setup.sh + +## HEP_HPC: this is not available in LGC: it's a Fermilab/NOvA-specific package +## DUNE cvmfs: compiled with an older HDF5 than LGC... +## Current solution: build HEP_HPC locally on SDCC with HDF5 from LGC. +## Use env variable to point to it later in CMakeLists.txt +export hep_hpc_path="/lbne/u/mvicenzi/hep_hpc/bin" + +echo "Setup completed!" diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index a80d687..7c1e96b 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -2,25 +2,10 @@ #include #include #include -#include +#include #include #include -#include "AnalysisManager.hh" -#include "geometry/GeometricalParameters.hh" -#include "LArBoxSD.hh" -#include "FASER2TrackerSD.hh" -#include "LArBoxHit.hh" -#include "PrimaryParticleInformation.hh" -#include "reco/PCAAnalysis3D.hh" -#include "reco/Cluster3D.hh" -#include "reco/LinearFit.hh" -#include "reco/ShowerLID.hh" -#include "reco/CircleFit.hh" -#include "reco/Barcode.hh" -#include "FPFParticle.hh" -#include "FPFNeutrino.hh" - #include #include #include @@ -28,6 +13,7 @@ #include #include +#include #include #include #include @@ -35,370 +21,382 @@ #include #include -using namespace hep_hpc::hdf5; +#include "AnalysisManager.hh" +#include "geometry/GeometricalParameters.hh" +#include "LArBoxSD.hh" +#include "FASER2TrackerSD.hh" +#include "LArBoxHit.hh" +#include "PrimaryParticleInformation.hh" +#include "reco/PCAAnalysis3D.hh" +#include "reco/Cluster3D.hh" +#include "reco/LinearFit.hh" +#include "reco/ShowerLID.hh" +#include "reco/Barcode.hh" +#include "FPFParticle.hh" +#include "FPFNeutrino.hh" -AnalysisManager* AnalysisManager::instance = 0; +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- +// AnalysisManager "singleton" instance +// once initialized, can be used to point to AnalysisManager +// from anywhere else in the codebase +AnalysisManager *AnalysisManager::fInstance = 0; -AnalysisManager* AnalysisManager::GetInstance() { - if (!instance) { - std::cout<<"AnalysisManager: Re-initialization"<Branch("evtID" , &evtID , "evtID/I"); - evt->Branch("FPFNeutrino" , &neutrino, 96000, 0); - evt->Branch("FPFParticle" , &primaries, 96000, 0); - evt->Branch("TotalDedxLongitudinal" , TotalDedxLongitudinal , "TotalDedxLongitudinal[3000]/D"); - evt->Branch("TrueTotalDedxLongitudinal" , TrueTotalDedxLongitudinal , "TrueTotalDedxLongitudinal[3000]/D"); - evt->Branch("nPrimaryParticle" , &nPrimaryParticle , "nPrimaryParticle/I"); - evt->Branch("primaryParentPDG" , primaryParentPDG , "primaryParentPDG[nPrimaryParticle]/I"); - evt->Branch("primaryTrackLength" , primaryTrackLength , "primaryTrackLength[nPrimaryParticle]/D"); - evt->Branch("primaryTrackLengthInTPC" , primaryTrackLengthInTPC , "primaryTrackLengthInTPC[nPrimaryParticle]/D"); - evt->Branch("ProngEInDetector" , ProngEInDetector , "ProngEInDetector[nPrimaryParticle]/D"); - evt->Branch("ProngEInLAr" , ProngEInLAr , "ProngEInLAr[nPrimaryParticle]/D"); - evt->Branch("ProngEInHadCal" , ProngEInHadCal , "ProngEInHadCal[nPrimaryParticle]/D"); - evt->Branch("ProngEInMuonFinder" , ProngEInMuonFinder , "ProngEInMuonFinder[nPrimaryParticle]/D"); - evt->Branch("ProngEInMuonFinderLayer1X" , ProngEInMuonFinderLayer1X , "ProngEInMuonFinderLayer1X[nPrimaryParticle]/D"); - evt->Branch("ProngEInMuonFinderLayer1Y" , ProngEInMuonFinderLayer1Y , "ProngEInMuonFinderLayer1Y[nPrimaryParticle]/D"); - evt->Branch("ProngEInMuonFinderLayer2X" , ProngEInMuonFinderLayer2X , "ProngEInMuonFinderLayer2X[nPrimaryParticle]/D"); - evt->Branch("ProngEInMuonFinderLayer2Y" , ProngEInMuonFinderLayer2Y , "ProngEInMuonFinderLayer2Y[nPrimaryParticle]/D"); - evt->Branch("ProngAngleToBeamDir" , ProngAngleToBeamDir , "ProngAngleToBeamDir[nPrimaryParticle]/D"); - evt->Branch("ShowerLength" , ShowerLength , "ShowerLength[nPrimaryParticle]/D"); - evt->Branch("ShowerLengthInLAr" , ShowerLengthInLAr , "ShowerLengthInLAr[nPrimaryParticle]/D"); - evt->Branch("ShowerWidth" , ShowerWidth , "ShowerWidth[nPrimaryParticle]/D"); - evt->Branch("ShowerWidthInLAr" , ShowerWidthInLAr , "ShowerWidthInLAr[nPrimaryParticle]/D"); - evt->Branch("ProngAvgdEdx" , ProngAvgdEdx , "ProngAvgdEdx[nPrimaryParticle]/D"); - evt->Branch("ProngAvgdEdxInLAr" , ProngAvgdEdxInLAr , "ProngAvgdEdxInLAr[nPrimaryParticle]/D"); - evt->Branch("ProngdEdxAlongTrack" , ProngdEdxAlongTrack , "ProngdEdxAlongTrack[nPrimaryParticle][100]/D"); - evt->Branch("ProngdEdxTrackLength" , ProngdEdxTrackLength , "ProngdEdxTrackLength[nPrimaryParticle][100]/I"); - evt->Branch("dir_pol_x" , dir_pol_x , "dir_pol_x[nPrimaryParticle]/D"); - evt->Branch("dir_pol_y" , dir_pol_y , "dir_pol_y[nPrimaryParticle]/D"); - evt->Branch("dir_pol_z" , dir_pol_z , "dir_pol_z[nPrimaryParticle]/D"); - evt->Branch("dir_coc_x" , dir_coc_x , "dir_coc_x[nPrimaryParticle]/D"); - evt->Branch("dir_coc_y" , dir_coc_y , "dir_coc_y[nPrimaryParticle]/D"); - evt->Branch("dir_coc_z" , dir_coc_z , "dir_coc_z[nPrimaryParticle]/D"); - - evt->Branch("nHits" , &nHits , "nHits/I"); - evt->Branch("sparseFractionMem" , &sparseFractionMem , "sparseFractionMem/D"); - evt->Branch("sparseFractionBins" , &sparseFractionBins , "sparseFractionBins/D"); - if (m_saveHit) { - evt->Branch("HitTID" , HitTID , "HitTID[nHits]/I"); - evt->Branch("HitPID" , HitPID , "HitPID[nHits]/I"); - evt->Branch("HitPDG" , HitPDG , "HitPDG[nHits]/I"); - evt->Branch("HitTrackStatus" , HitTrackStatus , "HitTrackStatus[nHits]/I"); - evt->Branch("HitPrePositionX" , HitPrePositionX , "HitPrePositionX[nHits]/D"); - evt->Branch("HitPrePositionY" , HitPrePositionY , "HitPrePositionY[nHits]/D"); - evt->Branch("HitPrePositionZ" , HitPrePositionZ , "HitPrePositionZ[nHits]/D"); - evt->Branch("HitPosPositionX" , HitPosPositionX , "HitPosPositionX[nHits]/D"); - evt->Branch("HitPosPositionY" , HitPosPositionY , "HitPosPositionY[nHits]/D"); - evt->Branch("HitPosPositionZ" , HitPosPositionZ , "HitPosPositionZ[nHits]/D"); - evt->Branch("HitEdep" , HitEdep , "HitEdep[nHits]/D"); - } + fMessenger = new AnalysisManagerMessenger(this); - evt->Branch("edepInLAr" , &edepInLAr , "edepInLAr/D"); - evt->Branch("edepInHadCalX" , &edepInHadCalX , "edepInHadCalX/D"); - evt->Branch("edepInHadCalY" , &edepInHadCalY , "edepInHadCalY/D"); - evt->Branch("edepInMuonFinderX" , &edepInMuonFinderX , "edepInMuonFinderX/D"); - evt->Branch("edepInMuonFinderY" , &edepInMuonFinderY , "edepInMuonFinderY/D"); - evt->Branch("edepInHadAborb" , &edepInHadAborb , "edepInHadAborb/D"); - evt->Branch("edepInMuonFinderAbsorb" , &edepInMuonFinderAbsorb , "edepInMuonFinderAbsorb/D"); - evt->Branch("missCountedEnergy" , &missCountedEnergy , "missCountedEnergy/D"); + fEvt = nullptr; + fTrk = nullptr; + fPrim = nullptr; + fFLArEHits = nullptr; + fFLArEHCALHits = nullptr; + fFLArEPseudoReco = nullptr; + fActsHitsTree = nullptr; + fActsParticlesTree = nullptr; - evt->Branch("nFromFSLParticles" , &nFromFSLParticles , "nFromFSLParticles/I"); - - if(m_circularFit){ - evt->Branch("circNhits" , &circNhits , "circNhits/I"); - evt->Branch("trkNhits" , &trkNhits , "trkNhits/I"); - evt->Branch("circXc" , &xc , "circXc/D"); - evt->Branch("circZc" , &zc , "circZc/D"); - evt->Branch("circRc" , &rc , "circRc/D"); - evt->Branch("circp0" , &p0 , "circp0/D"); - evt->Branch("circp1" , &p1 , "circp1/D"); - evt->Branch("circcosDip" , &cosDip , "circcosDip/D"); - evt->Branch("Nmagnets" , &Nmagnets , "Nmagnets/I"); - evt->Branch("trkp0" , &trkp0 , "trkp0/D"); - evt->Branch("trkp1" , &trkp1 , "trkp1/D"); - evt->Branch("trkcosDip" , &trkcosDip , "trkcosDip/D"); - evt->Branch("magnetZs" , &magzpos); - evt->Branch("trkXc" , &trkxc); - evt->Branch("trkZc" , &trkzc); - evt->Branch("trkRc" , &trkrc); - evt->Branch("trkqIn" , &trkqIn); - evt->Branch("trkmIn" , &trkmIn); - evt->Branch("trkqOut" , &trkqOut); - evt->Branch("trkmOut" , &trkmOut); - evt->Branch("circHitXFSL" , &hitXFSL); - evt->Branch("circHitYFSL" , &hitYFSL); - evt->Branch("circHitZFSL" , &hitZFSL); - evt->Branch("circHitPFSL" , &hitPFSL); - evt->Branch("trkHitXFSL" , &trkXFSL); - evt->Branch("trkHitYFSL" , &trkYFSL); - evt->Branch("trkHitZFSL" , &trkZFSL); - evt->Branch("trkHitPFSL" , &trkPFSL); - } + fSaveTrack = false; + fSave3DEvd = false; + fSave2DEvd = false; + fSavePseudoReco = false; + fEnableFLArE = true; + fSaveActs = true; + fAddDiffusion = "false"; +} - // Defaults to true if FASER2 is enabled. But give the option to disable output if required - m_saveActs = m_saveActs && GeometricalParameters::Get()->GetAddFASER2(); - if (m_saveActs) { - //* Acts Hits Tree [i == unsigned int; F == float; l == Long unsigned 64 int] - acts_hits_tree = new TTree("hits", "ActsHitsTree"); - acts_hits_tree->Branch("event_id" , &ActsHitsEventID, "event_id/i"); - acts_hits_tree->Branch("geometry_id" , &ActsHitsGeometryID, "geometryid/l"); - acts_hits_tree->Branch("particle_id" , &ActsHitsParticleID, "particle_id/l"); - acts_hits_tree->Branch("tx" , &ActsHitsX, "tx/F"); - acts_hits_tree->Branch("ty" , &ActsHitsY, "ty/F"); - acts_hits_tree->Branch("tz" , &ActsHitsZ, "tz/F"); - acts_hits_tree->Branch("tt" , &ActsHitsT, "tt/F"); - acts_hits_tree->Branch("tpx" , &ActsHitsPx, "tpx/F"); - acts_hits_tree->Branch("tpy" , &ActsHitsPy, "tpy/F"); - acts_hits_tree->Branch("tpz" , &ActsHitsPz, "tpz/F"); - acts_hits_tree->Branch("te" , &ActsHitsE, "tpe/F"); - acts_hits_tree->Branch("deltapx" , &ActsHitsDeltaPx, "deltapx/F"); - acts_hits_tree->Branch("deltapy" , &ActsHitsDeltaPy, "deltapy/F"); - acts_hits_tree->Branch("deltapz" , &ActsHitsDeltaPz, "deltapz/F"); - acts_hits_tree->Branch("deltae" , &ActsHitsDeltaE, "deltae/F"); - acts_hits_tree->Branch("index" , &ActsHitsIndex, "index/I"); - acts_hits_tree->Branch("volume_id" , &ActsHitsVolumeID, "volume_id/i"); - acts_hits_tree->Branch("boundary_id" , &ActsHitsBoundaryID, "boundary_id/i"); - acts_hits_tree->Branch("layer_id" , &ActsHitsLayerID, "layer_id/i"); - acts_hits_tree->Branch("approach_id" , &ActsHitsApproachID, "approach_id/i"); - acts_hits_tree->Branch("sensitive_id" , &ActsHitsSensitiveID, "sensitive_id/i"); - - //* Acts truth particle tree - acts_particles_tree = new TTree("particles", "ActsParticlesTree"); - acts_particles_tree->Branch("event_id" , &ActsHitsEventID, "event_id/i"); - acts_particles_tree->Branch("particle_id" , &ActsParticlesParticleId); - acts_particles_tree->Branch("particle_type" , &ActsParticlesParticleType); - acts_particles_tree->Branch("process" , &ActsParticlesProcess); - acts_particles_tree->Branch("vx" , &ActsParticlesVx); - acts_particles_tree->Branch("vy" , &ActsParticlesVy); - acts_particles_tree->Branch("vz" , &ActsParticlesVz); - acts_particles_tree->Branch("vt" , &ActsParticlesVt); - acts_particles_tree->Branch("px" , &ActsParticlesPx); - acts_particles_tree->Branch("py" , &ActsParticlesPy); - acts_particles_tree->Branch("pz" , &ActsParticlesPz); - acts_particles_tree->Branch("m" , &ActsParticlesM); - acts_particles_tree->Branch("q" , &ActsParticlesQ); - acts_particles_tree->Branch("eta" , &ActsParticlesEta); - acts_particles_tree->Branch("phi" , &ActsParticlesPhi); - acts_particles_tree->Branch("pt" , &ActsParticlesPt); - acts_particles_tree->Branch("p" , &ActsParticlesP); - acts_particles_tree->Branch("vertex_primary" , &ActsParticlesVertexPrimary); - acts_particles_tree->Branch("vertex_secondary" , &ActsParticlesVertexSecondary); - acts_particles_tree->Branch("particle" , &ActsParticlesParticle); - acts_particles_tree->Branch("generation" , &ActsParticlesGeneration); - acts_particles_tree->Branch("sub_particle" , &ActsParticlesSubParticle); - acts_particles_tree->Branch("e_loss" , &ActsParticlesELoss); - acts_particles_tree->Branch("total_x0" , &ActsParticlesPathInX0); - acts_particles_tree->Branch("total_l0" , &ActsParticlesPathInL0); - acts_particles_tree->Branch("number_of_hits" , &ActsParticlesNumberOfHits); - acts_particles_tree->Branch("outcome" , &ActsParticlesOutcome); - } +AnalysisManager::~AnalysisManager() {} + +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::bookEvtTree() +{ + fEvt = new TTree("event", "event info"); + fEvt->Branch("evtID", &evtID, "evtID/I"); + + //TODO EXPAND!!!! +} +void AnalysisManager::bookPrimTree() +{ + fPrim = new TTree("primaries", "primaries info"); + fPrim->Branch("evtID", &evtID, "evtID/I"); + fPrim->Branch("vtxID", &primVtxID, "vtxID/I"); + fPrim->Branch("PDG", &primPDG, "PDG/I"); + fPrim->Branch("trackID", &primTrackID, "trackID/I"); + fPrim->Branch("barcode", &primParticleID, "bardcode/I"); + fPrim->Branch("mass", &primM, "mass/F"); + fPrim->Branch("charge", &primQ, "charge/F"); + fPrim->Branch("Vx", &primVx, "Vx/F"); // position + fPrim->Branch("Vy", &primVy, "Vy/F"); + fPrim->Branch("Vz", &primVz, "Vz/F"); + fPrim->Branch("Vt", &primVt, "Vt/F"); + fPrim->Branch("Px", &primPx, "Px/F"); // momentum + fPrim->Branch("Py", &primPy, "Py/F"); + fPrim->Branch("Pz", &primPz, "Pz/F"); + fPrim->Branch("E", &primE, "E/F"); // initial total energy + fPrim->Branch("KE", &primKE, "KE/F"); // initial kinetic energy + fPrim->Branch("Eta", &primEta, "Eta/F"); + fPrim->Branch("Phi", &primPhi, "Phi/F"); + fPrim->Branch("Pt", &primPt, "Pt/F"); + fPrim->Branch("P", &primP, "P/F"); } -void AnalysisManager::bookTrkTree() { +void AnalysisManager::bookTrkTree() +{ + fTrk = new TTree("trajectories", "trajectories info"); + fTrk->Branch("evtID", &evtID, "evtID/I"); + fTrk->Branch("trackTID", &trackTID, "trackTID/I"); + fTrk->Branch("trackPID", &trackPID, "trackPID/I"); + fTrk->Branch("trackPDG", &trackPDG, "trackPDG/I"); + fTrk->Branch("trackKinE", &trackKinE, "trackKinE/D"); + fTrk->Branch("trackNPoints", &trackNPoints, "trackNPoints/I"); + fTrk->Branch("trackPointX", &trackPointX); + fTrk->Branch("trackPointY", &trackPointY); + fTrk->Branch("trackPointZ", &trackPointZ); +} - trk = new TTree("trk", "trkTreeInfo"); - trk->Branch("evtID", &evtID, "evtID/I"); - trk->Branch("trackTID", &trackTID, "trackTID/I"); - trk->Branch("trackPID", &trackPID, "trackPID/I"); - trk->Branch("trackPDG", &trackPDG, "trackPDG/I"); - trk->Branch("trackKinE", &trackKinE, "trackKinE/D"); - trk->Branch("trackNPoints", &trackNPoints, "trackNPoints/I"); - trk->Branch("trackPointX", &trackPointX); - trk->Branch("trackPointY", &trackPointY); - trk->Branch("trackPointZ", &trackPointZ); +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- +void AnalysisManager::bookFLArETrees() +{ + // create subdirectory in file + fFLArEDir = fFile->mkdir("flare","FLArE output",kTRUE); + fFile->cd(fFLArEDir->GetName()); + + // FLArE hits tree + fFLArEHits = new TTree("flare_hits", "FLArE hits info"); + fFLArEHits->Branch("flareEvtID", &evtID, "flareEvtID/I"); + fFLArEHits->Branch("flareTrackID", &flareTrackID, "flareTrackID/I"); + fFLArEHits->Branch("flareBarcode", &flareParticleID, "flareParticleID/I"); + fFLArEHits->Branch("flareParentID", &flareParentID, "flareParentID/I"); + fFLArEHits->Branch("flarePDG", &flarePDG, "flarePDG/I"); + fFLArEHits->Branch("flareCopyNum", &flareCopyNum, "flareCopyNum/I"); + fFLArEHits->Branch("flareT", &flareT, "flareT/I"); + fFLArEHits->Branch("flareX", &flareX, "flareX/F"); // Pre-position + fFLArEHits->Branch("flareY", &flareY, "flareY/F"); + fFLArEHits->Branch("flareZ", &flareZ, "flareZ/F"); + fFLArEHits->Branch("flarePx", &flarePx, "flarePx/F"); // momentum + fFLArEHits->Branch("flarePy", &flarePy, "flarePy/F"); + fFLArEHits->Branch("flarePz", &flarePz, "flarePz/F"); + fFLArEHits->Branch("flareDeltaPx", &flareDeltaPx, "flareDeltaPx/F"); + fFLArEHits->Branch("flareDeltaPy", &flareDeltaPy, "flareDeltaPy/F"); + fFLArEHits->Branch("flareDeltaPz", &flareDeltaPz, "flareDeltaPz/F"); + fFLArEHits->Branch("flareEdep", &flareEdep, "flareEdep/F"); + + // FLArE HCAL hits + fFLArEHCALHits = new TTree("hcal_hits", "FLArE HCAL hits info"); + + fFLArEHCALHits->Branch("hadEvtID", &evtID, "hadEvtID/I"); + fFLArEHCALHits->Branch("hadTrackID", &flareTrackID, "hadTrackID/I"); + fFLArEHCALHits->Branch("hadBarcode", &flareParticleID, "hadParticleID/I"); + fFLArEHCALHits->Branch("hadParentID", &flareParentID, "hadParentID/I"); + fFLArEHCALHits->Branch("hadPDG", &flarePDG, "hadPDG/I"); + fFLArEHCALHits->Branch("hadCopyNum", &flareCopyNum, "hadCopyNum/I"); + fFLArEHCALHits->Branch("hadT", &flareT, "hadT/I"); + fFLArEHCALHits->Branch("hadX", &flareX, "hadX/F"); // Pre-position + fFLArEHCALHits->Branch("hadY", &flareY, "hadY/F"); + fFLArEHCALHits->Branch("hadZ", &flareZ, "hadZ/F"); + fFLArEHCALHits->Branch("hadPx", &flarePx, "hadPx/F"); // momentum + fFLArEHCALHits->Branch("hadPy", &flarePy, "hadPy/F"); + fFLArEHCALHits->Branch("hadPz", &flarePz, "hadPz/F"); + fFLArEHCALHits->Branch("hadDeltaPx", &flareDeltaPx, "hadDeltaPx/F"); + fFLArEHCALHits->Branch("hadDeltaPy", &flareDeltaPy, "hadDeltaPy/F"); + fFLArEHCALHits->Branch("hadDeltaPz", &flareDeltaPz, "hadDeltaPz/F"); + fFLArEHCALHits->Branch("hadEdep", &flareEdep, "hadEdep/F"); + fFLArEHCALHits->Branch("hadIsZX", &flareIsZX, "hadIsZX/B"); + + // FLArE pseudo-reco + if (fSavePseudoReco) + { + fFLArEPseudoReco = new TTree("pseudo_reco", "FLArE pseudo-reco info"); + fFLArEPseudoReco->Branch("evtID", &evtID, "evtID/I"); + fFLArEPseudoReco->Branch("edepInLAr", &edepInLAr, "edepInLAr/F"); + fFLArEPseudoReco->Branch("edepInHCAlX", &edepInHCALX, "edepInHCALX/F"); + fFLArEPseudoReco->Branch("edepInHCALY", &edepInHCALY, "edepInHCALY/F"); + fFLArEPseudoReco->Branch("sparseFractionMem", &sparseFractionMem, "sparseFractionMem/F"); + fFLArEPseudoReco->Branch("sparseFractionBins", &sparseFractionBins, "sparseFractionBins/F"); + fFLArEPseudoReco->Branch("TotalDedxLongitudinal", &TotalDedxLongitudinal); + fFLArEPseudoReco->Branch("TrueTotalDedxLongitudinal", &TrueTotalDedxLongitudinal); + fFLArEPseudoReco->Branch("nPrimaryParticle", &nprimaries, "nPrimaryParticle/I"); + fFLArEPseudoReco->Branch("primaryPDG", &primaryPDG); + fFLArEPseudoReco->Branch("primaryTrackLength", &primaryTrackLength); + fFLArEPseudoReco->Branch("primaryTrackLengthInTPC", &primaryTrackLengthInTPC); + fFLArEPseudoReco->Branch("ProngEInLAr", &ProngEInLAr); + fFLArEPseudoReco->Branch("ProngEInHadCal", &ProngEInHadCal); + fFLArEPseudoReco->Branch("ProngAngleToBeamDir", &ProngAngleToBeamDir); + fFLArEPseudoReco->Branch("ShowerLength", &ShowerLength); + fFLArEPseudoReco->Branch("ShowerLengthInLAr", &ShowerLengthInLAr); + fFLArEPseudoReco->Branch("ShowerWidth", &ShowerWidth); + fFLArEPseudoReco->Branch("ShowerWidthInLAr", &ShowerWidthInLAr); + fFLArEPseudoReco->Branch("ProngAvgdEdx", &ProngAvgdEdx); + fFLArEPseudoReco->Branch("ProngAvgdEdxInLAr", &ProngAvgdEdxInLAr); + fFLArEPseudoReco->Branch("dir_pol_x", &dir_pol_x); + fFLArEPseudoReco->Branch("dir_pol_y", &dir_pol_y); + fFLArEPseudoReco->Branch("dir_pol_z", &dir_pol_z); + fFLArEPseudoReco->Branch("dir_coc_x", &dir_coc_x); + fFLArEPseudoReco->Branch("dir_coc_y", &dir_coc_y); + fFLArEPseudoReco->Branch("dir_coc_z", &dir_coc_z); + } + + fFile->cd(); +} + +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::bookFASER2Trees() +{ + // create subdirectory in file + fFASER2Dir = fFile->mkdir("faser2","FASER2 output",kTRUE); + fFile->cd(fFASER2Dir->GetName()); + + //* Acts Hits Tree [i == unsigned int; F == float; l == Long unsigned 64 int] + fActsHitsTree = new TTree("hits", "ActsHitsTree"); + fActsHitsTree->Branch("event_id", &ActsHitsEventID, "event_id/i"); + fActsHitsTree->Branch("geometry_id", &ActsHitsGeometryID, "geometryid/l"); + fActsHitsTree->Branch("particle_id", &ActsHitsParticleID, "particle_id/l"); + fActsHitsTree->Branch("tx", &ActsHitsX, "tx/F"); + fActsHitsTree->Branch("ty", &ActsHitsY, "ty/F"); + fActsHitsTree->Branch("tz", &ActsHitsZ, "tz/F"); + fActsHitsTree->Branch("tt", &ActsHitsT, "tt/F"); + fActsHitsTree->Branch("tpx", &ActsHitsPx, "tpx/F"); + fActsHitsTree->Branch("tpy", &ActsHitsPy, "tpy/F"); + fActsHitsTree->Branch("tpz", &ActsHitsPz, "tpz/F"); + fActsHitsTree->Branch("te", &ActsHitsE, "tpe/F"); + fActsHitsTree->Branch("deltapx", &ActsHitsDeltaPx, "deltapx/F"); + fActsHitsTree->Branch("deltapy", &ActsHitsDeltaPy, "deltapy/F"); + fActsHitsTree->Branch("deltapz", &ActsHitsDeltaPz, "deltapz/F"); + fActsHitsTree->Branch("deltae", &ActsHitsDeltaE, "deltae/F"); + fActsHitsTree->Branch("index", &ActsHitsIndex, "index/I"); + fActsHitsTree->Branch("volume_id", &ActsHitsVolumeID, "volume_id/i"); + fActsHitsTree->Branch("boundary_id", &ActsHitsBoundaryID, "boundary_id/i"); + fActsHitsTree->Branch("layer_id", &ActsHitsLayerID, "layer_id/i"); + fActsHitsTree->Branch("approach_id", &ActsHitsApproachID, "approach_id/i"); + fActsHitsTree->Branch("sensitive_id", &ActsHitsSensitiveID, "sensitive_id/i"); + + //* Acts truth particle tree + fActsParticlesTree = new TTree("particles", "ActsParticlesTree"); + fActsParticlesTree->Branch("event_id", &ActsHitsEventID, "event_id/i"); + fActsParticlesTree->Branch("particle_id", &ActsParticlesParticleId); + fActsParticlesTree->Branch("particle_type", &ActsParticlesParticleType); + fActsParticlesTree->Branch("process", &ActsParticlesProcess); + fActsParticlesTree->Branch("vx", &ActsParticlesVx); + fActsParticlesTree->Branch("vy", &ActsParticlesVy); + fActsParticlesTree->Branch("vz", &ActsParticlesVz); + fActsParticlesTree->Branch("vt", &ActsParticlesVt); + fActsParticlesTree->Branch("px", &ActsParticlesPx); + fActsParticlesTree->Branch("py", &ActsParticlesPy); + fActsParticlesTree->Branch("pz", &ActsParticlesPz); + fActsParticlesTree->Branch("m", &ActsParticlesM); + fActsParticlesTree->Branch("q", &ActsParticlesQ); + fActsParticlesTree->Branch("eta", &ActsParticlesEta); + fActsParticlesTree->Branch("phi", &ActsParticlesPhi); + fActsParticlesTree->Branch("pt", &ActsParticlesPt); + fActsParticlesTree->Branch("p", &ActsParticlesP); + fActsParticlesTree->Branch("vertex_primary", &ActsParticlesVertexPrimary); + fActsParticlesTree->Branch("vertex_secondary", &ActsParticlesVertexSecondary); + fActsParticlesTree->Branch("particle", &ActsParticlesParticle); + fActsParticlesTree->Branch("generation", &ActsParticlesGeneration); + fActsParticlesTree->Branch("sub_particle", &ActsParticlesSubParticle); + fActsParticlesTree->Branch("e_loss", &ActsParticlesELoss); + fActsParticlesTree->Branch("total_x0", &ActsParticlesPathInX0); + fActsParticlesTree->Branch("total_l0", &ActsParticlesPathInL0); + fActsParticlesTree->Branch("number_of_hits", &ActsParticlesNumberOfHits); + fActsParticlesTree->Branch("outcome", &ActsParticlesOutcome); + + fFile->cd(); } -void AnalysisManager::BeginOfRun() { +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::BeginOfRun() +{ + G4cout << "Run has been started, preparing output" << G4endl; + + if (fFile) + delete fFile; + + // Preparing output file + fFile = new TFile(fFilename.c_str(), "RECREATE"); - G4cout<<"TTree is booked and the run has been started"<GetSDNamelist(); + G4cout << "Number of SDs : " << fSDNamelist.size() << G4endl; + + // Loop through active sentive detectors + // book trees as necessary + for (auto sdname : fSDNamelist) + { + G4cout << sdname.first << " " << sdname.second << G4endl; - fH5Filename = m_filename; - if(fH5Filename.find(".root") != std::string::npos) { - const size_t pos = fH5Filename.find(".root"); - fH5Filename.resize(pos); + // if FLArE is enabled, book trees + if ( fEnableFLArE && sdname.second.find("FLArE") != std::string::npos ) + { + fFlareSDs.push_back(sdname.first); + bookFLArETrees(); + } + + // if FASER2 is enabled, book ACTS trees + // but give the option to disable output if required + else if (fSaveActs && sdname.second.find("FASER2") != std::string::npos ) + { + fFaser2SDs.push_back(sdname.first); + bookFASER2Trees(); + } } - fH5Filename += ".h5"; - fH5file = hep_hpc::hdf5::File(fH5Filename, H5F_ACC_TRUNC); + // if FLArE is enabled & its geometry was found + // prepare .h5 output file + if( fFlareSDs.size()>0 ) + { + fH5Filename = fFilename; + if (fH5Filename.find(".root") != std::string::npos) + { + const size_t pos = fH5Filename.find(".root"); + fH5Filename.resize(pos); + } + fH5Filename += ".h5"; + fH5file = hep_hpc::hdf5::File(fH5Filename, H5F_ACC_TRUNC); + } - SDNamelist = GeometricalParameters::Get()->GetSDNamelist(); - NumberOfSDs = SDNamelist.size(); - G4cout<<"Number of SDs : "<cd(); - evt->Write(); - if(m_saveTrack) trk->Write(); - if (m_saveActs) acts_hits_tree->Write(); - if (m_saveActs) acts_particles_tree->Write(); - thefile->Close(); - fH5file.close(); +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::EndOfRun() +{ + // save common trees at the top of the output file + fFile->cd(); + fEvt->Write(); + fPrim->Write(); + if (fSaveTrack) fTrk->Write(); + + // save detector-specif trees in their directories + if (fFlareSDs.size()>0) + { + fFile->cd(fFLArEDir->GetName()); + fFLArEHits->Write(); + fFLArEHCALHits->Write(); + if(fSavePseudoReco) fFLArEPseudoReco->Write(); + fFile->cd(); // go back to top + fH5file.close(); + } + if (fFaser2SDs.size()>0) + { + fFile->cd(fFASER2Dir->GetName()); + fActsHitsTree->Write(); + fActsParticlesTree->Write(); + fFile->cd(); // go back to top + } + + fFile->Close(); } -void AnalysisManager::BeginOfEvent() { - neutrino = FPFNeutrino(); +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::BeginOfEvent() +{ + // TODO: REVIEW, CLEAN-UP UNUSED VARIABLES + // TODO: change arrays into std::vectors + // reset vectors that need to be cleared for a new event + // only reset arrays or vectors, no need for other defaults?? + + neutrino = FPFNeutrino(); // remove??? primaries.clear(); primaryIDs.clear(); - nHits = 0; - sparseFractionMem = -1; - sparseFractionBins = -1; - edepInLAr = 0; - edepInHadCalX = 0; - edepInHadCalY = 0; - edepInMuonFinderX = 0; - edepInMuonFinderY = 0; - edepInHadAborb = 0; - edepInMuonFinderAbsorb = 0; - missCountedEnergy = 0; - nPrimaryParticle = 0; - nTestNPrimaryTrack = 0; - nFromFSLParticles = 0; - nFromFSPizeroParticles = 0; - nFromFSLDecayPizeroParticles = 0; - countPrimaryParticle = 0; - for (G4int j= 0; j< 3000; ++j) { - TotalDedxLongitudinal[j] = 0; - TrueTotalDedxLongitudinal[j] = 0; - } - for (G4int i= 0; i< 1000; ++i) { - primaryParentPDG[i] = 0; - primaryTrackLength[i] = 0; - primaryTrackLengthInTPC[i] = 0; - ProngEInDetector[i] = 0; - ProngEInLAr[i] = 0; - ProngEInHadCal[i] = 0; - ProngEInMuonFinder[i] = 0; - ProngEInMuonFinderLayer1X[i] = 0; - ProngEInMuonFinderLayer1Y[i] = 0; - ProngEInMuonFinderLayer2X[i] = 0; - ProngEInMuonFinderLayer2Y[i] = 0; - ProngAngleToBeamDir[i] = -1; - ShowerLength[i] = -1; - ShowerLengthInLAr[i] = -1; - ShowerWidth[i] = 0; - ShowerWidthInLAr[i] = 0; - ProngAvgdEdx[i] = -1; - ProngAvgdEdxInLAr[i] = -1; - for (G4int j= 0; j< 100; ++j) { - ProngdEdxAlongTrack[i][j] = 0; - ProngdEdxTrackLength[i][j] = -1; - } - dir_pol_x[i] = -999; - dir_pol_y[i] = -999; - dir_pol_z[i] = -999; - dir_coc_x[i] = -999; - dir_coc_y[i] = -999; - dir_coc_z[i] = -999; - fromFSLParticlePDG[i] = 0; - } - if (m_saveHit) { - for (G4int i= 0; i< 40000000; ++i) { - HitTID[i] = -1; - HitPID[i] = -1; - HitPDG[i] = 0; - HitTrackStatus[i] = -1; - HitPrePositionX[i] = -999; - HitPrePositionY[i] = -999; - HitPrePositionZ[i] = -999; - HitPosPositionX[i] = -999; - HitPosPositionY[i] = -999; - HitPosPositionZ[i] = -999; - HitEdep[i] = 0; - } - } - // vectors that need to be cleared for a new event - allTracksPTPair.clear(); - trackClusters.clear(); - tracksFromFSL.clear(); - tracksFromFSLSecondary.clear(); - tracksFromFSPizeroSecondary.clear(); - tracksFromFSLDecayPizeroSecondary.clear(); - fPrimIdxFSL = -1; - - magzpos.clear(); - trkxc.clear(); - trkzc.clear(); - trkrc.clear(); - trkqIn.clear(); - trkmIn.clear(); - trkqOut.clear(); - trkmOut.clear(); - - hitXFSL.clear(); - hitYFSL.clear(); - hitZFSL.clear(); - hitPFSL.clear(); - trkXFSL.clear(); - trkYFSL.clear(); - trkZFSL.clear(); - trkPFSL.clear(); - - trackTID = -1; - trackPID = -1; - trackPDG = -1; - trackKinE = -1; - trackNPoints = -1; - trackPointX.clear(); - trackPointY.clear(); - trackPointZ.clear(); + // track ID to primary ancestor association + trackToPrimaryAncestor.clear(); - ActsHitsEventID = 0; - ActsHitsGeometryID = 0; - ActsHitsParticleID = 0; - ActsHitsX = 0; - ActsHitsY = 0; - ActsHitsZ = 0; - ActsHitsT = 0; - ActsHitsPx = 0; - ActsHitsPy = 0; - ActsHitsPz = 0; - ActsHitsE = 0; - ActsHitsDeltaPx = 0; - ActsHitsDeltaPy = 0; - ActsHitsDeltaPz = 0; - ActsHitsDeltaE = 0; - ActsHitsIndex = 0; - ActsHitsVolumeID = 0; - ActsHitsBoundaryID = 0; - ActsHitsLayerID = 0; - ActsHitsApproachID = 0; - ActsHitsSensitiveID = 0; + trackPointX.clear(); + trackPointY.clear(); + trackPointZ.clear(); ActsParticlesParticleId.clear(); ActsParticlesParticleType.clear(); @@ -426,272 +424,376 @@ void AnalysisManager::BeginOfEvent() { ActsParticlesPathInL0.clear(); ActsParticlesNumberOfHits.clear(); ActsParticlesOutcome.clear(); + + // these are used to accumulate + // so need to be reset + edepInLAr = 0; + edepInHCALX = 0; + edepInHCALY = 0; + + TotalDedxLongitudinal.clear(); + TrueTotalDedxLongitudinal.clear(); + TotalDedxLongitudinal.resize(3000, 0.0);; + TrueTotalDedxLongitudinal.resize(3000, 0.0); + + primaryPDG.clear(); + primaryTrackLength.clear(); + primaryTrackLengthInTPC.clear(); + ProngEInLAr.clear(); + ProngEInHadCal.clear(); + ProngAngleToBeamDir.clear(); + ShowerLength.clear(); + ShowerLengthInLAr.clear(); + ShowerWidth.clear(); + ShowerWidthInLAr.clear(); + ProngAvgdEdx.clear(); + ProngAvgdEdxInLAr.clear(); + dir_pol_x.clear(); + dir_pol_y.clear(); + dir_pol_z.clear(); + dir_coc_x.clear(); + dir_coc_y.clear(); + dir_coc_z.clear(); + } -void AnalysisManager::EndOfEvent(const G4Event* event) { +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::EndOfEvent(const G4Event *event) +{ + /// evtID evtID = event->GetEventID(); + // FILL EVENT TREE + // nothing else except evtID for now... + // TODO: need to find good way to propagate generator-level info + fEvt->Fill(); + + //----------------------------------------------------------- + + // FILL PRIMARIES/TRAJECTORIES TREE + FillPrimariesTree(event); + if(fSaveTrack) FillTrajectoriesTree(event); + + //----------------------------------------------------------- + + // Get the hit collections + // If there is no hit collection, there is nothing to be done + fHCofEvent = event->GetHCofThisEvent(); + if (!fHCofEvent) + { + G4cout << "No hits recorded in any sensitive volume --> nothing to save!" << G4endl; + return; + } + + //----------------------------------------------------------- + + // FILL DETECTOR HITS + + if( fFlareSDs.size() > 0 ) FillFLArEOutput(); + + if( fFaser2SDs.size() > 0 ) FillFASER2Output(); + +} + +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::FillPrimariesTree(const G4Event *event) +{ + nPrimaryVertex = event->GetNumberOfPrimaryVertex(); + G4cout << "\nNumber of primary vertices : " << nPrimaryVertex << G4endl; + /// loop over the vertices, and then over primary particles, /// neutrino truth info from event generator. - for (G4int ivtx = 0; ivtx < event->GetNumberOfPrimaryVertex(); ++ivtx) { - for (G4int ipp = 0; ipp < event->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); ++ipp) { - G4PrimaryParticle* primary_particle = - event->GetPrimaryVertex(ivtx)->GetPrimary(ipp); - if (primary_particle) { - PrimaryParticleInformation* primary_particle_info = - dynamic_cast(primary_particle->GetUserInformation()); + for (G4int ivtx = 0; ivtx < event->GetNumberOfPrimaryVertex(); ++ivtx) + { + G4cout << "=== Vertex " << ivtx+1 << " of " << nPrimaryVertex << " -> " + << event->GetPrimaryVertex(ivtx)->GetNumberOfParticle() << " primaries ===" << G4endl; + for (G4int ipp = 0; ipp < event->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); ++ipp) + { + G4PrimaryParticle *primary_particle = event->GetPrimaryVertex(ivtx)->GetPrimary(ipp); + if (primary_particle) + { + PrimaryParticleInformation *primary_particle_info = dynamic_cast(primary_particle->GetUserInformation()); primary_particle_info->Print(); - // add to list of primary particles - // this list is then appended further if necessary - countPrimaryParticle++; - primaryIDs.push_back(primary_particle_info->GetPartID()); //store to avoid duplicates - primaries.push_back(FPFParticle(primary_particle_info->GetPDG(), 0, - primary_particle_info->GetPartID(), - countPrimaryParticle-1, 1, - primary_particle_info->GetMass(), - primary_particle_info->GetVertexMC().x(), - primary_particle_info->GetVertexMC().y(), - primary_particle_info->GetVertexMC().z(), - 0., - primary_particle_info->GetMomentumMC().x(), - primary_particle_info->GetMomentumMC().y(), - primary_particle_info->GetMomentumMC().z(), - GetTotalEnergy(primary_particle_info->GetMomentumMC().x(),primary_particle_info->GetMomentumMC().y(), - primary_particle_info->GetMomentumMC().z(), primary_particle_info->GetMass()) - )); - - // if it's the first, also store the nuetrino info (if available) - if (primary_particle_info->GetPartID()==0) { - nuIdx = primary_particle_info->GetNeuIdx(); - nuPDG = primary_particle_info->GetNeuPDG(); - nuE = primary_particle_info->GetNeuP4().T(); - nuX = primary_particle_info->GetNeuX4().X(); - nuY = primary_particle_info->GetNeuX4().Y(); - nuZ = primary_particle_info->GetNeuX4().Z(); - nuIntType = primary_particle_info->GetInteractionTypeId(); - nuScatteringType = primary_particle_info->GetScatteringTypeId(); - nuW = primary_particle_info->GetW(); - nuFSLPDG = primary_particle_info->GetFSLPDG(); - nuFSLPx = primary_particle_info->GetFSLP4().X(); - nuFSLPy = primary_particle_info->GetFSLP4().Y(); - nuFSLPz = primary_particle_info->GetFSLP4().Z(); - nuFSLE = primary_particle_info->GetFSLP4().T(); - neutrino = FPFNeutrino(nuPDG, nuX, nuY, nuZ, 0, 0, 0, 0, nuE, - nuIdx, nuIntType, nuScatteringType, nuW, nuFSLPDG, nuFSLPx, nuFSLPy, nuFSLPz, nuFSLE); - continue; - } + primVtxID = ivtx; + primTrackID = primary_particle_info->GetPartID(); // confirm matches track's id later? + + auto particleId = ActsFatras::Barcode(); + particleId.setVertexPrimary(ivtx); + particleId.setGeneration(0); + particleId.setSubParticle(0); + particleId.setParticle(primTrackID - 1); + + primParticleID = particleId.value(); + primPDG = primary_particle_info->GetPDG(); + primVx = primary_particle_info->GetVertexMC().x(); + primVy = primary_particle_info->GetVertexMC().y(); + primVz = primary_particle_info->GetVertexMC().z(); + primVt = 0; + primPx = primary_particle_info->GetMomentumMC().x(); + primPy = primary_particle_info->GetMomentumMC().y(); + primPz = primary_particle_info->GetMomentumMC().z(); + primM = primary_particle_info->GetMass()/MeV; + primQ = primary_particle_info->GetCharge(); + + G4double energy = GetTotalEnergy(primPx, primPy, primPz, primM); + TLorentzVector p4(primPx,primPy,primPz,energy); + primEta = p4.Eta(); + primPhi = p4.Phi(); + primPt = p4.Pt(); + primP = p4.P(); + primE = energy; + primKE = energy - primM; + + // store a copy as a FPFParticle for further processing + primaryIDs.push_back(primTrackID); //store to avoid duplicates + primaries.push_back(FPFParticle(primPDG, 0, + primTrackID, primaryIDs.size()-1, 1, + primM, + primVx, primVy, primVz, primVt, + primPx, primPy, primPz,energy)); + fPrim->Fill(); } } } - nPrimaryVertex = event->GetNumberOfPrimaryVertex(); - std::cout<<"\nnumber of primary vertices : "<"<GetHCofThisEvent(); - if (!hcofEvent) return; - // loop over all sensitive detectors to get all the hit collections - // fill the truth tree for primary particles - for (auto sdname : SDNamelist) { - FillPrimaryTruthTree(sdname.first, sdname.second); + G4cout << "\nNumber of primaries : " << primaryIDs.size() << G4endl; +} + +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::FillTrajectoriesTree(const G4Event* event) +{ + int count_tracks = 0; + + G4cout << "==== Saving track information to tree ====" << G4endl; + auto trajectoryContainer = event->GetTrajectoryContainer(); + if (!trajectoryContainer) + { + G4cout << "No tracks found: did you enable their storage with '/tracking/storeTrajectory 1'?" << G4endl; + return; } - // update number of primary particles after FillPrimaryTruthTree - // including decay products from primary tau and pizero - nPrimaryParticle = countPrimaryParticle; // TODO: they're equivalent, should only store one of them - nFromFSLParticles = tracksFromFSLSecondary.size(); - nFromFSPizeroParticles = tracksFromFSPizeroSecondary.size(); - nFromFSLDecayPizeroParticles = tracksFromFSLDecayPizeroSecondary.size(); - - // find all the tracks originate from the final state lepton, include FSL itself (TID=1) - // should only work with neutrino interaction generator - // exception to single particle generator: tau, mu - if (neutrino.NuPDG()!=0 || abs(primaries[0].PDG())==15 || abs(primaries[0].PDG())==13) { - tracksFromFSL.insert(1); - for (auto x : allTracksPTPair) { - if (tracksFromFSL.find(x.first) != tracksFromFSL.end()) { - tracksFromFSL.insert(x.second); - } + + for (size_t i = 0; i < trajectoryContainer->entries(); ++i) + { + auto trajectory = static_cast((*trajectoryContainer)[i]); + trackTID = trajectory->GetTrackID(); + trackPID = trajectory->GetParentID(); + trackPDG = trajectory->GetPDGEncoding(); + trackKinE = trajectory->GetInitialKineticEnergy(); + trackNPoints = trajectory->GetPointEntries(); + count_tracks++; + for (size_t j = 0; j < trackNPoints; ++j) + { + G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition(); + trackPointX.push_back( pos.x() ); + trackPointY.push_back( pos.y() ); + trackPointZ.push_back( pos.z() ); } + fTrk->Fill(); + trackPointX.clear(); + trackPointY.clear(); + trackPointZ.clear(); } + G4cout << "Total number of recorded track: " << count_tracks << std::endl; +} + +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::FillFLArEOutput() +{ + G4cout << "===== Filling FLArE output trees =====" << G4endl; + + // prepare the pixel map + const double_t res_tpc[3] = {1, 5, 5}; // mm - // tracksFromFSL includes all the tracks orginating from the fsl - // tracksFromFSLSecondary only inclues the tracks directly decayed from the fsl - std::cout<<"Recorded tracks : "<SetPMBoundary(GeometricalParameters::Get()->GetFLArEPosition()/mm - + GeometricalParameters::Get()->GetTPCSizeXYZ()/mm/2, + GeometricalParameters::Get()->GetFLArEPosition()/mm + + GeometricalParameters::Get()->GetTPCSizeXYZ()/mm/2); + pm3D->InitializePM(); + + // accumulate values per primary particle + ProngEInLAr.resize(primaries.size(),0.); + ProngEInHadCal.resize(primaries.size(),0.); + ShowerLength.resize(primaries.size(),0.); + ShowerLengthInLAr.resize(primaries.size(),0.); + ShowerWidth.resize(primaries.size(),0.); + ShowerWidthInLAr.resize(primaries.size(),0.); + primaryTrackLength.resize(primaries.size(),0.); + primaryTrackLengthInTPC.resize(primaries.size(),0.); + + // loop over the detected FLArE sensitive volumes + // for now, everything is lar_box + int nLArHits = 0; + int nHCALHits = 0; + for(const int sdId : fFlareSDs ) + { + // Get and cast hit collection with LArBoxHits + std::string sdName = fSDNamelist.at(sdId); + auto hitCollection = dynamic_cast(fHCofEvent->GetHC(sdId)); + if (!hitCollection) + { + G4cout << "No hits recorded by " << sdName << G4endl; + continue; } - } - if (fPrimIdxFSL>=0) primaries[fPrimIdxFSL].SetProngType(0); - for (auto x : allTracksPTPair) { - // if this track is the fsl (TID=1) and it decays (nFromFSLParticles>0), - // then it forms a single cluster by itself, this is mainly for studying the tau decay. - if ((x.second==1) && (nFromFSLParticles>0)) continue; - // if this track is the decay product of the fsl, it should already been added to the trackClusters - if ((x.first==1) && (nFromFSLParticles>0) && (tracksFromFSLSecondary.find(x.second) != tracksFromFSLSecondary.end())) continue; - // if this is the decay product of the final state pizero, it should already been added to the trackClusters - if ((tracksFromFSPizeroSecondary.find(x.second) != tracksFromFSPizeroSecondary.end())) continue; - // if this is the decay product of the tau decay pizero, it should already been added to the trackClusters - if ((tracksFromFSLDecayPizeroSecondary.find(x.second) != tracksFromFSLDecayPizeroSecondary.end())) continue; - // add the track to the corresponding cluster if its parent is in the cluster. - // one track can have only one parent, break the loop once its parent is found. - for (int iPrim= 0; iPrim< nPrimaryParticle; ++iPrim) { - if (trackClusters[iPrim].find(x.first) != trackClusters[iPrim].end()) { - trackClusters[iPrim].insert(x.second); - break; - } - } - } + for (auto hit : *hitCollection->GetVector()) + { + flareTrackID = hit->GetTID(); + flareParentID = hit->GetPID(); + flarePDG = hit->GetPDG(); + flareCopyNum = hit->GetCopyNum(); + flareT = hit->GetTime(); + flareX = hit->GetPostStepPosition().x(); + flareY = hit->GetPostStepPosition().y(); + flareZ = hit->GetPostStepPosition().z(); + flarePx = hit->GetInitMomentum().x(); + flarePy = hit->GetInitMomentum().y(); + flarePz = hit->GetInitMomentum().z(); + flareDeltaPx = hit->GetDeltaMom().x(); + flareDeltaPy = hit->GetDeltaMom().y(); + flareDeltaPz = hit->GetDeltaMom().z(); + flareEdep = hit->GetEdep(); - if (GeometricalParameters::Get()->GetAddFLArE()) - { - G4cout << "Adding FLArE pixel map..." << G4endl; - const Double_t res_tpc[3] = {1, 5, 5}; // mm - if (neutrino.NuPDG()!=0) { - pm3D = new PixelMap3D(evtID, nPrimaryParticle, neutrino.NuPDG(), res_tpc); - } else { - pm3D = new PixelMap3D(evtID, nPrimaryParticle, primaries[0].PDG(), res_tpc); - } - // boundary in global coord. - pm3D->SetPMBoundary(GeometricalParameters::Get()->GetFLArEPosition()/mm - - GeometricalParameters::Get()->GetTPCSizeXYZ()/mm/2, - GeometricalParameters::Get()->GetFLArEPosition()/mm + - GeometricalParameters::Get()->GetTPCSizeXYZ()/mm/2); - pm3D->InitializePM(); - - /// FillTrueEdep must run after FillPrimaryTruthTree, - /// otherwise tracksFromFSL and tracksFromFSLSecondary are invalid - /// Pixel map is also filled here - for (auto sdname : SDNamelist) { - FillTrueEdep(sdname.first, sdname.second); - } + auto particleId = ActsFatras::Barcode(); + particleId.setVertexPrimary(1); // fix this value + particleId.setGeneration(flareParentID); + particleId.setSubParticle(0); + particleId.setParticle(flareTrackID); + flareParticleID = particleId.value(); - if (m_save2DEvd) pm3D->Write2DPMToFile(thefile); - - pm3D->Process3DPM(fH5file, neutrino, m_save3DEvd); - sparseFractionMem = pm3D->GetSparseFractionMem(); - sparseFractionBins = pm3D->GetSparseFractionBins(); - - if (m_circularFit){ - - circNhits = hitXFSL.size(); - trkNhits = trkXFSL.size(); - - // apply circular fitting to FLArE hadCath/muonFinder - if ( circNhits > 0 ){ - circularfitter::CircleFit* circFit = new circularfitter::CircleFit(hitXFSL,hitZFSL); - circularfitter::LineFit* lineFit = new circularfitter::LineFit(hitZFSL,hitYFSL); - xc = circFit->GetXc(); - zc = circFit->GetZc(); - rc = circFit->GetR(); - p0 = lineFit->GetP0(); - p1 = lineFit->GetP1(); - cosDip = lineFit->GetCosDip(); - } + double hit_position_xyz[3] = { flareX, flareY, flareZ }; + double vtx_xyz[3] = { primaries[0].Vx(), primaries[0].Vy(), primaries[0].Vz() }; - // apply circular fitting for FASER2 spectrometer magnet - if( trkNhits > 0 ){ - - Nmagnets = (GeometricalParameters::Get()->GetFASER2MagnetOption() == GeometricalParameters::magnetOption::SAMURAI) ? 1 : - GeometricalParameters::Get()->GetNFASER2Magnets(); - G4cout << "Number of FASER2 magnets: " << Nmagnets << G4endl; - - circularfitter::CircleExtractor* circExtract = new circularfitter::CircleExtractor(trkXFSL,trkYFSL,trkZFSL); - magzpos = circExtract->GetMagnetZs(); - trkxc = circExtract->GetXc(); - trkzc = circExtract->GetZc(); - trkrc = circExtract->GetR(); - - std::vector pre = circExtract->GetPreLine(); - std::vector post = circExtract->GetPostLine(); - for(int i=0; iGetP0(); - trkp1 = lineFit2->GetP1(); - trkcosDip = lineFit2->GetCosDip(); + // pseudo reco: track/shower length and width + primaryTrackLength[whichPrim] += hit->GetStepLength(); + double ShowerP = primaries[whichPrim].P(); + double dsquare_hit_vtx = TMath::Power((flareX-primaries[whichPrim].Vx()),2)+ + TMath::Power((flareY-primaries[whichPrim].Vy()),2)+ + TMath::Power((flareZ-primaries[whichPrim].Vz()),2); + double product_hit_p = (flareX-primaries[whichPrim].Vx())*primaries[whichPrim].Px()+ + (flareY-primaries[whichPrim].Vy())*primaries[whichPrim].Py()+ + (flareZ-primaries[whichPrim].Vz())*primaries[whichPrim].Pz(); + double len_hit = TMath::Abs(product_hit_p)/ShowerP; + double width_hit = TMath::Sqrt((dsquare_hit_vtx - product_hit_p*product_hit_p/ShowerP/ShowerP)); + ShowerLength[whichPrim] = (ShowerLength[whichPrim]>len_hit) ? ShowerLength[whichPrim] : len_hit; + double weighted_width_hit = width_hit*flareEdep; + if (!std::isnan(weighted_width_hit)) ShowerWidth[whichPrim] += weighted_width_hit; + + if (sdName == "FLArEBoxSD/lar_box") + { + nLArHits++; + fFLArEHits->Fill(); + + if (fAddDiffusion == "toy") + pm3D->FillEntryWithToyElectronTransportation(hit_position_xyz, vtx_xyz, flareEdep, whichPrim); + else if (fAddDiffusion == "single") + pm3D->FillEntryWithToySingleElectronTransportation(hit_position_xyz, vtx_xyz, flareEdep, whichPrim); + else if (sdName == "lArBoxSD/lar_box") + pm3D->FillEntry(hit_position_xyz, vtx_xyz, flareEdep, whichPrim); + + // accumulate Edep in LAr + edepInLAr += flareEdep; + ProngEInLAr[whichPrim] += flareEdep; + primaryTrackLengthInTPC[whichPrim] += hit->GetStepLength(); + ShowerLengthInLAr[whichPrim] = (ShowerLengthInLAr[whichPrim]>len_hit) ? ShowerLengthInLAr[whichPrim] : len_hit; + if (!std::isnan(weighted_width_hit)) ShowerWidthInLAr[whichPrim] += weighted_width_hit; + + // accumulate dE/dx in LAr + float_t longitudinal_distance_to_vtx = ((flareX-vtx_xyz[0])*primaries[0].Px()+ + (flareY-vtx_xyz[1])*primaries[0].Py()+ + (flareZ-vtx_xyz[2])*primaries[0].Pz())/primaries[0].P(); + if (longitudinal_distance_to_vtx>=0 && longitudinal_distance_to_vtx<3000) // within 3000 mm + TrueTotalDedxLongitudinal[Int_t(longitudinal_distance_to_vtx)] += hit->GetEdep(); + } + else if (sdName == "FLArEHadCalXSD/lar_box" || + sdName == "FLArEMuonFinderXSD/lar_box" || + sdName == "FLArEBabyMINDHorBarSD/lar_box" ) + { + // specify this is an ZX hit + nHCALHits++; + flareIsZX = true; + fFLArEHCALHits->Fill(); + + // accumulate Edep in HCAL + edepInHCALX += flareEdep; + ProngEInHadCal[whichPrim] += flareEdep; + } + else if (sdName == "FLArEHadCalYSD/lar_box" || + sdName == "FLArEMuonFinderYSD/lar_box" || + sdName == "FLArEBabyMINDVerBarSD/lar_box" || + sdName == "FLArEHadAbsorbSD/lar_box" || + sdName == "FLArEMuonFinderAbsorbSD/lar_box" ) + { + nHCALHits++; + flareIsZX = false; + fFLArEHCALHits->Fill(); + + // accumulate Edep in HCAL + edepInHCALY += flareEdep; + ProngEInHadCal[whichPrim] += flareEdep; } } - - // FillPseudoRecoVar must run after FillTrueEdep, otherwise some of the variables won't be filled - FillPseudoRecoVar(); - - delete pm3D; } - int count_tracks = 0; - if(m_saveTrack){ - G4cout << "---> Saving track information to tree..." << G4endl; - auto trajectoryContainer = event->GetTrajectoryContainer(); - if (trajectoryContainer) { - for (size_t i = 0; i < trajectoryContainer->entries(); ++i) { - auto trajectory = static_cast((*trajectoryContainer)[i]); - trackTID = trajectory->GetTrackID(); - trackPID = trajectory->GetParentID(); - trackPDG = trajectory->GetPDGEncoding(); - trackKinE = trajectory->GetInitialKineticEnergy(); - trackNPoints = trajectory->GetPointEntries(); - count_tracks++; - for (size_t j = 0; j < trackNPoints; ++j) { - G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition(); - trackPointX.push_back( pos.x() ); - trackPointY.push_back( pos.y() ); - trackPointZ.push_back( pos.z() ); - } - trk->Fill(); - trackPointX.clear(); - trackPointY.clear(); - trackPointZ.clear(); - } - } else G4cout << "No tracks found: did you enable their storage with '/tracking/storeTrajectory 1'?" << G4endl; - G4cout << "---> Done!" << G4endl; - } - - evt->Fill(); + if (fSave2DEvd) pm3D->Write2DPMToFile(fFile,fFLArEDir); - G4cout << "Total number of recorded hits : " << nHits << std::endl; - if(m_saveTrack) G4cout << "Total number of recorded track: " << count_tracks << std::endl; + pm3D->Process3DPM(fH5file, neutrino, fSave3DEvd); - for (int iPrim= 0; iPrim< nPrimaryParticle; ++iPrim) { - trackClusters[iPrim].clear(); - } - trackClusters.clear(); - trackClusters.shrink_to_fit(); - + sparseFractionMem = pm3D->GetSparseFractionMem(); + sparseFractionBins = pm3D->GetSparseFractionBins(); + + if( fSavePseudoReco ) FillFLArEPseudoReco(); + + G4cout << "Total FLArE recorded hits: " << nLArHits << G4endl; + G4cout << "Total FLArE HCAL recorded hits: " << nHCALHits << G4endl; + + delete pm3D; } -void AnalysisManager::FillPrimaryTruthTree(G4int sdId, std::string sdName) +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + +void AnalysisManager::FillFASER2Output() { + G4cout << "==== Filling FASER2 output trees ====" << G4endl; - if (m_saveActs) + // loop over the detected FASER2 sensitive volumes + int nHits = 0; + for(const int sdId : fFaser2SDs ) { - auto hitCollection = dynamic_cast(hcofEvent->GetHC(sdId)); - if (!hitCollection) return; + // Get and cast hit collection with LArBoxHits + std::string sdName = fSDNamelist.at(sdId); + auto hitCollection = dynamic_cast(fHCofEvent->GetHC(sdId)); + if (!hitCollection) + { + G4cout << "No hits recorded by " << sdName << G4endl; + continue; + } + std::map sub_part_map{}; - for (auto hit: *hitCollection->GetVector()) + for (auto hit : *hitCollection->GetVector()) { - if (hit->GetCharge() == 0) continue; // skip neutral particles, they don't hit + if (hit->GetCharge() == 0) + continue; // skip neutral particles, they don't hit /* * A note on the ActsHitsGeometryID variable @@ -701,9 +803,10 @@ void AnalysisManager::FillPrimaryTruthTree(G4int sdId, std::string sdName) As a result I set `geometry_id` to zero and give the the resposibility of assigning this variable to the user during the reading of the `hits` tree. */ + nHits++; ActsHitsEventID = evtID; ActsHitsGeometryID = 0; - + int hitID = hit->GetTrackID(); int nPrimaries = ActsParticlesParticleId.size(); @@ -713,10 +816,10 @@ void AnalysisManager::FillPrimaryTruthTree(G4int sdId, std::string sdName) particleId.setParticle(hit->GetTrackID() - 1); // The track ID is the primary particle index plus one particleId.setGeneration(hit->GetParentID()); - sub_part_map.try_emplace(hit->GetTrackID()-1, sub_part_map.size()); - + sub_part_map.try_emplace(hit->GetTrackID() - 1, sub_part_map.size()); + // This is a fudge - assumes that that the secondary particles are always sub-particles of the primary particle - particleId.setSubParticle(hit->GetParentID() == 0 ? 0 : sub_part_map[hit->GetTrackID()-1]); + particleId.setSubParticle(hit->GetParentID() == 0 ? 0 : sub_part_map[hit->GetTrackID() - 1]); ActsHitsParticleID = particleId.value(); ActsHitsX = hit->GetX(); @@ -737,20 +840,22 @@ void AnalysisManager::FillPrimaryTruthTree(G4int sdId, std::string sdName) // In principle with the right headers from Acts we could construct the geometry ID value here ActsHitsVolumeID = 1; ActsHitsBoundaryID = 0; - ActsHitsLayerID = (hit->GetCopyNumSensor()+1)*2; // Acts specfic layer ID, goes 2, 4, 6, ... + ActsHitsLayerID = (hit->GetCopyNumSensor() + 1) * 2; // Acts specfic layer ID, goes 2, 4, 6, ... ActsHitsApproachID = 0; ActsHitsSensitiveID = 1; - acts_hits_tree->Fill(); + fActsHitsTree->Fill(); // Now fill the Acts particles tree bool isDuplicate = false; - for (const auto& id : ActsParticlesParticleId) { - if (id == particleId.value()) { + for (const auto &id : ActsParticlesParticleId) + { + if (id == particleId.value()) + { isDuplicate = true; } } if (isDuplicate) continue; // Skip this particle if it's already been added - + ActsParticlesParticleId.push_back(particleId.value()); ActsParticlesParticleType.push_back(hit->GetPDGID()); ActsParticlesProcess.push_back(0); @@ -763,394 +868,108 @@ void AnalysisManager::FillPrimaryTruthTree(G4int sdId, std::string sdName) ActsParticlesPz.push_back(hit->GetTrackP4().pz()); ActsParticlesM.push_back(hit->GetTrackP4().m()); ActsParticlesQ.push_back(hit->GetCharge()); - + ActsParticlesEta.push_back(hit->GetTrackP4().eta()); ActsParticlesPhi.push_back(hit->GetTrackP4().phi()); ActsParticlesPt.push_back(pow(pow(hit->GetTrackP4().px(), 2) + pow(hit->GetTrackP4().py(), 2), 0.5)); ActsParticlesP.push_back(pow(pow(hit->GetTrackP4().px(), 2) + pow(hit->GetTrackP4().py(), 2) + pow(hit->GetTrackP4().pz(), 2), 0.5)); - ActsParticlesVertexPrimary.push_back(hit->GetIsPrimaryTrack()); //? These variables need to be filled, but are unused by Acts - ActsParticlesVertexSecondary.push_back(hit->GetIsSecondaryTrack()); //? These variables need to be filled, but are unused by Acts - ActsParticlesParticle.push_back(1); //? These variables need to be filled, but are unused by Acts - ActsParticlesGeneration.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesSubParticle.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesELoss.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesPathInX0.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesPathInL0.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesNumberOfHits.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesOutcome.push_back(0); //? These variables need to be filled, but are unused by Acts + ActsParticlesVertexPrimary.push_back(hit->GetIsPrimaryTrack()); //? These variables need to be filled, but are unused by Acts + ActsParticlesVertexSecondary.push_back(hit->GetIsSecondaryTrack()); //? These variables need to be filled, but are unused by Acts + ActsParticlesParticle.push_back(1); //? These variables need to be filled, but are unused by Acts + ActsParticlesGeneration.push_back(0); //? These variables need to be filled, but are unused by Acts + ActsParticlesSubParticle.push_back(0); //? These variables need to be filled, but are unused by Acts + ActsParticlesELoss.push_back(0); //? These variables need to be filled, but are unused by Acts + ActsParticlesPathInX0.push_back(0); //? These variables need to be filled, but are unused by Acts + ActsParticlesPathInL0.push_back(0); //? These variables need to be filled, but are unused by Acts + ActsParticlesNumberOfHits.push_back(0); //? These variables need to be filled, but are unused by Acts + ActsParticlesOutcome.push_back(0); //? These variables need to be filled, but are unused by Acts } // end of loop over hits - acts_particles_tree->Fill(); + fActsParticlesTree->Fill(); } - - // Get and cast hit collection with LArBoxHits - // LArBoxHitsCollection* hitCollection = dynamic_cast(hcofEvent->GetHC(sdId)); - // if (hitCollection) { - // for (auto hit: *hitCollection->GetVector()) { - // nHits++; - - // double pre_x = hit->GetPreStepPosition().x(); - // double pre_y = hit->GetPreStepPosition().y(); - // double pre_z = hit->GetPreStepPosition().z(); - // double post_x = hit->GetPostStepPosition().x(); - // double post_y = hit->GetPostStepPosition().y(); - // double post_z = hit->GetPostStepPosition().z(); - - // if (m_saveHit) { - // if (nHits<=40000000) { - // HitTID[nHits-1] = hit->GetTID(); - // HitPID[nHits-1] = hit->GetPID(); - // HitPDG[nHits-1] = hit->GetParticle(); - // HitTrackStatus[nHits-1] = hit->GetTrackStatus(); - // HitPrePositionX[nHits-1] = pre_x; - // HitPrePositionY[nHits-1] = pre_y; - // HitPrePositionZ[nHits-1] = pre_z; - // HitPosPositionX[nHits-1] = post_x; - // HitPosPositionY[nHits-1] = post_y; - // HitPosPositionZ[nHits-1] = post_z; - // HitEdep[nHits-1] = hit->GetEdep(); - // } - // } - - // // energy deposition in different volumes of the detector - // if (sdName == "lArBoxSD/lar_box") - // edepInLAr += hit->GetEdep(); - // else if (sdName == "HadCalXSD/lar_box") - // edepInHadCalX += hit->GetEdep(); - // else if (sdName == "HadCalYSD/lar_box") - // edepInHadCalY += hit->GetEdep(); - // else if (sdName == "MuonFinderXSD/lar_box") - // edepInMuonFinderX += hit->GetEdep(); - // else if (sdName == "MuonFinderYSD/lar_box") - // edepInMuonFinderY += hit->GetEdep(); - // else if (sdName == "HadAbsorbSD/lar_box") - // edepInHadAborb += hit->GetEdep(); - // else if (sdName == "MuonFinderAbsorbSD/lar_box") - // edepInMuonFinderAbsorb += hit->GetEdep(); - - // // save FSL (only muons!) hits for circle fitting - // // check that particle is a primary muon - // if( TMath::Abs(hit->GetParticle())==13 && hit->GetPID() == 0 ){ - // double px = hit->GetInitMomentum().x(); - // double pz = hit->GetInitMomentum().z(); - // double p_perp = TMath::Sqrt(px*px+pz*pz); - // // save hits in FLArE hadCather+muonFinder - // if( (sdName == "HadCalXSD/lar_box") || (sdName == "HadCalYSD/lar_box") || - // (sdName == "MuonFinderXSD/lar_box") || (sdName == "MuonFinderYSD/lar_box") || - // (sdName == "HadAbsorbSD/lar_box") || (sdName == "MuonFinderAbsorbSD/lar_box") || - // (sdName == "BabyMINDHorBarSD/lar_box") || (sdName == "BabyMINDVerBarSD/lar_box")){ - // hitXFSL.push_back(post_x); - // hitYFSL.push_back(post_y); - // hitZFSL.push_back(post_z); - // hitPFSL.push_back(p_perp); - // } - // else if ((sdName == "TrkHorScinSD/lar_box") || (sdName == "TrkVerScinSD/lar_box")){ - // trkXFSL.push_back(post_x); - // trkYFSL.push_back(post_y); - // trkZFSL.push_back(post_z); - // trkPFSL.push_back(p_perp); - // } - // } - - // //allTracksPTPair.insert(std::make_pair(hit->GetPID(), hit->GetTID())); - - // // stable final state particles in GENIE, primary particles in Geant4 - // if (hit->GetCreatorProcess()=="PrimaryParticle") { // i.e. PID==0 - // if ( std::find(primaryIDs.begin(), primaryIDs.end(), hit->GetTID()) == primaryIDs.end() ) { - // // the following line excludes final state lepton tau from the primary particle list - // //if (abs(nuPDG)==16 && abs(nuFSLPDG)==15 && abs(hit->GetParticle()==15)) continue; - // countPrimaryParticle++; - // primaryIDs.push_back(hit->GetTID()); - // primaries.push_back(FPFParticle(hit->GetParticle(), - // hit->GetPID(), hit->GetTID(), countPrimaryParticle-1, 1, hit->GetParticleMass(), - // hit->GetTrackVertex().x(), hit->GetTrackVertex().y(), hit->GetTrackVertex().z(), 0, - // hit->GetInitMomentum().x(), hit->GetInitMomentum().y(), hit->GetInitMomentum().z(), - // GetTotalEnergy(hit->GetInitMomentum().x(), hit->GetInitMomentum().y(), - // hit->GetInitMomentum().z(), hit->GetParticleMass()))); - // } - // } - // // in case of the fsl decay, the decay products are counted as primary particles - // // * tau- decay (dominant) - // // * mu- decay - // //if (hit->GetPID()==1 && hit->GetCreatorProcess()=="Decay") { - // if (hit->GetIsTrackFromPrimaryLepton()) { - // tracksFromFSLSecondary.insert(hit->GetTID()); - // if (std::find(primaryIDs.begin(), primaryIDs.end(), hit->GetTID()) == primaryIDs.end()) { - // countPrimaryParticle++; - // primaryIDs.push_back(hit->GetTID()); - // primaries.push_back(FPFParticle(hit->GetParticle(), - // hit->GetPID(), hit->GetTID(), countPrimaryParticle-1, 2, hit->GetParticleMass(), - // hit->GetTrackVertex().x(), hit->GetTrackVertex().y(), hit->GetTrackVertex().z(), 0, - // hit->GetInitMomentum().x(), hit->GetInitMomentum().y(), hit->GetInitMomentum().z(), - // GetTotalEnergy(hit->GetInitMomentum().x(), hit->GetInitMomentum().y(), - // hit->GetInitMomentum().z(), hit->GetParticleMass()))); - // } - // } - // // in case of pizero in the list of primary track - // // its decay products are also counted as primary particles, mostly 2 gammas - // if (hit->GetIsTrackFromPrimaryPizero()) { - // tracksFromFSPizeroSecondary.insert(hit->GetTID()); - // if (std::find(primaryIDs.begin(), primaryIDs.end(), hit->GetTID()) == primaryIDs.end()) { - // countPrimaryParticle++; - // primaryIDs.push_back(hit->GetTID()); - // primaries.push_back(FPFParticle(hit->GetParticle(), - // hit->GetPID(), hit->GetTID(), countPrimaryParticle-1, 3, hit->GetParticleMass(), - // hit->GetTrackVertex().x(), hit->GetTrackVertex().y(), hit->GetTrackVertex().z(), 0, - // hit->GetInitMomentum().x(), hit->GetInitMomentum().y(), hit->GetInitMomentum().z(), - // GetTotalEnergy(hit->GetInitMomentum().x(), hit->GetInitMomentum().y(), - // hit->GetInitMomentum().z(), hit->GetParticleMass()))); - // } - // } - // // in case of tau decay pizero - // // decay products of this pizero are also counted as primary particles, mostly 2 gammas - // if (hit->GetIsTrackFromFSLPizero()) { - // tracksFromFSLDecayPizeroSecondary.insert(hit->GetTID()); - // if (std::find(primaryIDs.begin(), primaryIDs.end(), hit->GetTID()) == primaryIDs.end()) { - // countPrimaryParticle++; - // primaryIDs.push_back(hit->GetTID()); - // primaries.push_back(FPFParticle(hit->GetParticle(), - // hit->GetPID(), hit->GetTID(), countPrimaryParticle-1, 4, hit->GetParticleMass(), - // hit->GetTrackVertex().x(), hit->GetTrackVertex().y(), hit->GetTrackVertex().z(), 0, - // hit->GetInitMomentum().x(), hit->GetInitMomentum().y(), hit->GetInitMomentum().z(), - // GetTotalEnergy(hit->GetInitMomentum().x(), hit->GetInitMomentum().y(), - // hit->GetInitMomentum().z(), hit->GetParticleMass()))); - // } - // } - // } // end of hit loop - // } + G4cout << "Total FASER2 recorded hits: " << nHits << G4endl; } -void AnalysisManager::FillTrueEdep(G4int sdId, std::string sdName) +float_t AnalysisManager::GetTotalEnergy(float_t px, float_t py, float_t pz, float_t m) { - std::map map_tracksFromFSLSecondary; - int _idx = 0; - for (auto _tid : tracksFromFSLSecondary) { - map_tracksFromFSLSecondary[_tid] = _idx; - _idx++; - } - - LArBoxHitsCollection* hitCollection = dynamic_cast(hcofEvent->GetHC(sdId)); - if (hitCollection) { - for (auto hit: *hitCollection->GetVector()) { - // Particles decay from the final state lepton in GENIE, or decay from the primary particles in G4 - if (tracksFromFSLSecondary.find(hit->GetTID()) != tracksFromFSLSecondary.end()) { - int whichTrackFromFSL = map_tracksFromFSLSecondary[hit->GetTID()]; - if (hit->GetStepNo()==1) { - G4cout<<"TID : "<GetTID() <<", PID : " <GetPID() - <<", PDG : "<GetParticle() <<", CreatorProcess : "<GetCreatorProcess() - <<", Ek : " <GetInitKinEnergy()<<" MeV" <GetParticle(); - } - } - - int whichPrim = -1; - for (int iPrim= 0; iPrim< nPrimaryParticle; ++iPrim) { - if (trackClusters[iPrim].find(hit->GetTID()) != trackClusters[iPrim].end()) { - whichPrim = iPrim; - break; - } - } - if (whichPrim< 0) { - if (sdName == "lArBoxSD/lar_box") { - std::cout<<"Can't find the primary particle of the hit in TPC volume, something is wrong " - <GetParticle()<<" edep("<GetEdep()<<") TID("<GetTID()<<") PID(" - <GetPID()<<") creator process ("<GetCreatorProcess()<<") position" - <GetEdepPosition()<GetEdep(); - continue; - } - } - - double pos_x = hit->GetEdepPosition().x(); - double pos_y = hit->GetEdepPosition().y(); - double pos_z = hit->GetEdepPosition().z(); - double hit_position_xyz[3] = {pos_x, pos_y, pos_z}; - double vtx_xyz[3]; - if (neutrino.NuPDG()!=0) { - vtx_xyz[0] = neutrino.NuVx(); - vtx_xyz[1] = neutrino.NuVy(); - vtx_xyz[2] = neutrino.NuVz(); - } - else { - vtx_xyz[0] = primaries[0].Vx(); - vtx_xyz[1] = primaries[0].Vy(); - vtx_xyz[2] = primaries[0].Vz(); - } - - if ((sdName == "lArBoxSD/lar_box") && (m_addDiffusion == "toy")) { - pm3D->FillEntryWithToyElectronTransportation(hit_position_xyz, vtx_xyz, hit->GetEdep(), whichPrim); - } else if ((sdName == "lArBoxSD/lar_box") && (m_addDiffusion == "single")) { - pm3D->FillEntryWithToySingleElectronTransportation(hit_position_xyz, vtx_xyz, hit->GetEdep(), whichPrim); - } else if (sdName == "lArBoxSD/lar_box") { - pm3D->FillEntry(hit_position_xyz, vtx_xyz, hit->GetEdep(), whichPrim); - } - - if (sdName == "lArBoxSD/lar_box") { - double longitudinal_distance_to_vtx; // in mm - if (neutrino.NuPDG()!=0) { - longitudinal_distance_to_vtx = (pos_z-vtx_xyz[2]); - } else { - longitudinal_distance_to_vtx = ((pos_x-vtx_xyz[0])*primaries[0].Px()+ - (pos_y-vtx_xyz[1])*primaries[0].Py()+ - (pos_z-vtx_xyz[2])*primaries[0].Pz())/primaries[0].P(); - } - if (Int_t(longitudinal_distance_to_vtx)>=0 && Int_t(longitudinal_distance_to_vtx)<3000) { // within 3000 mm - TrueTotalDedxLongitudinal[Int_t(longitudinal_distance_to_vtx)] += hit->GetEdep(); - } - } - // calculate dEdx along the track - // combine the tracks if they come from the final state lepton, namely tau- - double ShowerP = primaries[whichPrim].P(); - if ((hit->GetPID()==0) | - (tracksFromFSLSecondary.find(hit->GetTID()) != tracksFromFSLSecondary.end()) | - (tracksFromFSPizeroSecondary.find(hit->GetTID()) != tracksFromFSPizeroSecondary.end()) | - (tracksFromFSLDecayPizeroSecondary.find(hit->GetTID()) != tracksFromFSLDecayPizeroSecondary.end())) { - primaryTrackLength[whichPrim] += hit->GetStepLength(); - if (sdName == "lArBoxSD/lar_box") { - primaryTrackLengthInTPC[whichPrim] += hit->GetStepLength(); - if ((hit->GetPID()==0) | - (tracksFromFSPizeroSecondary.find(hit->GetTID()) != tracksFromFSPizeroSecondary.end())) { - double longitudinal_distance_to_vtx = ((pos_x-primaries[whichPrim].Vx())*primaries[whichPrim].Px()+ - (pos_y-primaries[whichPrim].Vy())*primaries[whichPrim].Py()+ - (pos_z-primaries[whichPrim].Vz())*primaries[whichPrim].Pz())/ShowerP; - if (int(longitudinal_distance_to_vtx)>=0 && int(longitudinal_distance_to_vtx)<100) { // within 100 mm - ProngdEdxAlongTrack[whichPrim][int(longitudinal_distance_to_vtx)] += hit->GetEdep(); - ProngdEdxTrackLength[whichPrim][int(longitudinal_distance_to_vtx)] = int(longitudinal_distance_to_vtx); - } - } else { - if (fPrimIdxFSL>=0) { - double ShowerP_FSL = primaries[fPrimIdxFSL].P(); - double longitudinal_distance_to_vtx = ((pos_x-primaries[fPrimIdxFSL].Vx())*primaries[fPrimIdxFSL].Px()+ - (pos_y-primaries[fPrimIdxFSL].Vy())*primaries[fPrimIdxFSL].Py()+ - (pos_z-primaries[fPrimIdxFSL].Vz())*primaries[fPrimIdxFSL].Pz())/ShowerP_FSL; - if (int(longitudinal_distance_to_vtx)>=0 && int(longitudinal_distance_to_vtx)<100) { // within 100 mm - ProngdEdxAlongTrack[fPrimIdxFSL][int(longitudinal_distance_to_vtx)] += hit->GetEdep(); - ProngdEdxTrackLength[fPrimIdxFSL][int(longitudinal_distance_to_vtx)] = int(longitudinal_distance_to_vtx); - } - } - } - } - } - // calculate the shower/track length and width - // length: defined as the longest projection distance at the true direction between vertex and hits - // length = |\vector{hit_position}|\cdot\cos(theta) = \vertor{hit_position}\cdot\vector{P} / |\vector{P}| - // width: defined as the weighted average of the least distance of the hits to the true direction - double dsquare_hit_vtx = TMath::Power((pos_x-primaries[whichPrim].Vx()),2)+ - TMath::Power((pos_y-primaries[whichPrim].Vy()),2)+ - TMath::Power((pos_z-primaries[whichPrim].Vz()),2); - double product_hit_p = (pos_x-primaries[whichPrim].Vx())*primaries[whichPrim].Px()+ - (pos_y-primaries[whichPrim].Vy())*primaries[whichPrim].Py()+ - (pos_z-primaries[whichPrim].Vz())*primaries[whichPrim].Pz(); - double len_hit = TMath::Abs(product_hit_p)/ShowerP; - double width_hit = TMath::Sqrt((dsquare_hit_vtx - product_hit_p*product_hit_p/ShowerP/ShowerP)); - // exclude zero hit when calculating showerlength of the primary particle - // exclude hits from the cryo gap (detID=8) - if (hit->GetEdep()>0 && - (sdName=="lArBoxSD/lar_box" || sdName=="HadAbsorbSD/lar_box" || sdName=="MuonFinderAbsorbSD/lar_box")) { - ProngEInDetector[whichPrim] += hit->GetEdep(); - ShowerLength[whichPrim] = std::max({ShowerLength[whichPrim], len_hit}); - //double square_weighted_width_hit = TMath::Power(width_hit*hit->GetEdep(),2); - double weighted_width_hit = width_hit*hit->GetEdep(); - if (!std::isnan(weighted_width_hit)) ShowerWidth[whichPrim] += weighted_width_hit; - if (sdName=="lArBoxSD/lar_box") { - ProngEInLAr[whichPrim] += hit->GetEdep(); - ShowerLengthInLAr[whichPrim] = std::max({ShowerLengthInLAr[whichPrim], len_hit}); - if (!std::isnan(weighted_width_hit)) ShowerWidthInLAr[whichPrim] += weighted_width_hit; - } - } - if ((sdName=="HadCalXSD/lar_box") || - (sdName=="HadCalYSD/lar_box") || - (sdName=="HadAbsorbSD/lar_box")) { - ProngEInHadCal[whichPrim] += hit->GetEdep(); - } - if ((sdName=="MuonFinderXSD/lar_box") || - (sdName=="MuonFinderYSD/lar_box") || - (sdName=="MuonFinderAbsorbSD/lar_box")) { - ProngEInMuonFinder[whichPrim] += hit->GetEdep(); - if (sdName=="MuonFinderXSD/lar_box") { - if (pos_z < GeometricalParameters::Get()->GetFLArEPosition().z()/mm+ - GeometricalParameters::Get()->GetTPCInsulationThickness()/mm+ - GeometricalParameters::Get()->GetHadCalLength()/mm+ - GeometricalParameters::Get()->GetMuonCatcherLength()/mm/2) { - ProngEInMuonFinderLayer1X[whichPrim] += hit->GetEdep(); - } else { - ProngEInMuonFinderLayer2X[whichPrim] += hit->GetEdep(); - } - } else if (sdName=="MuonFinderYSD/lar_box") { - if (pos_z < GeometricalParameters::Get()->GetFLArEPosition().z()/mm+ - GeometricalParameters::Get()->GetTPCInsulationThickness()/mm+ - GeometricalParameters::Get()->GetHadCalLength()/mm+ - GeometricalParameters::Get()->GetMuonCatcherLength()/mm/2) { - ProngEInMuonFinderLayer1Y[whichPrim] += hit->GetEdep(); - } else { - ProngEInMuonFinderLayer2Y[whichPrim] += hit->GetEdep(); - } - } - } - - } // end of hit loop - } -} - -double AnalysisManager::GetTotalEnergy(double px, double py, double pz, double m) { - return TMath::Sqrt(px*px+py*py+pz*pz+m*m); + return TMath::Sqrt(px * px + py * py + pz * pz + m * m); } -void AnalysisManager::FillPseudoRecoVar() { +void AnalysisManager::FillFLArEPseudoReco() +{ // AngleToBeamDir, dEdx, dEdxInLAr ProngType - std::cout<0) { - ShowerWidth[iPrim] = ShowerWidth[iPrim]/ProngEInDetector[iPrim]; + std::cout << std::fixed << std::setw(10) << "PDG" + << std::setw(12) << "Angle" + << std::setw(13) << "TrackLength" + << std::setw(13) << "ShowerLength" + << std::setw(18) << "ShowerWidthInLAr" + << std::setw(12) << "EInLAr" + << std::setw(12) << "EInHadCal" + << std::setw(12) << "dEdxInLAr" + << std::setw(10) << "ProngType" + << std::setw(12) << "Pz" << std::endl; + + nprimaries = primaries.size(); + primaryPDG.resize(nprimaries,0.); + ProngAngleToBeamDir.resize(nprimaries,0.); + ProngAvgdEdx.resize(nprimaries,0.); + ProngAvgdEdxInLAr.resize(nprimaries,0.); + dir_pol_x.resize(nprimaries,0.); + dir_pol_y.resize(nprimaries,0.); + dir_pol_z.resize(nprimaries,0.); + dir_coc_x.resize(nprimaries,0.); + dir_coc_y.resize(nprimaries,0.); + dir_coc_z.resize(nprimaries,0.); + + for (auto iPrim : primaryIDs ) + { + primaryPDG[iPrim] = primaries[iPrim].PDG(); + + float_t totProngE = ProngEInLAr[iPrim]+ProngEInHadCal[iPrim]; + if (totProngE>0) + { + ShowerWidth[iPrim] = ShowerWidth[iPrim] / totProngE; } - if (ProngEInLAr[iPrim]>0) { - ShowerWidthInLAr[iPrim] = ShowerWidthInLAr[iPrim]/ProngEInLAr[iPrim]; + if (ProngEInLAr[iPrim] > 0) + { + ShowerWidthInLAr[iPrim] = ShowerWidthInLAr[iPrim] / ProngEInLAr[iPrim]; } double ShowerP = primaries[iPrim].P(); - double costheta = primaries[iPrim].Pz()/ShowerP; + double costheta = primaries[iPrim].Pz() / ShowerP; ProngAngleToBeamDir[iPrim] = TMath::ACos(costheta); - ProngAvgdEdx[iPrim] = (ProngEInLAr[iPrim] + - ProngEInHadCal[iPrim] + - ProngEInMuonFinder[iPrim])/ShowerLength[iPrim]; - ProngAvgdEdxInLAr[iPrim] = ProngEInLAr[iPrim]/ShowerLengthInLAr[iPrim]; - - std::cout<Get3DPixelMap(), - neutrino.NuVx(), neutrino.NuVy(), neutrino.NuVz(), 0., 0., 1.); - Double_t* ptr_dedx = shwlid->GetTotalDedxLongitudinal(); - std::copy(ptr_dedx, ptr_dedx+3000, TotalDedxLongitudinal); - - for(int iPrim= 0; iPrim< nPrimaryParticle; ++iPrim) { - directionfitter::LinearFit* linFit = new directionfitter::LinearFit( - pm3D->Get2DPixelMapZX(iPrim+1), - pm3D->Get2DPixelMapZY(iPrim+1), - pm3D->Get2DVtxPixelMapZX(iPrim+1), - pm3D->Get2DVtxPixelMapZY(iPrim+1), - neutrino.NuVx(), neutrino.NuVy(), neutrino.NuVz(), + slid::ShowerLID *shwlid = new slid::ShowerLID(pm3D->Get3DPixelMap(), + neutrino.NuVx(), neutrino.NuVy(), neutrino.NuVz(), 0., 0., 1.); + Double_t *ptr_dedx = shwlid->GetTotalDedxLongitudinal(); + std::copy(ptr_dedx, ptr_dedx + 3000, TotalDedxLongitudinal.begin()); + + for (int iPrim = 0; iPrim < nprimaries; ++iPrim) + { + directionfitter::LinearFit *linFit = new directionfitter::LinearFit( + pm3D->Get2DPixelMapZX(iPrim + 1), + pm3D->Get2DPixelMapZY(iPrim + 1), + pm3D->Get2DVtxPixelMapZX(iPrim + 1), + pm3D->Get2DVtxPixelMapZY(iPrim + 1), + neutrino.NuVx(), neutrino.NuVy(), neutrino.NuVz(), primaries[iPrim].Vx(), primaries[iPrim].Vy(), primaries[iPrim].Vz()); dir_pol_x[iPrim] = linFit->GetDir().X(); dir_pol_y[iPrim] = linFit->GetDir().Y(); @@ -1160,6 +979,7 @@ void AnalysisManager::FillPseudoRecoVar() { dir_coc_z[iPrim] = linFit->GetCOCDir().Z(); delete linFit; } - delete shwlid; + + fFLArEPseudoReco->Fill(); } diff --git a/src/AnalysisManagerMessenger.cc b/src/AnalysisManagerMessenger.cc index 8485a58..428eed5 100644 --- a/src/AnalysisManagerMessenger.cc +++ b/src/AnalysisManagerMessenger.cc @@ -1,4 +1,4 @@ -// +// // ******************************************************************** // * License and Disclaimer * // * * @@ -40,83 +40,89 @@ //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... AnalysisManagerMessenger::AnalysisManagerMessenger(AnalysisManager* manager) -:anamanager (manager) + :fAnalysisManager (manager) { - histoDir = new G4UIdirectory("/histo/"); - histoDir->SetGuidance("histograms control"); - - factoryCmd = new G4UIcmdWithAString("/histo/fileName", this); - factoryCmd->SetGuidance("set name for the histograms file"); - factoryCmd->AvailableForStates(G4State_PreInit,G4State_Idle); - - saveHitCmd = new G4UIcmdWithABool("/histo/saveHit", this); - saveHitCmd->SetGuidance("whether save the information of all hits"); - saveHitCmd->SetParameterName("saveHit", true); - saveHitCmd->SetDefaultValue(false); - - saveTrackCmd = new G4UIcmdWithABool("/histo/saveTrack", this); - saveTrackCmd->SetGuidance("whether save the information of all tracks"); - saveTrackCmd->SetParameterName("saveTrack", true); - saveTrackCmd->SetDefaultValue(false); - - save3DEvdCmd = new G4UIcmdWithABool("/histo/save3DEvd", this); - save3DEvdCmd->SetGuidance("whether save 3D Pixel Map"); - save3DEvdCmd->SetParameterName("save3DEvd", true); - save3DEvdCmd->SetDefaultValue(false); - - save2DEvdCmd = new G4UIcmdWithABool("/histo/save2DEvd", this); - save2DEvdCmd->SetGuidance("whether save 2D Pixel Map"); - save2DEvdCmd->SetParameterName("save2DEvd", true); - save2DEvdCmd->SetDefaultValue(false); - - circleFitCmd = new G4UIcmdWithABool("/histo/circleFit", this); - circleFitCmd->SetGuidance("perform circular fit"); - circleFitCmd->SetParameterName("circleFit", true); - circleFitCmd->SetDefaultValue(false); - - addDiffusionCmd = new G4UIcmdWithAString("/histo/addDiffusion", this); - addDiffusionCmd->SetGuidance("add toy diffusion effect"); - addDiffusionCmd->AvailableForStates(G4State_PreInit,G4State_Idle); - - saveActsCmd = new G4UIcmdWithABool("/histo/actsHits", this); - saveActsCmd->SetGuidance("save hits in Acts format"); - saveActsCmd->SetParameterName("actsHits", true); - saveActsCmd->SetDefaultValue(true); - - // histoCmd = new G4UIcmdWithAnInteger("/histo/setSeed",this); - // histoCmd->SetGuidance("Set random seed :"); - // histoCmd->SetDefaultValue(1); + fOutDir = new G4UIdirectory("/out/"); + fOutDir->SetGuidance("output control"); + + fFileCmd = new G4UIcmdWithAString("/out/fileName", this); + fFileCmd->SetGuidance("set name for the histograms file"); + fFileCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + + fSaveTrackCmd = new G4UIcmdWithABool("/out/saveTrack", this); + fSaveTrackCmd->SetGuidance("whether save the information of all tracks"); + fSaveTrackCmd->SetParameterName("saveTrack", true); + fSaveTrackCmd->SetDefaultValue(false); + + fFLArEDir = new G4UIdirectory("/out/flare/"); + fFLArEDir->SetGuidance("flare output control"); + + fEnableFLArEOutCmd = new G4UIcmdWithABool("/out/flare/enableOutput", this); + fEnableFLArEOutCmd->SetGuidance("whether to save FLArE output"); + fEnableFLArEOutCmd->SetParameterName("enableOutput", true); + fEnableFLArEOutCmd->SetDefaultValue(true); + + fSave3DEvdCmd = new G4UIcmdWithABool("/out/flare/save3DEvd", this); + fSave3DEvdCmd->SetGuidance("whether save 3D Pixel Map"); + fSave3DEvdCmd->SetParameterName("save3DEvd", true); + fSave3DEvdCmd->SetDefaultValue(false); + + fSave2DEvdCmd = new G4UIcmdWithABool("/out/flare/save2DEvd", this); + fSave2DEvdCmd->SetGuidance("whether save 2D Pixel Map"); + fSave2DEvdCmd->SetParameterName("save2DEvd", true); + fSave2DEvdCmd->SetDefaultValue(false); + + fAddDiffusionCmd = new G4UIcmdWithAString("/out/flare/addDiffusion", this); + fAddDiffusionCmd->SetGuidance("add toy diffusion effect"); + fAddDiffusionCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + + fPseudoRecoCmd = new G4UIcmdWithABool("/out/flare/pseudoReco", this); + fPseudoRecoCmd->SetGuidance("whether saving pseudo reco variables"); + fPseudoRecoCmd->SetParameterName("pseudoReco", true); + fPseudoRecoCmd->SetDefaultValue(false); + + fFASER2Dir = new G4UIdirectory("/out/faser/"); + fFASER2Dir->SetGuidance("flare output control"); + + fSaveActsCmd = new G4UIcmdWithABool("/out/faser/actsHits", this); + fSaveActsCmd->SetGuidance("save hits in Acts format"); + fSaveActsCmd->SetParameterName("actsHits", true); + fSaveActsCmd->SetDefaultValue(true); + } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... AnalysisManagerMessenger::~AnalysisManagerMessenger() { - delete factoryCmd; - delete saveHitCmd; - delete saveTrackCmd; - delete save3DEvdCmd; - delete save2DEvdCmd; - delete circleFitCmd; - delete saveActsCmd; - delete addDiffusionCmd; - delete histoDir; + delete fFileCmd; + delete fSaveTrackCmd; + delete fSave3DEvdCmd; + delete fSave2DEvdCmd; + delete fAddDiffusionCmd; + delete fPseudoRecoCmd; + delete fEnableFLArEOutCmd; + delete fSaveActsCmd; + delete fOutDir; + delete fFLArEDir; + delete fFASER2Dir; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void AnalysisManagerMessenger::SetNewValue(G4UIcommand* command,G4String newValues) { - if (command == factoryCmd) anamanager->setFileName(newValues); - if (command == saveHitCmd) anamanager->saveHit(saveHitCmd->GetNewBoolValue(newValues)); - if (command == saveTrackCmd) anamanager->saveTrack(saveTrackCmd->GetNewBoolValue(newValues)); - if (command == save3DEvdCmd) anamanager->save3DEvd(save3DEvdCmd->GetNewBoolValue(newValues)); - if (command == save2DEvdCmd) anamanager->save2DEvd(save2DEvdCmd->GetNewBoolValue(newValues)); - if (command == circleFitCmd) anamanager->circleFit(circleFitCmd->GetNewBoolValue(newValues)); - if (command == saveActsCmd) anamanager->saveActs(saveActsCmd->GetNewBoolValue(newValues)); - if (command == addDiffusionCmd) anamanager->addDiffusion(newValues); - - // if(command == histoCmd) histo->setSeed(newValues); + if (command == fFileCmd) fAnalysisManager->setFileName(newValues); + if (command == fSaveTrackCmd) fAnalysisManager->saveTrack(fSaveTrackCmd->GetNewBoolValue(newValues)); + + if (command == fEnableFLArEOutCmd) fAnalysisManager->enableFLArE(fEnableFLArEOutCmd->GetNewBoolValue(newValues)); + if (command == fSave3DEvdCmd) fAnalysisManager->save3DEvd(fSave3DEvdCmd->GetNewBoolValue(newValues)); + if (command == fSave2DEvdCmd) fAnalysisManager->save2DEvd(fSave2DEvdCmd->GetNewBoolValue(newValues)); + if (command == fPseudoRecoCmd) fAnalysisManager->savePseudoReco(fPseudoRecoCmd->GetNewBoolValue(newValues)); + if (command == fAddDiffusionCmd) fAnalysisManager->addDiffusion(newValues); + + if (command == fSaveActsCmd) fAnalysisManager->saveActs(fSaveActsCmd->GetNewBoolValue(newValues)); + } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/src/DetectorConstruction.cc b/src/DetectorConstruction.cc index db5517a..6d11af7 100644 --- a/src/DetectorConstruction.cc +++ b/src/DetectorConstruction.cc @@ -132,7 +132,6 @@ G4VPhysicalVolume* DetectorConstruction::Construct() // FLArE TPC volume if (m_addFLArE) { - GeometricalParameters::Get()->SetAddFLArE(m_addFLArE); FLArETPCDetectorConstruction *FLArETPCAssembler = new FLArETPCDetectorConstruction(); G4double lArSizeZ = GeometricalParameters::Get()->GetTPCSizeZ(); G4double TPCInsulationThickness = GeometricalParameters::Get()->GetTPCInsulationThickness(); @@ -203,7 +202,6 @@ G4VPhysicalVolume* DetectorConstruction::Construct() // FORMOSA if (m_addFORMOSA) { - GeometricalParameters::Get()->SetAddFORMOSA(m_addFORMOSA); FORMOSADetectorConstruction *FORMOSAAssembler = new FORMOSADetectorConstruction(); G4LogicalVolume* FORMOSAAssembly = FORMOSAAssembler->GetFORMOSAAssembly(); FORMOSAScintillatorBarLogical = FORMOSAAssembler->GetScintillatorBar(); @@ -222,7 +220,6 @@ G4VPhysicalVolume* DetectorConstruction::Construct() // FASERnu2 Emulsion Detector if (m_addFASERnu2) { - GeometricalParameters::Get()->SetAddFASERnu2(m_addFASERnu2); FASERnu2DetectorConstruction *FASERnu2Assembler = new FASERnu2DetectorConstruction(); FASERnu2EmulsionLogical = FASERnu2Assembler->GetEmulsionFilm(); FASERnu2VetoInterfaceLogical = FASERnu2Assembler->GetVetoInterfaceDetector(); @@ -242,7 +239,6 @@ G4VPhysicalVolume* DetectorConstruction::Construct() // FASER2 Magnet + Tracking stations if (m_addFASER2) { - GeometricalParameters::Get()->SetAddFASER2(m_addFASER2); FASER2DetectorConstruction *magnetAssembler = new FASER2DetectorConstruction(); FASER2MagnetLogical = magnetAssembler->GetMagneticVolume(); //need to assign B field FASER2TrackingLogical = magnetAssembler->GetTrackingStations(); @@ -290,24 +286,24 @@ void DetectorConstruction::ConstructSDandField() { if (m_addFLArE) { - LArBoxSD* TPCModuleSD = new LArBoxSD("lArBoxSD"); + LArBoxSD* TPCModuleSD = new LArBoxSD("FLArEBoxSD"); TPCModuleLogical->SetSensitiveDetector(TPCModuleSD); sdManager->AddNewDetector(TPCModuleSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "lArBoxSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx,"FLArEBoxSD/lar_box"); SDIdx++; if (m_useBabyMIND) { - LArBoxSD* BabyMINDHorBarSD = new LArBoxSD("BabyMINDHorBarSD"); + LArBoxSD* BabyMINDHorBarSD = new LArBoxSD("FLArEBabyMINDHorBarSD"); BabyMINDHorizontalBar->SetSensitiveDetector(BabyMINDHorBarSD); sdManager->AddNewDetector(BabyMINDHorBarSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "BabyMINDHorBarSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx,"FLArEBabyMINDHorBarSD/lar_box"); SDIdx++; - LArBoxSD* BabyMINDVerBarSD = new LArBoxSD("BabyMINDVerBarSD"); + LArBoxSD* BabyMINDVerBarSD = new LArBoxSD("FLArEBabyMINDVerBarSD"); BabyMINDVerticalBar->SetSensitiveDetector(BabyMINDVerBarSD); sdManager->AddNewDetector(BabyMINDVerBarSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "BabyMINDVerBarSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx,"FLArEBabyMINDVerBarSD/lar_box"); SDIdx++; // magnetic field for BabyMIND @@ -319,40 +315,40 @@ void DetectorConstruction::ConstructSDandField() { } else { - LArBoxSD* HadCalXSD = new LArBoxSD("HadCalXSD"); + LArBoxSD* HadCalXSD = new LArBoxSD("FLArEHadCalXSD"); HadCalXCellLogical->SetSensitiveDetector(HadCalXSD); sdManager->AddNewDetector(HadCalXSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "HadCalXSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx, "FLArEHadCalXSD/lar_box"); SDIdx++; - LArBoxSD* HadCalYSD = new LArBoxSD("HadCalYSD"); + LArBoxSD* HadCalYSD = new LArBoxSD("FLArEHadCalYSD"); HadCalYCellLogical->SetSensitiveDetector(HadCalYSD); sdManager->AddNewDetector(HadCalYSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "HadCalYSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx, "FLArEHadCalYSD/lar_box"); SDIdx++; - LArBoxSD* MuonFinderXSD = new LArBoxSD("MuonFinderXSD"); + LArBoxSD* MuonFinderXSD = new LArBoxSD("FLArEMuonFinderXSD"); MuonFinderXCellLogical->SetSensitiveDetector(MuonFinderXSD); sdManager->AddNewDetector(MuonFinderXSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "MuonFinderXSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx, "FLArEMuonFinderXSD/lar_box"); SDIdx++; - LArBoxSD* MuonFinderYSD = new LArBoxSD("MuonFinderYSD"); + LArBoxSD* MuonFinderYSD = new LArBoxSD("FLArEMuonFinderYSD"); MuonFinderYCellLogical->SetSensitiveDetector(MuonFinderYSD); sdManager->AddNewDetector(MuonFinderYSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "MuonFinderYSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx, "FLArEMuonFinderYSD/lar_box"); SDIdx++; - LArBoxSD* HadAbsorbSD = new LArBoxSD("HadAbsorbSD"); + LArBoxSD* HadAbsorbSD = new LArBoxSD("FLArEHadAbsorbSD"); HadAbsorLayersLogical->SetSensitiveDetector(HadAbsorbSD); sdManager->AddNewDetector(HadAbsorbSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "HadAbsorbSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx, "FLArEHadAbsorbSD/lar_box"); SDIdx++; - LArBoxSD* MuonFinderAbsorbSD = new LArBoxSD("MuonFinderAbsorbSD"); + LArBoxSD* MuonFinderAbsorbSD = new LArBoxSD("FLArEMuonFinderAbsorbSD"); MuonFinderAbsorLayersLogical->SetSensitiveDetector(MuonFinderAbsorbSD); sdManager->AddNewDetector(MuonFinderAbsorbSD); - GeometricalParameters::Get()->AddSD2List(SDIdx, "MuonFinderAbsorbSD/lar_box"); + GeometricalParameters::Get()->AddSD2List(SDIdx, "FLArEMuonFinderAbsorbSD/lar_box"); SDIdx++; // magnetic field for HadCatcher + MuonFinder diff --git a/src/EventAction.cc b/src/EventAction.cc index 0f0ebb5..97e2d43 100644 --- a/src/EventAction.cc +++ b/src/EventAction.cc @@ -11,17 +11,19 @@ EventAction::EventAction() : G4UserEventAction(), fNPrimaryTrack("NPrimaryTrack", 0), fNSecondaryTrack("NSecondaryTrack", 0), - fNSecondaryTrackNotGamma("NSecondaryTrackNotGamma", 0) { - // Register created accumulables - G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance(); - accumulableManager->RegisterAccumulable(fNPrimaryTrack); - accumulableManager->RegisterAccumulable(fNSecondaryTrack); - accumulableManager->RegisterAccumulable(fNSecondaryTrackNotGamma); - } + fNSecondaryTrackNotGamma("NSecondaryTrackNotGamma", 0) +{ + // Register created accumulables + G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance(); + accumulableManager->RegisterAccumulable(fNPrimaryTrack); + accumulableManager->RegisterAccumulable(fNSecondaryTrack); + accumulableManager->RegisterAccumulable(fNSecondaryTrackNotGamma); +} EventAction::~EventAction() {;} -void EventAction::BeginOfEventAction(const G4Event* event) { +void EventAction::BeginOfEventAction(const G4Event* event) +{ // Reset all accumulables to their initial values G4AccumulableManager* accumulableManager = G4AccumulableManager::Instance(); accumulableManager->Reset(); @@ -30,37 +32,45 @@ void EventAction::BeginOfEventAction(const G4Event* event) { ana->BeginOfEvent(); } -void EventAction::EndOfEventAction(const G4Event* event) { - G4cout<<"This is the "<GetEventID()<<"th event"<GetEventID() << "th event"<EndOfEvent(event); + } -void EventAction::AddPrimaryTrack() { +void EventAction::AddPrimaryTrack() +{ fNPrimaryTrack += 1; } -void EventAction::AddSecondaryTrack() { +void EventAction::AddSecondaryTrack() +{ fNSecondaryTrack += 1; } -void EventAction::AddSecondaryTrackNotGamma() { +void EventAction::AddSecondaryTrackNotGamma() +{ fNSecondaryTrackNotGamma += 1; } diff --git a/src/FPFParticle.cc b/src/FPFParticle.cc index 1c05d2f..f61890a 100644 --- a/src/FPFParticle.cc +++ b/src/FPFParticle.cc @@ -18,7 +18,7 @@ FPFParticle::FPFParticle() , fPx(-999) , fPy(-999) , fPz(-999) - , fE(-999) + , fE(-999) { } @@ -49,7 +49,7 @@ FPFParticle::FPFParticle(const int pdg, , fPx(Px) , fPy(Py) , fPz(Pz) - , fE(E) + , fE(E) { } diff --git a/src/LArBoxSD.cc b/src/LArBoxSD.cc index 77bb537..cdc88d4 100644 --- a/src/LArBoxSD.cc +++ b/src/LArBoxSD.cc @@ -17,6 +17,12 @@ LArBoxSD::LArBoxSD(G4String name) : G4VSensitiveDetector(name) { G4bool LArBoxSD::ProcessHits(G4Step* aStep, G4TouchableHistory* R0hist) { G4Track* aTrack = aStep->GetTrack(); + G4float edep = aStep->GetTotalEnergyDeposit(); + if(edep <= 0.02*MeV) //no energy deposited, skip + { + return false; + } + // enumerator: fAlive, fStopButAlive, fStopAndKill, fKillTrackAndSecondaries, fSuspend, fPostponeToNextEvent // https://apc.u-paris.fr/~franco/g4doxy/html/G4TrackStatus_8hh.html#734825af9cdc612606614fdce0545157 // http://geant4.in2p3.fr/2005/Workshop/ShortCourse/session1/M.Asai.pdf @@ -24,10 +30,10 @@ G4bool LArBoxSD::ProcessHits(G4Step* aStep, G4TouchableHistory* R0hist) { //G4cout<<"debug (track status): "<GetVertexPosition(); - G4double TrackLength = aTrack->GetTrackLength(); + G4float TrackLength = aTrack->GetTrackLength(); G4int ParticleName = aTrack->GetParticleDefinition()->GetPDGEncoding(); - G4double ParticleMass = aTrack->GetParticleDefinition()->GetPDGMass(); + G4float ParticleMass = aTrack->GetParticleDefinition()->GetPDGMass(); G4int PID = aTrack->GetParentID(); G4int TID = aTrack->GetTrackID(); G4int Stepno = aTrack->GetCurrentStepNumber(); @@ -60,15 +66,18 @@ G4bool LArBoxSD::ProcessHits(G4Step* aStep, G4TouchableHistory* R0hist) { // } // extra info - G4double edep = aStep->GetTotalEnergyDeposit(); G4ThreeVector pos = 0.5*(PreStepPosition + PostStepPosition); G4String PosVolume = PostStep->GetPhysicalVolume()->GetName(); G4int StepStatus = PostStep->GetStepStatus(); + G4int copyNum = aStep->GetPreStepPoint()->GetTouchableHandle()->GetCopyNumber(); + G4ThreeVector delta_momentum = (PostStep->GetMomentum())-InitMomentum; TrackInformation* aTrackInfo = (TrackInformation*)(aTrack->GetUserInformation()); G4int trackIsFromPrimaryPizero = 0; G4int trackIsFromFSLPizero = 0; G4int trackIsFromPrimaryLepton = 0; + G4int trackPDG = aTrack->GetParticleDefinition()->GetPDGEncoding(); + G4float Time = aTrack->GetLocalTime(); if (aTrackInfo) { trackIsFromPrimaryPizero = aTrackInfo->IsTrackFromPrimaryPizero(); trackIsFromFSLPizero = aTrackInfo->IsTrackFromFSLPizero(); @@ -85,6 +94,12 @@ G4bool LArBoxSD::ProcessHits(G4Step* aStep, G4TouchableHistory* R0hist) { hit->SetPID(PID); hit->SetTID(TID); hit->SetStepNo(Stepno); + + hit->SetCopyNum(copyNum); + hit->SetPDG(trackPDG); + hit->SetDeltaMom(delta_momentum); + hit->SetTime(Time); + hit->SetPreStepPosition(PreStepPosition); hit->SetPostStepPosition(PostStepPosition); hit->SetInitMomentum(InitMomentum); diff --git a/src/PixelMap3D.cc b/src/PixelMap3D.cc index 3eee24d..5e6feb6 100644 --- a/src/PixelMap3D.cc +++ b/src/PixelMap3D.cc @@ -17,7 +17,7 @@ using namespace hep_hpc::hdf5; PixelMap3D::PixelMap3D(const Int_t evtID, const Int_t nPrim, const Int_t PDG, const Double_t* res) : fEvtID(evtID), fNPrim(nPrim), fGeneratorPDG(PDG) { - G4cout<<"Creating PixelMap3D for event "<cd(thedir->GetName()); std::string dirname = "edep2D/evt_"+std::to_string(fEvtID)+"/"; - thefile->mkdir(dirname.c_str()); - thefile->cd(dirname.c_str()); + thedir->mkdir(dirname.c_str()); + thedir->cd(dirname.c_str()); for (int i=0; iGetEntries()==0) hitClusterZX[i]->SetEntries(1); if (hitClusterZY[i]->GetEntries()==0) hitClusterZY[i]->SetEntries(1); hitClusterZX[i]->Write(); hitClusterZY[i]->Write(); } + thefile->cd(thedir->GetName()); dirname = "edep2Dvtx/evt_"+std::to_string(fEvtID)+"/"; - thefile->mkdir(dirname.c_str()); - thefile->cd(dirname.c_str()); + thedir->mkdir(dirname.c_str()); + thedir->cd(dirname.c_str()); for (int i=0; iGetEntries()==0) vtxHitClusterZX[i]->SetEntries(1); if (vtxHitClusterZY[i]->GetEntries()==0) vtxHitClusterZY[i]->SetEntries(1); vtxHitClusterZX[i]->Write(); vtxHitClusterZY[i]->Write(); } + thefile->cd(); } G4double PixelMap3D::DistanceToAnode(G4double x) { diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc index 8701f4d..cd659b0 100644 --- a/src/PrimaryGeneratorAction.cc +++ b/src/PrimaryGeneratorAction.cc @@ -15,18 +15,18 @@ #include "G4Exception.hh" -PrimaryGeneratorAction::PrimaryGeneratorAction() -{ +PrimaryGeneratorAction::PrimaryGeneratorAction() +{ // create a messenger for this class fGenMessenger = new PrimaryGeneratorMessenger(this); // start with default generator - fGenerator = new GPSGenerator(); + fGenerator = new GPSGenerator(); fInitialized = false; } -PrimaryGeneratorAction::~PrimaryGeneratorAction() +PrimaryGeneratorAction::~PrimaryGeneratorAction() { delete fGenerator; delete fGenMessenger; @@ -50,14 +50,14 @@ void PrimaryGeneratorAction::SetGenerator(G4String name) "UnknownOption", FatalErrorInArgument, err.c_str()); - } + } } void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) -{ +{ // load generator data at first event - // this function opens files, reads trees, etc (if required) + // this function opens files, reads trees, etc (if required) if(!fInitialized){ fGenerator->LoadData(); fInitialized = true; @@ -71,8 +71,8 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) // save additional truth information alongside primary particles // this makes it available for the output tree (e.g: neutrino info for genie) - - // TODO: make it more generic for other generators? + + // TODO: make it more generic for other generators? bool isGenie = (fGenerator->GetGeneratorName() == "genie"); // downcast: if it fails, we won't use it anyway... GENIEGenerator *genieGen = dynamic_cast(fGenerator); @@ -91,12 +91,13 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) G4int count_particles = 0; for (G4int ivtx = 0; ivtx < anEvent->GetNumberOfPrimaryVertex(); ++ivtx) { for (G4int ipp = 0; ipp < anEvent->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); ++ipp) { - + G4PrimaryParticle* primary_particle = anEvent->GetPrimaryVertex(ivtx)->GetPrimary(ipp); - + if (primary_particle) { primary_particle->SetUserInformation(new PrimaryParticleInformation( - count_particles, primary_particle->GetPDGcode(), primary_particle->GetMass(), primary_particle->GetCharge(), + count_particles, primary_particle->GetPDGcode(), primary_particle->GetMass(), + primary_particle->GetCharge(), primary_particle->GetMomentum(), anEvent->GetPrimaryVertex(ivtx)->GetPosition(), neuidx, neupdg, neup4, neux4, int_type, scattering_type, w, fslpdg, fslp4)); diff --git a/src/PrimaryParticleInformation.cc b/src/PrimaryParticleInformation.cc index 872c0b6..2413429 100644 --- a/src/PrimaryParticleInformation.cc +++ b/src/PrimaryParticleInformation.cc @@ -1,12 +1,12 @@ #include "PrimaryParticleInformation.hh" PrimaryParticleInformation::PrimaryParticleInformation(G4int aID, G4int aPDG, G4double aMass, G4double aCharge, - G4ThreeVector aMomentum, G4ThreeVector aVertex, + G4ThreeVector aMomentum, G4ThreeVector aVertex, G4int aneuIdx, G4int aneuPDG, TLorentzVector aneuP4, TLorentzVector aneuX4, G4int aInttype, G4int aScatteringtype, G4double aW, G4int afslPDG, TLorentzVector afslP4 ) : fPartID(aID), fPDG(aPDG), fMass(aMass), fCharge{aCharge}, - fMomentumMC(aMomentum), fVertexMC(aVertex), + fMomentumMC(aMomentum), fVertexMC(aVertex), fNeuIdx(aneuIdx), fNeuPDG(aneuPDG), fNeuP4(aneuP4), fNeuX4(aneuX4), fInteractionType(aInttype), fScatteringType(aScatteringtype), fW(aW), fFSLPDG(afslPDG), fFSLP4(afslP4) diff --git a/src/StackingAction.cc b/src/StackingAction.cc index 8aaa3af..ba00128 100644 --- a/src/StackingAction.cc +++ b/src/StackingAction.cc @@ -10,23 +10,37 @@ StackingAction::StackingAction(RunAction* aRunAction, EventAction* aEventAction) G4UserStackingAction(), fRunAction(aRunAction), fEventAction(aEventAction) {;} -G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack (const G4Track* aTrack) { - - AnalysisManager::GetInstance()->SetTrackPTPair(aTrack->GetParentID(), aTrack->GetTrackID()); +G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack (const G4Track* aTrack) +{ + // for each track, build track ID to primary ancestor + // primaries have themselves as ancestor + // go up the tree until original primary for everything else + G4int trackID = aTrack->GetTrackID(); + G4int parentID = aTrack->GetParentID(); + G4int ancestorID = -1; // Register primary tracks - if (aTrack->GetParentID()==0) { + if (parentID==0) + { fEventAction->AddPrimaryTrack(); + ancestorID = trackID; // primary is its own ancestor } // Register only secondaries, i.e. tracks having ParentID > 0 - if (aTrack->GetParentID()) { + if (parentID > 0) + { fEventAction->AddSecondaryTrack(); - if (aTrack->GetParticleDefinition()->GetPDGEncoding() != 22) { + if (aTrack->GetParticleDefinition()->GetPDGEncoding() != 22) + { fEventAction->AddSecondaryTrackNotGamma(); } + + // if you are a secondary, ancestor is your parent's ancestor + ancestorID = AnalysisManager::GetInstance()->GetTrackPrimaryAncestor(parentID); } + // add track with its ancestor!!! + AnalysisManager::GetInstance()->SetTrackPrimaryAncestor(trackID,ancestorID); // Do not affect track classification. Just return what would have // been returned by the base class diff --git a/src/SteppingAction.cc b/src/SteppingAction.cc index 9aa70cb..2bbc4f0 100644 --- a/src/SteppingAction.cc +++ b/src/SteppingAction.cc @@ -15,30 +15,30 @@ SteppingAction::SteppingAction(RunAction* runAction) } void SteppingAction::UserSteppingAction(const G4Step* aStep) { - + //TrackLiveDebugging(aStep); - + G4Track* aTrack = aStep->GetTrack(); G4ThreeVector post_pos = aStep->GetPostStepPoint()->GetPosition(); - + // if the track is out of the active volumes/area, kill this track G4VPhysicalVolume* volume = aStep->GetPostStepPoint()->GetTouchable()->GetVolume(); - + if( volume->GetName() == "worldPV" ) aTrack->SetTrackStatus(G4TrackStatus::fStopAndKill); - //else if( volume->GetName() == "hallPV" ){ + //else if( volume->GetName() == "hallPV" ){ // G4double endFLArE = GeometricalParameters::Get()->GetFLArEPosition()+ // GeometricalParameters::Get()->GetTPCSizeZ()/2.; - // G4double totYokeX = 4.5 * m; + // G4double totYokeX = 4.5 * m; // G4double totYokeY = 3 *m; // // //if before end of FLArE, kill - // if( post_pos.z() < endFLArE ) aTrack->SetTrackStatus(G4TrackStatus::fStopAndKill); - // + // if( post_pos.z() < endFLArE ) aTrack->SetTrackStatus(G4TrackStatus::fStopAndKill); + // // //if after enf of FLArE, kill if outside FASER magnet cross-section // else if ( TMath::Abs(post_pos.x()) > totYokeX/2. || TMath::Abs(post_pos.y()) > totYokeY/2. ) // aTrack->SetTrackStatus(G4TrackStatus::fStopAndKill); //} - + /*// if the track is out of the active volumes, kill this track G4VPhysicalVolume* volume = aStep->GetPostStepPoint()->GetTouchable()->GetVolume(); std::string active_volumes[9] = {"lArBox", @@ -63,11 +63,11 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) { // for tracks in HadCatcher and MuonFinder absorber if (std::find(std::begin(active_logical), std::end(active_logical), mother_logical->GetName()) ==std::end(active_logical)) { - aTrack->SetTrackStatus(G4TrackStatus::fStopAndKill); + aTrack->SetTrackStatus(G4TrackStatus::fStopAndKill); } } else { // for tracks in world - aTrack->SetTrackStatus(G4TrackStatus::fStopAndKill); + aTrack->SetTrackStatus(G4TrackStatus::fStopAndKill); } } }*/ @@ -87,12 +87,12 @@ void SteppingAction::TrackLiveDebugging(const G4Step* step){ G4ThreeVector pre_mom = step->GetPreStepPoint()->GetMomentum(); G4ThreeVector post_mom = step->GetPostStepPoint()->GetMomentum(); G4double edep = step->GetTotalEnergyDeposit(); - + int PDG = track->GetParticleDefinition()->GetPDGEncoding(); G4String ParticleName = track->GetDynamicParticle()->GetParticleDefinition()->GetParticleName(); int TID = track->GetTrackID(); int SID = step->GetTrack()->GetCurrentStepNumber(); - + std::cout << "Track " << TID << " - " << "PDG " << PDG << " " << ParticleName << std::endl; std::cout << "stepping... " << SID << " edep" << edep << std::endl; std::cout << "(" << pre_pos.x() << "," << pre_pos.y() << "," << pre_pos.z() << ") in " << volume->GetName(); diff --git a/src/TrackingAction.cc b/src/TrackingAction.cc index dd9113e..900d7fa 100644 --- a/src/TrackingAction.cc +++ b/src/TrackingAction.cc @@ -13,19 +13,26 @@ void TrackingAction::PreUserTrackingAction(const G4Track* aTrack) void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) { - if (aTrack->GetParentID()==0) { + if (aTrack->GetParentID()==0) + { AnalysisManager::GetInstance()->AddOnePrimaryTrack(); } - if (aTrack->GetParentID()==0) { - if (aTrack->GetParticleDefinition()->GetPDGEncoding()==111) { + if (aTrack->GetParentID()==0) + { + if (aTrack->GetParticleDefinition()->GetPDGEncoding()==111) + { // in case of pizero in the list of primary track // its decay products are also counted as primary particles, mostly 2 gammas G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); - if (secondaries) { + if (secondaries) + { size_t nSeco = secondaries->size(); - if (nSeco>0) { - for (size_t i=0; iGetCreatorProcess()->GetProcessName()=="Decay") { + if (nSeco>0) + { + for (size_t i=0; iGetCreatorProcess()->GetProcessName()=="Decay") + { TrackInformation* info = new TrackInformation(); info->SetTrackIsFromPrimaryPizero(1); (*secondaries)[i]->SetUserInformation(info); @@ -33,21 +40,27 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) } } } - } + } } } + if (aTrack->GetTrackID()==1 && (abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==15 || - abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==13)) { + abs(aTrack->GetParticleDefinition()->GetPDGEncoding())==13)) + { // in case of the lepton decays, the decay products are counted as primary particles // * tau- decay (dominant) // * mu- decay G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); - if (secondaries) { + if (secondaries) + { size_t nSeco = secondaries->size(); - if (nSeco>0) { - for (size_t i=0; iGetCreatorProcess()->GetProcessName()=="Decay") { + if (nSeco>0) + { + for (size_t i=0; iGetCreatorProcess()->GetProcessName()=="Decay") + { TrackInformation* info = new TrackInformation(); info->SetTrackIsFromPrimaryLepton(1); (*secondaries)[i]->SetUserInformation(info); @@ -57,16 +70,23 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) } } } - if (aTrack->GetParentID()==1 && aTrack->GetCreatorProcess()->GetProcessName()=="Decay") { + + if (aTrack->GetParentID()==1 && aTrack->GetCreatorProcess()->GetProcessName()=="Decay") + { // in case of tau decay pizero // decay products of this pizero are also counted as primary particles, mostly 2 gammas - if (aTrack->GetParticleDefinition()->GetPDGEncoding()==111) { + if (aTrack->GetParticleDefinition()->GetPDGEncoding()==111) + { G4TrackVector* secondaries = fpTrackingManager->GimmeSecondaries(); - if (secondaries) { + if (secondaries) + { size_t nSeco = secondaries->size(); - if (nSeco>0) { - for (size_t i=0; iGetCreatorProcess()->GetProcessName()=="Decay") { + if (nSeco>0) + { + for (size_t i=0; iGetCreatorProcess()->GetProcessName()=="Decay") + { TrackInformation* info = new TrackInformation(); info->SetTrackIsFromFSLPizero(1); (*secondaries)[i]->SetUserInformation(info); @@ -77,6 +97,7 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) } } } + //TrackInformation* aTrackInfo = (TrackInformation*)(aTrack->GetUserInformation()); //if (aTrackInfo) { // if (aTrackInfo->IsTrackFromPrimaryTau() | aTrackInfo->IsTrackFromPrimaryPizero()) { diff --git a/src/geometry/FLArETPCDetectorConstruction.cc b/src/geometry/FLArETPCDetectorConstruction.cc index 5de38a5..bebf3c2 100644 --- a/src/geometry/FLArETPCDetectorConstruction.cc +++ b/src/geometry/FLArETPCDetectorConstruction.cc @@ -9,7 +9,7 @@ #include "G4SystemOfUnits.hh" #include "G4ThreeVector.hh" #include "G4VisAttributes.hh" -#include "G4Colour.hh" +#include "G4Colour.hh" #include "G4PVReplica.hh" #include "G4UserLimits.hh" #include "G4PVPlacement.hh" @@ -28,9 +28,9 @@ FLArETPCDetectorConstruction::FLArETPCDetectorConstruction() } else if (fDetMaterialName == GeometricalParameters::tpcMaterialOption::LiquidKrypton) { detectorMaterial = fMaterials->Material("LiquidKrypton"); G4cout<<"**** FLArE TPC Material : Liquid Krypton ****"<SetUserLimits(new G4UserLimits(0.5*mm)); } else if (fDetGeomOption == GeometricalParameters::tpcConfigOption::ThreeBySeven) { G4cout << "TPC module configuration: 3x7" << G4endl; - lArBoxLog = new G4LogicalVolume(lArBox, detectorMaterial, "TPCModuleLogical"); + lArBoxLog = new G4LogicalVolume(lArBox, detectorMaterial, "lArLogical"); auto lArBoxVis = new G4VisAttributes(G4Colour(86./255, 152./255, 195./255)); lArBoxVis->SetVisibility(false); lArBoxLog->SetVisAttributes(lArBoxVis); - G4double TPCLayerWidth = fLArSizeX; - G4double TPCLayerHeight = fLArSizeY; - G4double TPCLayerLength = fLArSizeZ / 7.0; - G4double TPCModuleWidth = TPCLayerWidth / 3.0; - G4double TPCModuleHeight = TPCLayerHeight; - G4double TPCModuleLength = TPCLayerLength; - auto TPCLayerSolid - = new G4Box("TPCLayerBox", TPCLayerWidth/2, TPCLayerHeight/2, TPCLayerLength/2); - auto TPCLayerLogical - = new G4LogicalVolume(TPCLayerSolid, detectorMaterial, "TPCLayerLogical"); - auto TPCModuleSolid - = new G4Box("TPCModuleBox", TPCModuleWidth/2, TPCModuleHeight/2, TPCModuleLength/2); - fFLArETPCLog = new G4LogicalVolume(TPCModuleSolid, detectorMaterial, "TPCModuleLog"); - new G4PVReplica("TPCModulePhysical", fFLArETPCLog, TPCLayerLogical, kXAxis, 3, TPCModuleWidth); - new G4PVReplica("TPC", TPCLayerLogical, lArBoxLog, kZAxis, 7, TPCLayerLength); + G4double TPCModuleWidth = fLArSizeX / 3.0; + G4double TPCModuleHeight = fLArSizeY; + G4double TPCModuleLength = fLArSizeZ / 7.0; + + auto TPCModuleSolid = new G4Box("TPCModuleBox", TPCModuleWidth/2, TPCModuleHeight/2, TPCModuleLength/2); + fFLArETPCLog = new G4LogicalVolume(TPCModuleSolid, detectorMaterial, "TPCModuleLogical"); + + //making indiv boxes below to replace the replicas (more copy nums saved) + int boxCt = 0; + double startEndX = (fLArSizeX/2) - 0.5*TPCModuleWidth; + double startEndZ = (fLArSizeZ/2) - 0.5*TPCModuleLength; + + for (int boxNumZ = 0; boxNumZ < 7; boxNumZ++) { + for (int boxNumX = 0; boxNumX < 3; boxNumX++){ + G4ThreeVector pos (boxNumX*TPCModuleWidth - startEndX, 0, boxNumZ*TPCModuleLength - startEndZ ); + new G4PVPlacement(0, pos, fFLArETPCLog, "TPCModulePhysical", lArBoxLog, false, boxCt); + boxCt++; + } + } + G4VisAttributes* TPCModuleVis = new G4VisAttributes(G4Colour(86./255, 152./255, 195./255)); TPCModuleVis->SetVisibility(true); TPCModuleVis->SetForceWireframe(true); TPCModuleVis->SetForceAuxEdgeVisible(true); - TPCLayerLogical->SetVisAttributes(lArBoxVis); fFLArETPCLog->SetVisAttributes(TPCModuleVis); - fFLArETPCLog->SetUserLimits(new G4UserLimits(0.5*mm)); + //fFLArETPCLog->SetUserLimits(new G4UserLimits(1.0*mm)); // FIXME? } } @@ -122,8 +126,8 @@ void FLArETPCDetectorConstruction::BuildCryostatInsulation() //----------------------------------- // insulation auto lArBox = new G4Box("lArBox", fLArSizeX/2., fLArSizeY/2., fLArSizeZ/2.); - auto CryoInsulationBlockSolid = new G4Box("CryoInsulationBlock", - fLArSizeX/2.+fThicknessInsulation, + auto CryoInsulationBlockSolid = new G4Box("CryoInsulationBlock", + fLArSizeX/2.+fThicknessInsulation, fLArSizeY/2.+fThicknessInsulation, fLArSizeZ/2.+fThicknessInsulation); auto CryoInsulationSolid = new G4SubtractionSolid("CryoInsulation", @@ -131,9 +135,9 @@ void FLArETPCDetectorConstruction::BuildCryostatInsulation() lArBox, 0, //no rotation G4ThreeVector(0,0,0)); // no translation - cryoInsulationLog = new G4LogicalVolume(CryoInsulationSolid, - fMaterials->Material("R_PUF"), - "CryoInsulationLogical"); + cryoInsulationLog = new G4LogicalVolume(CryoInsulationSolid, + fMaterials->Material("R_PUF"), + "CryoInsulationLogical"); auto CryoInsulationVis = new G4VisAttributes(G4Colour(86./255, 152./255, 195./255)); CryoInsulationVis->SetVisibility(true); CryoInsulationVis->SetForceWireframe(true); diff --git a/src/geometry/GeometricalParameters.cc b/src/geometry/GeometricalParameters.cc index cfbe1e9..d63fa22 100644 --- a/src/geometry/GeometricalParameters.cc +++ b/src/geometry/GeometricalParameters.cc @@ -12,9 +12,6 @@ GeometricalParameters::GeometricalParameters() fHallOffsetX = 1.44*m; fHallOffsetY = 2.21*m; - // FLArE TPC volume - fAddFLArE = false; //default - // rock envelope fEnableRockEnvelope = false; fRockFrontThickness = 10*m; @@ -34,7 +31,6 @@ GeometricalParameters::GeometricalParameters() fFLArEPos = G4ThreeVector(0., 0., 4300.*mm); // BabyMIND - fUseBabyMIND = false; //default fBabyMINDMagnetPlateThickness = 30*mm; fBabyMINDMagnetPlateSizeX = 3.5*m; fBabyMINDMagnetPlateSizeY = 2.*m; @@ -57,7 +53,6 @@ GeometricalParameters::GeometricalParameters() fBabyMINDBlockSequence = "|MMMMD||DMMMD||DMMMMD||MMDMMD||MMDMMD||MDMDMD||DMMMD|"; // FASER2 magnet - fAddFASER2 = false; //default fFASER2MagnetOption = magnetOption::SAMURAI; fFASER2MagnetField = 1.0*tesla; fMagnetTotalSizeZ = 4*m; //updates during construction @@ -96,7 +91,6 @@ GeometricalParameters::GeometricalParameters() fFASER2VetoShieldThickness = 10*cm; //! These values are placeholders and require dedicated optmisation studies! // FASERnu2 emulsion detector - fAddFASERnu2 = false; //default fFASERnu2TotalSizeZ = 8.5*m; //updates during construction fNEmulsionTungstenLayers = 3300; fTungstenThickness = 2*mm; @@ -109,7 +103,6 @@ GeometricalParameters::GeometricalParameters() fFASERnu2Pos = G4ThreeVector(0., 0., 22023.*mm); // FORMOSA - fAddFORMOSA = false; //default fFORMOSATotalSizeZ = 5*m; //updates during construction fNScinBarsX = 20; fNScinBarsY = 20; diff --git a/src/reco/CircleFit.cc b/src/reco/CircleFit.cc deleted file mode 100644 index e26fb5c..0000000 --- a/src/reco/CircleFit.cc +++ /dev/null @@ -1,331 +0,0 @@ -//////////////////////////////////////////////////////////////////////// -// \file CircularFit.cc -// \brief Classes to perform circular fits using analytical methods -// for momentum estimation in magnetic fields -// \author M. Vicenzi (mvicenzi@bnl.gov) -//////////////////////////////////////////////////////////////////////// - -#include "reco/CircleFit.hh" -#include "geometry/GeometricalParameters.hh" -#include "G4SystemOfUnits.hh" - -#include -#include -#include - -#include -#include -#include - -namespace circularfitter { - CircleFit::CircleFit(const std::vector x, const std::vector z) - { - fXc = -9999; - fZc = -9999; - fR = -9999; - fChi2 = -9999; - - // This code follows an analytical method for circle fitting (Modified Least Squares method) - // Ref: https://doi.org/10.1109/TIM.2003.820472 - // using paper formulas with y -> z - - if( x.size() != z.size() ) { fStatus = 1; return; } - int n = x.size(); - - double sumx = 0, sumz = 0; // linear terms - double sumx2 = 0, sumz2 = 0, sumxz = 0; // quadratic terms - double sumxz2 = 0, sumx2z = 0, sumx3 = 0, sumz3 = 0; // cubic terms - - for (int i = 0; i < n; i++) { - double xp = x.at(i); - double zp = z.at(i); - sumx += xp; - sumz += zp; - sumx2 += xp * xp; - sumz2 += zp * zp; - sumxz += xp * zp; - sumxz2 += xp * zp * zp; - sumx2z += xp * xp * zp; - sumx3 += xp * xp * xp; - sumz3 += zp * zp * zp; - } - - double a = n * sumx2 - sumx * sumx; - double b = n * sumxz - sumx * sumz; - double c = n * sumz2 - sumz * sumz; - double d = 0.5 * (n * sumxz2 - sumx * sumz2 + n * sumx3 - sumx * sumx2); - double e = 0.5 * (n * sumx2z - sumz * sumx2 + n * sumz3 - sumz * sumz2); - - if (a * c - b * b == 0.) { fStatus = 2; return; } - - //center of the circle - fXc = (d * c - b * e) / (a * c - b * b); - fZc = (a * e - b * d) / (a * c - b * b); - - //radius of the circle - fR = 0; - for (int i = 0; i < n; i++) { - double xp = x.at(i); - double zp = z.at(i); - double r2 = (xp - fXc) * (xp - fXc) + (zp - fZc) * (zp - fZc); - - fR += TMath::Sqrt(r2); - } - - fR /= n; - - //compute chi2 - fChi2 = 0.0; - for (int i = 0; i < n; i++) { - double dist = TMath::Sqrt((z.at(i) - fZc) * (z.at(i) - fZc) + (x.at(i) - fXc) * (x.at(i) - fXc) ) - fR; - fChi2 += dist*dist; - } - - fChi2 /= n; - fStatus = 0; - } - - CircleFit::~CircleFit() - { - } - -// ----------------------------------------------------------------------------------------------- -// ----------------------------------------------------------------------------------------------- - - LineFit::LineFit(const std::vector z, const std::vector y) - { - fP0 = -9999; - fP1 = -9999; - fCosDip = -9999; - - if( z.size() != y.size() ) { fStatus = 1; return; } - int n = z.size(); - - double S1 = 0.; - double SZ = 0.; - double SY = 0.; - double SYZ = 0.; - double SZZ = 0.; - - for (int i = 0; i < n; i++) { - S1 += 1.; - SZ += z.at(i); - SY += y.at(i); - SYZ += z.at(i) * y.at(i); - SZZ += z.at(i) * z.at(i); - } - - double D = S1 * SZZ - SZ * SZ; - - if (D == 0.) { fStatus = 2; return; } - - fP0 = (SY * SZZ - SZ * SYZ) / D; // i.p. at x = 0 - fP1 = (S1 * SYZ - SZ * SY) / D; // tg(theta) - fCosDip = 1./TMath::Sqrt(1+fP1*fP1); - - fChi2 = 0.; - for (int i = 0; i < n; i++) { - fChi2 += (y.at(i) - fP0 - fP1 * z.at(i)) * (y.at(i) - fP0 - fP1 * z.at(i)); - } - - fChi2 /= n; - fStatus = 0; - } - - line LineFit::GetLine() - { - line l; - l.q = fP0; - l.m = fP1; - return l; - } - - LineFit::~LineFit() - { - } - -// --------------------------------------------------------------------------------------------- -// --------------------------------------------------------------------------------------------- - - void CircleExtractor::splitHits(std::vector magnetZs, - std::vector x, std::vector y, std::vector z, - std::vector> &hits) - { - // first step is to sort the separators - std::sort(magnetZs.begin(), magnetZs.end()); - hits.resize(magnetZs.size()+1); - - //scan all recorded hits - for(unsigned int i=0; i x, const std::vector y, const std::vector z) - { - - // two cases depending on the magnet design: - // there can be many magnets, so many magnet positions - // get the magnet position(s) for each geometry (SAMURAI or CrystalPulling design) first - GeometricalParameters::magnetOption opt = GeometricalParameters::Get()->GetFASER2MagnetOption(); - G4int nmag = 0; - std::vector zin; - std::vector zout; - - // STEP 1: get the geometry - if( opt == GeometricalParameters::magnetOption::SAMURAI ){ - - // there is only one magnet, so one magnet position - nmag = 1; - fzpos.push_back(GeometricalParameters::Get()->GetMagnetZPosition()); - double size = GeometricalParameters::Get()->GetMagnetTotalSizeZ(); - zin.push_back(fzpos[0] - size/2.); //begin of magnet window - zout.push_back(fzpos[0] + size/2.); //end of magnet window - - }else if ( opt == GeometricalParameters::magnetOption::CrystalPulling ){ - - double zcenter = GeometricalParameters::Get()->GetMagnetZPosition(); - double size = GeometricalParameters::Get()->GetFASER2MagnetWindowZ(); - double spacing = GeometricalParameters::Get()->GetFASER2MagnetSpacing(); - nmag = GeometricalParameters::Get()->GetNFASER2Magnets(); - for(int i=0; i> hits; - splitHits( fzpos, x, y, z, hits); - - // STEP 3: apply the procedure for each magnet available - // Each magnet uses its own pre and post hits (i and i+1) - // field is along Y: look in z-x plan (for now: FIXME??) - for(int i=0; i zpre; std::vector xpre; - std::vector zpost; std::vector xpost; - - for( unsigned int j=0; j< hits[i].size(); j++){ - zpre.push_back( hits[i][j].z ); - xpre.push_back( hits[i][j].x ); - } - for( unsigned int j=0; j< hits[i+1].size(); j++){ - zpost.push_back( hits[i+1][j].z ); - xpost.push_back( hits[i+1][j].x ); - } - - //and now the actual circle extraction - computeCircle( zpre, xpre, zpost, xpost, zin[i], zout[i]); - } - } - - void CircleExtractor::computeCircle(std::vector zpre, std::vector xpre, std::vector zpost, - std::vector xpost, double Zin, double Zout) - { - LineFit lf_pre(zpre,xpre); - LineFit lf_post(zpost,xpost); - line pre = lf_pre.GetLine(); - line post = lf_post.GetLine(); - fpre.push_back( pre ); - fpost.push_back( post ); - - double Xin = extrapolateLine(Zin,pre); - double Xout = extrapolateLine(Zout,post); - std::pair p = getPerpLineIntersection(Zin,Xin,Zout,Xout,pre,post); - fXc.push_back( p.second ); - fZc.push_back( p.first ); - - double R1 = getR(Zin,Xin,p.first,p.second); - double R2 = getR(Zout,Xout,p.first,p.second); - fR1.push_back( R1 ); - fR2.push_back( R2 ); - } - - double CircleExtractor::extrapolateLine(double z, line l) - { - return l.q + z*l.m; - } - - std::pair CircleExtractor::getPerpLineIntersection(double zin, double xin, double zout, double xout, line l1, line l2) - { - double m1 = l1.m; - double m2 = l2.m; - - double z = m1*m2/(m1-m2)*(xout-xin) + m1*zout/(m1-m2) - m2*zin/(m1-m2); - double x1 = xin - 1./m1*(z-zin); - double x2 = xout - 1./m2*(z-zout); - - return std::make_pair(z,(x1+x2)/2.); - } - - double CircleExtractor::getR(double za, double xa, double zb, double xb) - { - return TMath::Sqrt( (za-zb)*(za-zb) + (xa-xb)*(xa-xb) ); - } - - std::vector CircleExtractor::GetR(){ - std::vector R; - for( unsigned int i=0; i z, std::vector x, double r0) - { - TGraph *g = new TGraph(); - - for(unsigned int i=0; iSetPoint(i,z.at(i),x.at(i)); } - - // parabolic in x-z plane: x = a + b*z + x*z2 - TF1 *f = new TF1("f","[0]+[1]*x+[2]*x*x",z.front()-100,z.back()+100); - f->SetParameter(0,0); - f->SetParameter(1, (x.at(1)-x.at(0))/(z.at(1)-z.at(0)) ); - f->SetParameter(2, 1./(2*r0) ); - - g->Fit("f", "QN"); - - fA = f->GetParameter(0); - fB = f->GetParameter(1); - double c = f->GetParameter(2); - - fR = 1./(2*c); - } - - ParabolicFit::~ParabolicFit() - { - } - -}