From 12a9da270cf04db13933fe22326603e8338dd1af Mon Sep 17 00:00:00 2001 From: Javier Galan Date: Mon, 20 Mar 2023 19:37:13 +0100 Subject: [PATCH] Adding a validation for cosmic muon generator --- examples/04.MuonScan/ValidateCosmicMuons.C | 29 ++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/examples/04.MuonScan/ValidateCosmicMuons.C b/examples/04.MuonScan/ValidateCosmicMuons.C index 02b87a8f..a13eca41 100644 --- a/examples/04.MuonScan/ValidateCosmicMuons.C +++ b/examples/04.MuonScan/ValidateCosmicMuons.C @@ -5,6 +5,7 @@ Int_t ValidateCosmicMuons(const char* filename) { TRestRun run(filename); TRestGeant4Event* event = run.GetInputEvent(); + TRestGeant4Metadata* md = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata"); if (run.GetRunTag() != "CosmicMuons") { cout << "Run tag: " << run.GetRunTag() << endl; @@ -29,11 +30,39 @@ Int_t ValidateCosmicMuons(const char* filename) { constexpr double tolerance = 0.1; + Int_t detDown = 0; + Int_t detUp = 0; + Int_t detBoth = 0; for (Long64_t i = 0; i < run.GetEntries(); i++) { run.GetEntry(i); averageTotalEnergy += event->GetTotalDepositedEnergy() / nEvents; averageSensitiveEnergy += event->GetSensitiveVolumeEnergy() / nEvents; + + Bool_t down = !TMath::IsNaN(event->GetFirstPositionInVolume(md->GetActiveVolumeID("det_dw_01")).X()); + Bool_t up = !TMath::IsNaN(event->GetFirstPositionInVolume(md->GetActiveVolumeID("det_up_01")).X()); + if (down) detDown++; + if (up) detUp++; + if (down && up) detBoth++; + } + + cout << "Number of events that crossed down detector : " << detDown << endl; + cout << "Number of events that crossed up detector : " << detUp << endl; + cout << "Number of events that crossed both detectors : " << detBoth << endl; + + if (detDown != 10000) { + cout << "The number of cosmics crossing the sensitive volume is not 10000!" << endl; + return 13; + } + + if (detUp != 462) { + cout << "The number of cosmics crossing the up volume is not 462!" << endl; + return 14; + } + + if (detBoth != 462) { + cout << "The number of cosmics crossing both volumes is not 462!" << endl; + return 15; } cout << "Average total energy: " << averageTotalEnergy << " keV" << endl;