diff --git a/doc/developers.md b/doc/developers.md index 8ec1f94b..3ced93b2 100644 --- a/doc/developers.md +++ b/doc/developers.md @@ -62,8 +62,9 @@ mu2e -c EventNtuple/fcl/from_mcs-mockdata.fcl -s test-art-file.art -n 10 8. Add to ```ntuplehelper``` * instructions [here](./ntuplehelper.md#Adding-a-branch) -9. Add to validation, RooUtil, and PyUtil - * [instructions to be completed... for the time being contact Andy and Sophie for this] +9. Add to validation, RooUtil, and pyutils + * instructions for validation and RooUtil are [here](../rooutil/README.md#For-Developers). Contact Andy for more info if needed + * contact Sophie for pyutils 10. Run validation * both [test_fcls.sh](../validation/README.md#Validating-EventNtuple-Runs), and * and [valCompare](../validation/README.md#Validating-EventNtuple-Contents) steps diff --git a/fcl/README.md b/fcl/README.md index dc1838b5..9fe79829 100644 --- a/fcl/README.md +++ b/fcl/README.md @@ -22,6 +22,7 @@ where ```tier``` is the data tier of the input dataset, ```type``` is the type o | from_mcs-mockdata_separateTrkBranches.fcl | mock datasets | example on how to separate the tracks into separate branches again| | from_mcs-mockdata_selectorExample.fcl | mock datasets | example on how to use a selector to select certain types of tracks before putting them into the EventNtuple | | from_mcs-mixed_trkQualCompare.fcl | reconstructed mixed (i.e. primary+background hits) datasets | shows how to output result of more than one TrkQual | +| from_mcs-primary_addVDSteps.fcl | reconstructed primary (i.e. no background hits) datasets | shows how to add the branch for virtual detector steps | | from_mcs-DeMCalib.fcl | reconstructed primary or mixed datasets | only writes one track per event | | from_mcs-OffSpill.fcl | off spill datasets | only contains ```CentralHelix``` tracks (i.e. field-on cosmics) | | from_dig-mockdata.fcl | mock datasets (digis) | runs reconstruction and creates EventNtuple in one job | diff --git a/fcl/from_mcs-primary_addVDSteps.fcl b/fcl/from_mcs-primary_addVDSteps.fcl new file mode 100644 index 00000000..85906f29 --- /dev/null +++ b/fcl/from_mcs-primary_addVDSteps.fcl @@ -0,0 +1,3 @@ +#include "EventNtuple/fcl/from_mcs-primary.fcl" + +physics.analyzers.EventNtuple.StepPointMCTags : [ "compressRecoMCs:virtualdetector" ] diff --git a/fcl/prolog.fcl b/fcl/prolog.fcl index 868f28f1..b18595c8 100644 --- a/fcl/prolog.fcl +++ b/fcl/prolog.fcl @@ -319,6 +319,7 @@ EventNtupleMaker : { ExtraMCStepCollectionTags : [] SurfaceStepCollectionTag : "compressRecoMCs" FitType : LoopHelix + StepPointMCTags : [ ] } # instance for processing trigger (digitization) output from simulation EventNtupleMakerTTMC: { diff --git a/helper/ntuplehelper.py b/helper/ntuplehelper.py index 4136d36b..7e1e67d3 100644 --- a/helper/ntuplehelper.py +++ b/helper/ntuplehelper.py @@ -3,14 +3,15 @@ class nthelper: single_object_branches = ['evtinfo', 'evtinfomc', 'hitcount', 'tcnt', 'crvsummary', 'crvsummarymc'] - vector_object_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'caloclusters', 'calohits', 'calorecodigis', 'calodigis', 'crvcoincs', 'crvcoincsmc', 'crvcoincsmcplane', 'trkqual', 'trkpid'] + vector_object_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'caloclusters', 'calohits', 'calorecodigis', 'calodigis', 'crvcoincs', 'crvcoincsmc', 'crvcoincsmcplane', 'trkqual', 'trkpid', 'mcsteps'] vector_vector_object_branches = ['trksegs', 'trksegpars_lh', 'trksegpars_ch', 'trksegpars_kl', 'trkmcsim', 'trkhits', 'trkhitsmc', 'trkmats', 'trkhitcalibs', 'trkmcsci', 'trkmcssi', 'trksegsmc' ] evt_branches = ['evtinfo','evtinfomc','hitcount','tcnt'] trk_branches = ['trk', 'trkmc', 'trkcalohit', 'trkcalohitmc', 'trkqual', 'trkpid'] trksegs_branches = ['trksegs', 'trksegpars_lh', 'trksegpars_ch', 'trksegpars_kl', 'trksegsmc'] straw_branches = ['trkhits', 'trkmats', 'trkhitsmc', 'trkhitcalibs'] - mc_branches = ['trkmcsim'] + trk_mc_branches = [ 'trkmcsim' ] + general_mc_branches = [ 'mcsteps' ] calo_branches = ['caloclusters', 'calohits', 'calorecodigis', 'calodigis'] crv_branches = ['crvsummary','crvsummarymc','crvcoincs','crvcoincsmc','crvcoincsmcplane'] deprecated_branches = ['trkmcsci','trkmcssi'] @@ -58,7 +59,8 @@ class nthelper: "trkqual" : "MVAResultInfo", "trkpid" : "MVAResultInfo", "helices" : "HelixInfo", - "trksegsmc" : "SurfaceStepInfo" + "trksegsmc" : "SurfaceStepInfo", + "mcsteps" : "MCStepInfo" } # @@ -170,16 +172,34 @@ def list_all_branches(self, export_to_md=False): print("| " + tokens[0] + " | " + tokens[1] + " | " + tokens[2] + "| [see " + struct_file + "](../inc/"+struct_file+")") if not export_to_md: - print("\nMonte Carlo Branches") + print("\nTrk Monte Carlo Branches") print("================") else: - print("## Monte Carlo Branches\n") + print("## Trk Monte Carlo Branches\n") print("These branches contain 4 elements per event corresponding to different Kalman fit hypotheses (see Track branches).\n") print("Within each Kalman fit element, there is a vector containing Monte Carlo truth information about the particle making the track and its parent particles.\n") print("The vector is sorted in reverse chronological order, such that the last element is the initial particle simulated in GEANT4, and each element before correspond to one of its daughter particles.\n") print("| branch | structure | explanation | leaf information |") print("|--------|-----------|-------------|------------------|") - for branch in self.mc_branches: + for branch in self.trk_mc_branches: + explanation = self.get_branch_explanation(branch) + struct = self.branch_struct_dict[branch] + struct_file = struct + ".hh"; + if not export_to_md: + print(explanation) + else: + tokens=explanation.split(":") + print("| " + tokens[0] + " | " + tokens[1] + " | " + tokens[2] + "| [see " + struct_file + "](../inc/"+struct_file+")") + + if not export_to_md: + print("\nGeneral Monte Carlo Branches") + print("================") + else: + print("## General Monte Carlo Branches\n") + print("These branches contain MC information in the event and are not required to be related to tracks") + print("| branch | structure | explanation | leaf information |") + print("|--------|-----------|-------------|------------------|") + for branch in self.general_mc_branches: explanation = self.get_branch_explanation(branch) struct = self.branch_struct_dict[branch] struct_file = struct + ".hh"; diff --git a/inc/InfoMCStructHelper.hh b/inc/InfoMCStructHelper.hh index 63d2d238..ba555ae6 100644 --- a/inc/InfoMCStructHelper.hh +++ b/inc/InfoMCStructHelper.hh @@ -75,7 +75,7 @@ namespace mu2e { void fillExtraMCStepInfos(KalSeedMC const& kseedmc, StepPointMCCollection const& mcsteps, std::vector& mcsics, std::vector& mcssis); void fillSurfaceStepInfos(KalSeedMC const& kseedmc, SurfaceStepCollection const& surfsteps,std::vector& ssic); - + void fillStepPointMCInfo(StepPointMCCollection const& mcsteps, MCStepInfos& mcstepinfos); }; } diff --git a/inc/MCStepInfo.hh b/inc/MCStepInfo.hh index b10a420d..dfaca703 100644 --- a/inc/MCStepInfo.hh +++ b/inc/MCStepInfo.hh @@ -14,12 +14,16 @@ namespace mu2e { float time = 0; // time of this step WRT MC primary proton (ns) float de = 0; // energy deposit through this step (MeV) float dp = 0; // momentum magnitude change throw this step (MeV/c) - bool early = false; - bool late = false; // flag if this is the earliest or latest step + bool early = false; // flag if this is the earliest step + bool late = false; // flag if this is the latest step XYZVectorF mom; // particle momentum at the start of this step XYZVectorF pos; // particle position at the start of this step void reset() {*this = MCStepInfo(); } bool valid() { return vid>=0; } + int pdg = -1; // true PDG code + int startCode = -1; // creation process code + int stopCode = -1; // stop process code + }; using MCStepInfos = std::vector; diff --git a/rooutil/README.md b/rooutil/README.md index e163e63a..7f600dd6 100644 --- a/rooutil/README.md +++ b/rooutil/README.md @@ -114,6 +114,7 @@ Some branches are not contained in any of the above classes: * ```crvcoincsmcplane``` * ```calohits```, ```calorecodigis```, ```calodigis``` * ```trig_``` branches +* ```mcsteps_virtualdetector``` Reach out to the developers on the #analysis-tools Slack channel if you need to have these added somewhere. @@ -224,23 +225,22 @@ Beware: there are a lot of debug messages. If a new branch is added to the EventNtuple, then the following needs to be done to so that RooUtil can access the branch: 1. In [Event.hh](inc/Event.hh) add the pointers at the bottom of the file -2. In [Event.hh](inc/Event.hh) constructor, set the branch address - - make sure to test that the branch exists -3. In [Event.hh](inc/Event.hh) add the #include to the underlying object +2. In [Event.hh](inc/Event.hh) constructor, set the branch address using the ```CheckForBranch()``` function +3. In [Event.hh](inc/Event.hh) add the #include to the underlying info struct 4. In [RooUtil.hh](inc/RooUtil.hh) add it to the ```CreateOutputEventNtuple()``` function 5. Add to validation places: - [PrintEvents.C](examples/PrintEvents.C): at least the first and last leaf in the struct - [create_val_file_rooutil.C](../../validation/create_val_file_rooutil.C) - copy contents of struct into place where histograms are defined - then copy and replace "type " with "TH1F* h_branchname_" (e.g. "float " with "TH1F* h_trkcalohit_") - - then copy in " = new TH1F("h_branchname_", "", 100,0,100); //" + - then copy in " = new TH1F("h_branchname_", "", 100,0,100); //" after the histogram name - delete line after "//" - add in leaf names to histname (and make separate x, y, z histograms for XYZVectorF leaves) - copy histogram lines into main loop - - search and replace "TH1F* " with "" - - search and replace " = new TH1F("h_" with "->Fill(" - - search and replace "", "", 100,0,100);" with ");" - - search and replace "(branchname_" with "(branchname." + - search and replace ```TH1F* ``` (note space) with nothing + - search and replace ``` = new TH1F("h_``` with ```->Fill(``` + - search and replace ```", "", 100,0,100);``` with ```);``` + - search and replace ```(branchname_``` with ```(branchname.``` in the ```Fill()``` commands - update histogram ranges 6. (Optional) If appropriate, add branches to other classes (e.g. Track.hh) and to ```Event::Update()``` - be sure to test it with an example script diff --git a/rooutil/examples/PlotVDSteps.C b/rooutil/examples/PlotVDSteps.C new file mode 100644 index 00000000..bd51978e --- /dev/null +++ b/rooutil/examples/PlotVDSteps.C @@ -0,0 +1,52 @@ +// +// An example of how to plot reco vs true +// This uses cut functions defined in common_cuts.hh +// + +#include "EventNtuple/rooutil/inc/RooUtil.hh" +#include "EventNtuple/rooutil/inc/common_cuts.hh" + +#include "TCanvas.h" +#include "TH1F.h" +#include "TH2F.h" + +using namespace rooutil; +void PlotVDSteps(std::string filename) { + + // Create the histogram you want to fill + TH2F* hVD13_XY = new TH2F("hVD13_XY", "XY position at VD13 (tracker entrance)", 100,-500,500, 100,-500,500); + hVD13_XY->SetXTitle("X (det. coords) [mm]"); + hVD13_XY->SetYTitle("Y (det. coords) [mm]"); + + TH1F* hVD13_Mom = new TH1F("hVD13_Mom", "Momentum at VD13 (tracker entrance)", 120,0,120); + hVD13_Mom->SetXTitle("Momentum [MeV/c]"); + + // Set up RooUtil + RooUtil util(filename); + + // Loop through the events + for (int i_event = 0; i_event < util.GetNEvents(); ++i_event) { + // Get the next event + auto& event = util.GetEvent(i_event); + + // Get the e_minus tracks from the event + if (event.mcsteps_virtualdetector != nullptr) { + for (const auto& vdstep : *(event.mcsteps_virtualdetector)) { + if (vdstep.vid == mu2e::VirtualDetectorId::TT_FrontHollow) { + // Fill the histogram + hVD13_XY->Fill(vdstep.pos.x(), vdstep.pos.y()); + + hVD13_Mom->Fill(vdstep.mom.R()); + } + } + } + } + + // Draw the histogram + TCanvas* c = new TCanvas(); + c->Divide(2); + c->cd(1); + hVD13_XY->Draw("COLZ"); + c->cd(2); + hVD13_Mom->Draw("HIST E"); +} diff --git a/rooutil/examples/PrintEvents.C b/rooutil/examples/PrintEvents.C index a4aa22dd..07d70321 100644 --- a/rooutil/examples/PrintEvents.C +++ b/rooutil/examples/PrintEvents.C @@ -263,5 +263,11 @@ void PrintEvents(std::string filename) { std::cout << pair.first << ": " << event.trigger.Fired(pair.second) << std::endl; } + // mcsteps_virtualdetector branch + if (event.mcsteps_virtualdetector != nullptr) { + for (const auto& vdstep : *(event.mcsteps_virtualdetector)) { + std::cout << "vdstep: " << vdstep.vid << ", " << vdstep.time << ", " << vdstep.mom << ", " << vdstep.pos << ", " << vdstep.pdg << ", " << vdstep.startCode << ", " << vdstep.stopCode << std::endl; + } + } } } diff --git a/rooutil/inc/Event.hh b/rooutil/inc/Event.hh index b400a987..349d113b 100644 --- a/rooutil/inc/Event.hh +++ b/rooutil/inc/Event.hh @@ -37,6 +37,7 @@ #include "EventNtuple/inc/MVAResultInfo.hh" #include "EventNtuple/inc/TrigInfo.hh" +#include "EventNtuple/inc/MCStepInfo.hh" #include "EventNtuple/rooutil/inc/Track.hh" #include "EventNtuple/rooutil/inc/CrvCoinc.hh" @@ -87,6 +88,7 @@ namespace rooutil { CheckForBranch(ntuple, "calorecodigis", &this->calorecodigis); CheckForBranch(ntuple, "calodigis", &this->calodigis); + CheckForBranch(ntuple, "mcsteps_virtualdetector", &this->mcsteps_virtualdetector); } // Check if a branch exists in the TChain, and optionally set its address @@ -342,6 +344,7 @@ namespace rooutil { std::vector* crvcoincsmcplane = nullptr; std::vector>* trkmcsim = nullptr; + std::vector* mcsteps_virtualdetector = nullptr; // TODO: EventNtuple could have other mcsteps_* branches but for the time being just hardcode for the virtualdetector ones }; } // namespace rooutil #endif diff --git a/rooutil/inc/RooUtil.hh b/rooutil/inc/RooUtil.hh index 2f6c9855..3d6d84af 100644 --- a/rooutil/inc/RooUtil.hh +++ b/rooutil/inc/RooUtil.hh @@ -148,6 +148,8 @@ namespace rooutil { output_ntuple->Branch(("trig_" + pair.first).c_str(), &event->triginfo._triggerArray[pair.second]); } + if (event->mcsteps_virtualdetector) { output_ntuple->Branch("mcsteps_virtualdetector", event->mcsteps_virtualdetector); } + // Write out histograms from input to output hVersion->Write(); } diff --git a/rooutil/inc/common_cuts.hh b/rooutil/inc/common_cuts.hh index ef22d70b..5f3928b2 100644 --- a/rooutil/inc/common_cuts.hh +++ b/rooutil/inc/common_cuts.hh @@ -10,6 +10,7 @@ #include "Offline/MCDataProducts/inc/ProcessCode.hh" #include "Offline/MCDataProducts/inc/GenId.hh" #include "Offline/MCDataProducts/inc/MCRelationship.hh" +#include "Offline/DataProducts/inc/VirtualDetectorId.hh" #include "EventNtuple/rooutil/inc/Track.hh" #include "EventNtuple/rooutil/inc/TrackSegment.hh" diff --git a/src/EventNtupleMaker_module.cc b/src/EventNtupleMaker_module.cc index e0562ec4..a9d4af71 100644 --- a/src/EventNtupleMaker_module.cc +++ b/src/EventNtupleMaker_module.cc @@ -182,9 +182,10 @@ namespace mu2e { fhicl::Atom primaryParticleTag{Name("PrimaryParticleTag"), Comment("Tag for PrimaryParticle"), art::InputTag()}; fhicl::Atom kalSeedMCTag{Name("KalSeedMCAssns"), Comment("Tag for KalSeedMCAssn"), art::InputTag()}; // extra MC - fhicl::OptionalSequence extraMCStepTags{Name("ExtraMCStepCollectionTags"), Comment("Input tags for any other StepPointMCCollections you want written out")}; + fhicl::OptionalSequence extraMCStepTags{Name("ExtraMCStepCollectionTags"), Comment("Input tags for any other StepPointMCCollections you want written out. This will only write out steps associated with the KalSeed")}; // passive elements and Virtual Detector MC information fhicl::OptionalAtom SurfaceStepsTag{Name("SurfaceStepCollectionTag"), Comment("SurfaceStep Collection")}; + fhicl::OptionalSequence stepPointMCTags{Name("StepPointMCTags"), Comment("Input tags for any other StepPointMCCollections you want written out. This will write out all steps in the collection")}; // Calo MC fhicl::Atom fillCaloMC{ Name("FillCaloMC"),Comment("Fill CaloMC information"), true}; fhicl::Atom caloClusterMCTag{Name("CaloClusterMCTag"), Comment("Tag for CaloClusterMCCollection") ,art::InputTag()}; @@ -249,7 +250,7 @@ namespace mu2e { std::map _tmap; // map between path and trigger ID. ID should come from trigger itself FIXME! // MC truth (fcl parameters) bool _fillmc; - // MC truth inputs + // MC steps (associated with KalSeed) std::vector _extraMCStepTags; std::vector> _extraMCStepCollections; std::map>> _extraMCStepInfos; @@ -258,6 +259,11 @@ namespace mu2e { art::InputTag _surfaceStepsTag; std::map>> _surfaceStepInfos; art::Handle _surfaceStepsHandle; + // MC steps (all) + std::vector _stepPointMCTags; + std::vector> _stepPointMCCollections; + std::map _stepPointMCInfos; + // art::Handle _pph; art::Handle _ksmcah; @@ -441,6 +447,11 @@ namespace mu2e { _surfaceStepInfos[i_branch] = std::vector>(); } } + if(_conf.stepPointMCTags(_stepPointMCTags)){ + for (StepCollIndex i_stepPointMCTag = 0; i_stepPointMCTag < _stepPointMCTags.size(); ++i_stepPointMCTag) { + _stepPointMCInfos[i_stepPointMCTag] = MCStepInfos(); + } + } } } @@ -566,6 +577,18 @@ namespace mu2e { } // helix info if(_conf.helices()) _ntuple->Branch("helices.",&_hinfos,_buffsize,_splitlevel); + + // all MC information + if (_fillmc) { + if(_conf.stepPointMCTags(_stepPointMCTags)){ + for(size_t icoll=0;icoll < _stepPointMCTags.size(); ++icoll) { + auto const& tag = _stepPointMCTags[icoll]; + auto inst = tag.instance(); + std::string branchname = "mcsteps_" + inst + "."; + _ntuple->Branch(branchname.c_str(),&_stepPointMCInfos[icoll],_buffsize,_splitlevel); + } + } + } } void EventNtupleMaker::beginSubRun(const art::SubRun & subrun ) { @@ -675,6 +698,16 @@ namespace mu2e { } event.getByLabel(_surfaceStepsTag,_surfaceStepsHandle); + // find extra MCStep collections + _stepPointMCCollections.clear(); + for(size_t icoll = 0; icoll < _stepPointMCTags.size(); icoll++){ + auto const& tag = _stepPointMCTags[icoll]; + art::Handle mcstepch; + event.getByLabel(tag,mcstepch); + _stepPointMCCollections.push_back(mcstepch); + } + + // loop through all track types unsigned ntrks(0); for (BranchIndex i_branch = 0; i_branch < _allBranches.size(); ++i_branch) { @@ -709,6 +742,10 @@ namespace mu2e { _extraMCStepSummaryInfos.at(i_branch).at(i_extraMCStepTag).clear(); } _surfaceStepInfos.at(i_branch).clear(); + for (StepCollIndex i_stepPointMCTag = 0; i_stepPointMCTag < _stepPointMCTags.size(); ++i_stepPointMCTag) { + _stepPointMCInfos.at(i_stepPointMCTag).clear(); + } + if(_fillcalomc) { _allMCTCHIs.at(i_branch).clear(); } @@ -875,6 +912,17 @@ namespace mu2e { } } + + // fill MC information unrelated to any reconstructed object + if (_fillmc) { + // fill MCStepCollection branch + for(size_t icoll = 0; icoll < _stepPointMCTags.size(); icoll++){ + auto const& mcsteps = *_stepPointMCCollections[icoll]; + auto& mcstepinfos = _stepPointMCInfos[icoll]; + _infoMCStructHelper.fillStepPointMCInfo(mcsteps, mcstepinfos); + } + } + // fill this row in the TTree bool fill = true; // default to fliling event if(_hastrks && ntrks == 0) { // if we require tracks (_hastrks) but we have none, don't write diff --git a/src/InfoMCStructHelper.cc b/src/InfoMCStructHelper.cc index 3acfe722..b6fd0946 100644 --- a/src/InfoMCStructHelper.cc +++ b/src/InfoMCStructHelper.cc @@ -409,6 +409,9 @@ namespace mu2e { mcsi.dp = mcstep.postMomentum().mag()- mcstep.momentum().mag(); mcsi.mom = mcstep.momentum(); mcsi.pos = spos; + mcsi.pdg = mcstep.simParticle()->pdgId(); + mcsi.startCode = mcstep.simParticle()->creationCode(); + mcsi.stopCode = mcstep.simParticle()->stoppingCode(); } } } @@ -434,4 +437,22 @@ namespace mu2e { // std::sort(ssic.begin(),ssic.end(),[](const auto& a, const auto& b){return a.time < b.time;}); } + void InfoMCStructHelper::fillStepPointMCInfo(StepPointMCCollection const& mcsteps, MCStepInfos& mcstepinfos) { + GeomHandle det; + MCStepInfo mcstepinfo; + for(auto const& mcstep : mcsteps) { + mcstepinfo.reset(); + mcstepinfo.vid = mcstep.volumeId(); + mcstepinfo.time = mcstep.time(); + mcstepinfo.de = mcstep.totalEDep(); + mcstepinfo.dp = mcstep.postMomentum().mag()- mcstep.momentum().mag(); + mcstepinfo.mom = mcstep.momentum(); + mcstepinfo.pos = XYZVectorF(det->toDetector(mcstep.position())); + const auto& simParticle = mcstep.simParticle(); + mcstepinfo.pdg = simParticle->pdgId(); + mcstepinfo.startCode = simParticle->creationCode(); + mcstepinfo.stopCode = simParticle->stoppingCode(); + mcstepinfos.emplace_back(mcstepinfo); + } + } } diff --git a/validation/create_val_file_rooutil.C b/validation/create_val_file_rooutil.C index c172bb26..87a62a25 100644 --- a/validation/create_val_file_rooutil.C +++ b/validation/create_val_file_rooutil.C @@ -505,6 +505,25 @@ void create_val_file_rooutil(std::string filename, std::string outfilename) { TH1F* h_calodigis_caloRecoDigiIdx_ = new TH1F("h_calodigis_caloRecoDigiIdx_", "", 100,0,100); TH1F* h_triginfo = new TH1F("h_triginfo", "", mu2e::TrigInfo::ntrig_,0,mu2e::TrigInfo::ntrig_); + + TH1F* h_mcsteps_virtualdetector_vid = new TH1F("h_mcsteps_virtualdetector_vid", "", 150,0,150); + TH1F* h_mcsteps_virtualdetector_sid = new TH1F("h_mcsteps_virtualdetector_sid", "", 100,0,100); + TH1F* h_mcsteps_virtualdetector_iinter = new TH1F("h_mcsteps_virtualdetector_iinter", "", 100,0,100); + TH1F* h_mcsteps_virtualdetector_time = new TH1F("h_mcsteps_virtualdetector_time", "", 100,0,2000); + TH1F* h_mcsteps_virtualdetector_de = new TH1F("h_mcsteps_virtualdetector_de", "", 100,0,200); + TH1F* h_mcsteps_virtualdetector_dp = new TH1F("h_mcsteps_virtualdetector_dp", "", 100,0,200); + TH1F* h_mcsteps_virtualdetector_early = new TH1F("h_mcsteps_virtualdetector_early", "", 2,0,2); + TH1F* h_mcsteps_virtualdetector_late = new TH1F("h_mcsteps_virtualdetector_late", "", 2,0,2); + TH1F* h_mcsteps_virtualdetector_mom_x = new TH1F("h_mcsteps_virtualdetector_mom_x", "", 100,-200,200); + TH1F* h_mcsteps_virtualdetector_mom_y = new TH1F("h_mcsteps_virtualdetector_mom_y", "", 100,-200,200); + TH1F* h_mcsteps_virtualdetector_mom_z = new TH1F("h_mcsteps_virtualdetector_mom_z", "", 320,-200,200); + TH1F* h_mcsteps_virtualdetector_pos_x = new TH1F("h_mcsteps_virtualdetector_pos_x", "", 100,-1000,1000); + TH1F* h_mcsteps_virtualdetector_pos_y = new TH1F("h_mcsteps_virtualdetector_pos_y", "", 100,-1000,1000); + TH1F* h_mcsteps_virtualdetector_pos_z = new TH1F("h_mcsteps_virtualdetector_pos_z", "", 320,-1600,1600); + TH1F* h_mcsteps_virtualdetector_pdg = new TH1F("h_mcsteps_virtualdetector_pdg", "", 100,0,100); + TH1F* h_mcsteps_virtualdetector_startCode = new TH1F("h_mcsteps_virtualdetector_startCode", "", 100,0,100); + TH1F* h_mcsteps_virtualdetector_stopCode = new TH1F("h_mcsteps_virtualdetector_stopCode", "", 100,0,100); + // Grab the first event to name the bins const auto& evt = util.GetEvent(0); for (const auto pair : evt.trigger.NameToIndexMap()) { @@ -1153,6 +1172,30 @@ void create_val_file_rooutil(std::string filename, std::string outfilename) { h_triginfo->AddBinContent(i_trig+1, event.triginfo._triggerArray[i_trig]); } + // Creating mcsteps_virtualdetector histogram + std::cout << "Creating mcsteps_virtualdetector histogram" << std::endl; + if (event.mcsteps_virtualdetector != nullptr) { + for (const auto& vdstep : *(event.mcsteps_virtualdetector)) { + h_mcsteps_virtualdetector_vid->Fill(vdstep.vid); + h_mcsteps_virtualdetector_sid->Fill(vdstep.sid); + h_mcsteps_virtualdetector_iinter->Fill(vdstep.iinter); + h_mcsteps_virtualdetector_time->Fill(vdstep.time); + h_mcsteps_virtualdetector_de->Fill(vdstep.de); + h_mcsteps_virtualdetector_dp->Fill(vdstep.dp); + h_mcsteps_virtualdetector_early->Fill(vdstep.early); + h_mcsteps_virtualdetector_late->Fill(vdstep.late); + h_mcsteps_virtualdetector_mom_x->Fill(vdstep.mom.x()); + h_mcsteps_virtualdetector_mom_y->Fill(vdstep.mom.y()); + h_mcsteps_virtualdetector_mom_z->Fill(vdstep.mom.z()); + h_mcsteps_virtualdetector_pos_x->Fill(vdstep.pos.x()); + h_mcsteps_virtualdetector_pos_y->Fill(vdstep.pos.y()); + h_mcsteps_virtualdetector_pos_z->Fill(vdstep.pos.z()); + h_mcsteps_virtualdetector_pdg->Fill(vdstep.pdg); + h_mcsteps_virtualdetector_startCode->Fill(vdstep.startCode); + h_mcsteps_virtualdetector_stopCode->Fill(vdstep.stopCode); + + } + } } file->Write(); diff --git a/validation/test_fcls.sh b/validation/test_fcls.sh index 6f05f89a..6792ad42 100644 --- a/validation/test_fcls.sh +++ b/validation/test_fcls.sh @@ -2,6 +2,10 @@ # test_fcls.sh - runs fcl files to make sure that they complete successfully # Note: requires relevant filelists in a directory above this one # +if [ $# -ne 1 ]; then + echo "Usage: . ./validation/test_fcls.sh MDC202X" + return +fi vMDC=$1 log_file="test_fcls.log" @@ -137,6 +141,16 @@ else echo "FAIL" fi +echo -n "from_mcs-primary_addVDSteps.fcl... " +echo "mu2e -c fcl/from_mcs-primary_addVDSteps.fcl -S ../filelists/$primary_dataset.list --TFileName nts.ntuple.primaryVDSteps.root -n 100" >> ${log_file} 2>&1 +mu2e -c fcl/from_mcs-primary_addVDSteps.fcl -S ../filelists/$primary_dataset.list --TFileName nts.ntuple.primaryVDSteps.root -n 100 >> ${log_file} 2>&1 +if [ $? == 0 ]; then + echo "OK" +else + echo "FAIL" +fi + + echo -n "from_dig-mockdata.fcl... " echo "mu2e -c fcl/from_dig-mockdata.fcl -S ../filelists/${digi_dataset}.list --TFileName nts.ntuple.dig.root -n 100" >> ${log_file} 2>&1 mu2e -c fcl/from_dig-mockdata.fcl -S ../filelists/${digi_dataset}.list --TFileName nts.ntuple.dig.root -n 100 >> ${log_file} 2>&1 diff --git a/validation/test_rooutil_examples.sh b/validation/test_rooutil_examples.sh index 852f6ffb..eb74676c 100644 --- a/validation/test_rooutil_examples.sh +++ b/validation/test_rooutil_examples.sh @@ -1,10 +1,10 @@ -filename="\"/pnfs/mu2e/tape/phy-nts/nts/mu2e/CeEndpointOnSpillTriggered/MDC2020aq_best_v1_3_v06_03_00/root/ed/df/nts.mu2e.CeEndpointOnSpillTriggered.MDC2020aq_best_v1_3_v06_03_00.001210_00000000.root\"" +filename="\"nts.ntuple.after.root\"" scripts=( "PrintEvents.C" "CreateNtuple.C" "CreateTrackNtuple.C" "PlotCRVPEs.C" "PlotCRVPEsVsMCEDep.C" "PlotEntranceFitPars.C" "PlotEntranceMomentum.C" "PlotEntranceMomentumCRVCut.C" "PlotEntranceMomentumResolution.C" "PlotEntranceMomentumResolution_TrkQualCut.C" "PlotMCParentPosZ.C" "PlotMCParticleMom.C" "PlotMuonPosZ.C" "PlotStoppingTargetFoilSegment.C" "PlotTrackNHits_RecoVsTrue.C" "PlotTrkCaloHitEnergy.C" "PrintEventsNoMC.C" "TrackCounting.C" "PlotTrackHitTimes.C" "PlotTrackHitTimesMC.C" "PlotStrawMaterials.C" - "PlotGenCosmicMom.C" "PlotCRVTotalPEs.C" "PlotEntranceMomentum_UpstreamDownstream.C" ) + "PlotGenCosmicMom.C" "PlotCRVTotalPEs.C" "PlotEntranceMomentum_UpstreamDownstream.C" "PlotVDSteps.C" ) for script in "${scripts[@]}" do