diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000000..e643e93fae --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,150 @@ +# Changelog + +All notable changes to this project will be documented in this file. +The latest changes are placed on top get assigned to the corresponding +software release. + +The software versioning scheme is vMAJOR.MINOR.PATCH+DATE-comment, where + - PATCH - quick bug fixes + - MINOR - new classes or functionality, updates to adapt to new external packages + - MAJOR - additions of the detectors, breaking changes - hopefully we won't have such + - DATE - to ease communication and mimic our CVMFS releases + - comment - a short description of changes if available + - e.g. v1.0.0+2025-07-updateScifi + + +The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/). +There could be several types of changes: +- `Added` for new features. +- `Changed` for changes in existing functionality. +- `Deprecated` for soon-to-be removed features. +- `Removed` for now removed features. +- `Fixed` for any bug fixes. +- `Security` in case of vulnerabilities. + +We start with the first sndsw release: v1.0.0+2025-07-updateScifi. +Shall there be a strong will/need, one can go back and create and +fill in the logs for previous stacks. + +## v1.4.0+2026-01-G4cutoff + +### Added + +- Geant4 range cut set to 2mm +- Record simulation run and event ID +- `is_mc` option for sndConfiguration, used in all sndPlane classes +- Veto and DS plane classes +- 2025 data file locations for the analysis tools + +### Changed + +- Geant4 physics list to FTFP_BERT_HP_EMZ +- retire small SiPMs: add a flag to skip(default) them in analysis + +### Removed + +- Geant4 production cut at 1MeV +- hard coded digi params in run_digi macro + +### Fixed + +- MeV to qdc conversion for US +- Scifi station 5 rotation for target 256 +- revised the drift velocity and attenuation length of Veto/US/DS using 2022-2025 data + + +## v1.3.0+2025-11-showerToolsAndP8Decayer + +### Added + +- shower tools: start, stop, shower direction +- 2dED: Add option to plot shower direction +- python script to make run lists with DB info +- SciFi and US anisotropy tools +- USPlane: GetNHitBars() returting N bars with fired SiPMs above threshold +- MuFilterHit: Add GetImpactXpos method +- add time window selection for US(test_beam_2023 configuration) +- alignment of targets 256-258 + +### Changed + +- G4conf: Remove obsolete neutron and hadronic process verbosity settings +- improve centroid error: use propagation of uncertainty respecting the weight(=qdc) per hit. +- USPlane: Utilize left/right measurements to determine X position along bar +- USPlane: use only small SiPMs to determine if plane HasShower() +- set default csv tables for geofiles and data ROOT files + - these changes in are not backward-compatible, but only affect analysis code that is not part of sndsw +- remove invalid SciFi hits from ScifiPlane's constructor +- 2dED: Plot tracks starting from most upstream Veto + +### Fixed + +- change pythia8 decayer configuration method + - heavy mesons and taus are decay’ed using external Pythia8 decayer +- USPlane: Fix FindCentroid method +- Fix reference SciFi time intervals to values used in prev analysis +- MuFilterHit: fix GetImpactT method +- multiple fixes due to drop of attribute pythonization for TDirectory and TClonesArray +- Disable web-viewer when checking for geo overlaps (workaround!) + +## v1.2.1+2025-09 + +### Added + +- alignment of target 255 + +### Changed + +- optimize the MCEventBuider: + - break loop over chunks after first one if using option saveFirst25nsOnly + - in NC neutrino events, the outgoing neutrino is now saved in the first chunk + +### Fixed + +- fix in SciFi and US plane classes to allow channel indexing in MC and data alike + +## v1.2.0+2025-09-MCEventBuilder + +### Added + +- MCEventBuilder Fair Task + +## v1.1.0+2025-08-BolognaTools + +### Added + +- analysys classes for SciFi and US planes, together with parameter configuraions for TI18 and testbeams setups +- SciFi tools: hit density + +## v1.0.1+2025-07-fixTB24Wtarget + +### Fixed + +- flag logic to select target geometry for testbeam 2024 +- unit fix for selectScifiHits's `time_lower_range` default + +### Added + +- options for the 2dEvent display: filter hits in time, plot collision axis and more. + +## v1.0.0+2025-07-updateScifi + +### Added + +- CHANGELOG +- Data file and geometry getters from run number + +### Changed + +- Major update on the TI18 and baby SciFi detector geometry +- Update on the testbeams system geometry to take care of air gaps + +### Fixed + +- Units for energy and snd_TDC2ns in the cpp units header +- storage of Emulsion(using a flag) and Veto points(always) +- Adapted to the latest FairRoot v19 + +### Removed + +- TI18 MC geometry config file is gone. We use a single one for data and simulation diff --git a/analysis/makeRunListDB.py b/analysis/makeRunListDB.py new file mode 100644 index 0000000000..4aef4bbc13 --- /dev/null +++ b/analysis/makeRunListDB.py @@ -0,0 +1,153 @@ +import pymongo +import argparse +import xml.etree.ElementTree as ET +from datetime import datetime +from collections import defaultdict +import os + +parser = argparse.ArgumentParser( + prog="makeRunListDB", + description="Extracts a list of runs from the SND@LHC DB, matching the conditions given in the arguments. Produces an XML file with the run list and a summary of the selection.") +parser.add_argument("--username", type=str, help="Database username. Contact database maintainer to obtain it", required=True) +parser.add_argument("--password", type=str, help="Database password. Contact database maintainer to obtain it", required=True) +parser.add_argument("--name", type=str, help="Run list name", required=True) +parser.add_argument("--years", nargs="+", type=int, help="Years to be included, e.g. 2022 2023", required=True) +parser.add_argument("--min_events", type=int, help="Minimum number of events in run", default=0) +parser.add_argument("--min_lumi", type=float, help="Minimum integrated luminosity for a run to be included, in fb-1", default=0.) +parser.add_argument("--min_stable_time", type=float, help="Minimum stable beams time of the corresponding LHC fill for a run to be included, in minutes", default=0.) +parser.add_argument("--particle_B1", type=str, help="Beam 1 particle, e.g. p+ or PB82", required=True) +parser.add_argument("--particle_B2", type=str, help="Beam 2 particle, e.g. p+ or PB82", required=True) +parser.add_argument("--min_energy", type=float, help="Minimum beam energy (in GeV)", default=0.) +parser.add_argument("--max_energy", type=float, help="Maximum beam energy (in GeV)", default=0.) +parser.add_argument("--min_bunches_IP1", type=int, help="Minimum number of bunches colliding at IP1", default=0) +parser.add_argument("--max_bunches_IP1", type=int, help="Maximum number of bunches colliding at IP1", default=0) +parser.add_argument("--include_runs", nargs="+", type=int, help="Runs to include, regardless of the other criteria", default=[]) +parser.add_argument("--exclude_runs", nargs="+", type=int, help="Runs to exclude from the list", default=[]) +args = parser.parse_args() + +client = pymongo.MongoClient('sndrundb.cern.ch', username=args.username, password=args.password, authSource='sndrundb', authMechanism='SCRAM-SHA-1') + +db = client.sndrundb + +pipeline = [] + +# Get the run list corresponding to the selected years +pipeline.append({"$match": {"$expr": {"$in": [{"$year": "$start"}, args.years]}}}) + +# Select runs with a minimum of min_events events: +if args.min_events > 0: + pipeline.append({"$match": {"events": {"$gt": args.min_events}}}) + +# Combine data fill from the LPC +pipeline.append({"$lookup": {"from": "FILL_LPC", "localField": "fill", "foreignField": "_id", "as": "LPC"}}) + +if args.min_lumi > 0.: + # Select runs with at least min_lumi integrated luminosity + pipeline.append({"$match": {"LPC.ATLAS_Int_Lumi": {"$gt": args.min_lumi*1e6}}}) + +if args.min_stable_time > 0.: + # Select runs with at least min_stable_time minutes of stable beams + pipeline.append({"$match": {"LPC.Stable_time": {"$gt": args.min_stable_time/60.}}}) + +# Select B1 particle +pipeline.append({"$match": {"LPC.Particle_B1": args.particle_B1}}) + +# Select B2 particle +pipeline.append({"$match": {"LPC.Particle_B2": args.particle_B2}}) + +if args.min_energy > 0: + # Select runs with args.min_energy energy + pipeline.append({"$match": {"LPC.Energy": {"$gt": args.min_energy}}}) + +if args.max_energy > 0: + # Select runs with args.max_energy energy + pipeline.append({"$match": {"LPC.Energy": {"$lt": args.max_energy}}}) + +if args.min_bunches_IP1 > 0: + # Select runs with at least min_bunches_IP1 bunches colliding at IP1 + pipeline.append({"$match": {"LPC.Coll_IP_1_5": {"$gt": args.min_bunches_IP1}}}) + +if args.max_bunches_IP1 > 0: + # Select runs with at most max_bunches_IP1 bunches colliding at IP1 + pipeline.append({"$match": {"LPC.Coll_IP_1_5": {"$lt": args.max_bunches_IP1}}}) + +if len(args.exclude_runs) > 0: + pipeline.append({"$match": {"runNumber": {"$nin" : args.exclude_runs}}}) + +# Expression for Calculating run length from start and stop datetimes +run_length_expr = {"$dateDiff": {"startDate": "$start", "endDate": "$stop", "unit": "minute"}} + +# Extract the following data from the DB +projection = {"$project":{"_id": 0, # Do not extract DB entry ID + "run_number": "$runNumber", # Run number + "n_events": "$events", # Number of events + "start": 1, # Start date + "end": 1, # End date + "duration": run_length_expr, # Run duration + "n_files": {"$size": "$files"}, # Number of files + "path": {"$first": "$files.file"}, # Path of the first file + "fill_number": "$fill", # Fill number + "fill_int_lumi": {"$first": "$LPC.ATLAS_Int_Lumi"}, # Integrated luminosity + "fill_stable_time" : {"$first": "$LPC.Stable_time"}} # Stable beams duration + } + +pipeline.append(projection) + +result = list(db["EcsData"].aggregate(pipeline)) + +if len(args.include_runs) > 0: + include_pipeline = [] + include_pipeline.append({"$match": {"runNumber": {"$in": args.include_runs}}}) + include_pipeline.append(projection) + + include_result = list(db["EcsData"].aggregate(include_pipeline)) + + result.extend(list(include_result)) + +result.sort(key=lambda x: x["run_number"]) + +# Get the time now +now = datetime.now() + +# Format into xml tree +root = ET.Element("runlist") + +meta_data = ET.SubElement(root, "meta") +ET.SubElement(meta_data, "name").text = args.name +ET.SubElement(meta_data, "datetime").text = now.strftime("%Y-%m-%dT%H:%M:%S.%f") +selection = ET.SubElement(meta_data, "selection") +ET.SubElement(selection, "years").text = ','.join([str(y) for y in args.years]) +for criterion in ["min_events", "min_lumi", "min_stable_time", "particle_B1", "particle_B2", "min_energy", "max_energy", "min_bunches_IP1", "max_bunches_IP1", "exclude_runs", "include_runs"]: + ET.SubElement(selection, criterion).text = str(getattr(args, criterion)) +runs = ET.SubElement(root, "runs") + +# Counters for summary statistics +n_runs = 0 +totals = defaultdict(int) + +# Run loop +for run in result: + this_run = ET.SubElement(runs, "run") + totals["n_runs"] += 1 + for field_name in ["run_number", "start", "end", "n_events", "duration", "n_files", "fill_number", "fill_int_lumi", "fill_stable_time", "path"]: + try: + data = run[field_name] + except KeyError: + continue + + if field_name == "path": + data = os.path.dirname(data) + + ET.SubElement(this_run, field_name).text = str(data) + + if field_name not in ["run_number", "start", "end", "fill_number", "path", "fill_int_lumi", "fill_stable_time"] and data is not None: + totals["tot_"+field_name] += data + +stats = ET.SubElement(meta_data, "statistics") +for key, data in totals.items(): + ET.SubElement(stats, key).text = str(data) + +# Write to xml file +tree = ET.ElementTree(root) +ET.indent(tree, space=" ") +tree.write(args.name+"_"+str(now.timestamp())+".xml", encoding="utf-8", xml_declaration=True) diff --git a/analysis/tools/AnalysisToolsLinkDef.h b/analysis/tools/AnalysisToolsLinkDef.h index f12356de71..36d9b7250e 100644 --- a/analysis/tools/AnalysisToolsLinkDef.h +++ b/analysis/tools/AnalysisToolsLinkDef.h @@ -7,6 +7,12 @@ #pragma link C++ nestedclasses; #pragma link C++ nestedtypedef; +#pragma link C++ class snd::Configuration+; +#pragma link C++ class snd::analysis_tools::VetoPlane+; +#pragma link C++ class snd::analysis_tools::ScifiPlane+; +#pragma link C++ class snd::analysis_tools::USPlane+; +#pragma link C++ class snd::analysis_tools::DSPlane+; + #pragma link C++ namespace snd::analysis_tools; #pragma link C++ defined_in namespace snd::analysis_tools; @@ -28,11 +34,23 @@ #pragma link C++ function snd::analysis_tools::densityCheck(const TClonesArray &, int, int, int, bool); #pragma link C++ function snd::analysis_tools::showerInteractionWall(const TClonesArray &,const std::map &, int, std::string); #pragma link C++ function snd::analysis_tools::showerInteractionWall(const TClonesArray &, int, std::string); -#pragma link C++ function snd::analysis_tools::GetDataBasePath(const std::string &, int); -#pragma link C++ function snd::analysis_tools::GetTChain(const std::string &, int, int); -#pragma link C++ function snd::analysis_tools::GetTChain(std::string); -#pragma link C++ function snd::analysis_tools::GetGeoPath(const std::string &,int); -#pragma link C++ function snd::analysis_tools::GetGeometry(std::string); -#pragma link C++ function snd::analysis_tools::GetGeometry(const std::string &,int); +#pragma link C++ function snd::analysis_tools::GetDataBasePath(int, std::string); +#pragma link C++ function snd::analysis_tools::GetTChain(int, int, const std::string &); +#pragma link C++ function snd::analysis_tools::GetTChain(const std::string &); +#pragma link C++ function snd::analysis_tools::GetGeoPath(int, std::string); +#pragma link C++ function snd::analysis_tools::GetGeometry(int, const std::string &); +#pragma link C++ function snd::analysis_tools::GetGeometry(const std::string &); +#pragma link C++ function snd::analysis_tools::FillVeto(const snd::Configuration &, TClonesArray *, MuFilter *); +#pragma link C++ function snd::analysis_tools::FillScifi(const snd::Configuration &, TClonesArray *, Scifi *); +#pragma link C++ function snd::analysis_tools::FillUS(const snd::Configuration &, TClonesArray *, MuFilter *, bool); +#pragma link C++ function snd::analysis_tools::FillDS(const snd::Configuration &, TClonesArray *, MuFilter *); +#pragma link C++ function snd::analysis_tools::GetScifiShowerStart(const std::vector &); +#pragma link C++ function snd::analysis_tools::GetUSShowerEnd(const std::vector &); +#pragma link C++ function snd::analysis_tools::GetShowerInterceptAndDirection(const snd::Configuration &, const std::vector &, const std::vector &); +#pragma link C++ function snd::analysis_tools::GetShoweringPlanes(const std::vector &, const std::vector &); +#pragma link C++ function snd::analysis_tools::GetSciFiSpatialAnisotropy(const std::vector &, bool); +#pragma link C++ function snd::analysis_tools::GetUSSpatialAnisotropy(const std::vector &, bool); #endif diff --git a/analysis/tools/CMakeLists.txt b/analysis/tools/CMakeLists.txt index aed9e2e5e8..6f7e12556e 100644 --- a/analysis/tools/CMakeLists.txt +++ b/analysis/tools/CMakeLists.txt @@ -13,6 +13,13 @@ set(SRCS sndSciFiTools.cxx sndTchainGetter.cxx sndGeometryGetter.cxx +sndConfiguration.cxx +sndVetoPlane.cxx +sndScifiPlane.cxx +sndUSPlane.cxx +sndDSPlane.cxx +sndPlaneTools.cxx +sndShowerTools.cxx ) set(LINK_DIRECTORIES @@ -29,4 +36,4 @@ Set(DEPENDENCIES shipLHC) GENERATE_LIBRARY() -target_link_libraries(snd_analysis_tools -lROOTTPython) \ No newline at end of file +target_link_libraries(snd_analysis_tools -lROOTTPython -lMatrix -lMathCore -lCore) \ No newline at end of file diff --git a/analysis/tools/data_paths.csv b/analysis/tools/data_paths.csv index f91e4da5f4..b68c64c326 100644 --- a/analysis/tools/data_paths.csv +++ b/analysis/tools/data_paths.csv @@ -1,5 +1,5 @@ min_run_number,max_run_number,path -0,5421,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2022/ +4361,5421,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2022/ 5422,7356,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2023/ 7649,8317,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2024/run_241/ 8318,8580,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2024/run_242/ @@ -13,5 +13,13 @@ min_run_number,max_run_number,path 9692,9880,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2024/run_2410/ 9881,10011,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2024/run_2411/ 10012,10422,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2024/run_2412/ +10919,11156,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/run_251/ +11158,11575,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/run_252/ +11576,11675,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/run_253/ +11676,11794,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/run_254/ +11795,12019,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/run_255/ +12021,12121,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/run_256/ +12123,12411,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/run_257/ +12412,12792,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/run_258/ 100238,100679,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/commissioning/testbeam_June2023_H8/ 100841,100985,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/commissioning/testbeam_24/ \ No newline at end of file diff --git a/analysis/tools/geo_paths.csv b/analysis/tools/geo_paths.csv index f682f2bbd4..4588ed00b0 100644 --- a/analysis/tools/geo_paths.csv +++ b/analysis/tools/geo_paths.csv @@ -1,6 +1,8 @@ min_run_number,max_run_number,path -0,7356,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2023/geofile_sndlhc_TI18_V4_2023.root +4361,5422,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2022/geofile_sndlhc_TI18_V4_2022.root +5482,7356,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2023/geofile_sndlhc_TI18_V3_2023.root 7357,10422,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2024/geofile_sndlhc_TI18_V12_2024.root +10919,12792,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/physics/2025/geofile_sndlhc_TI18_V8_2025.root 100238,100679,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/commissioning/testbeam_June2023_H8/geofile_sndlhc_H8_2023_3walls.root -100841,100953,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/commissioning/testbeam_24/geofile_sndlhc_H4_2024_W_2walls_v2.root +100841,100953,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/commissioning/testbeam_24/geofile_sndlhc_H4_2024_W_2walls.root 100954,100985,root://eospublic.cern.ch//eos/experiment/sndlhc/convertedData/commissioning/testbeam_24/geofile_sndlhc_H4_2024_Fe_1wall.root diff --git a/analysis/tools/sndConfiguration.cxx b/analysis/tools/sndConfiguration.cxx new file mode 100644 index 0000000000..398e44cffc --- /dev/null +++ b/analysis/tools/sndConfiguration.cxx @@ -0,0 +1,185 @@ +#include "sndConfiguration.h" + +#include +#include +#include + +#include "Scifi.h" +#include "MuFilter.h" + +snd::Configuration::Configuration(Option option, Scifi *scifi_geometry, MuFilter *muon_filter_geometry, bool is_MC) : is_mc(is_MC) +{ + // Parameters from geometry + scifi_n_stations = scifi_geometry->GetConfParI("Scifi/nscifi"); + scifi_boards_per_plane = scifi_geometry->GetConfParI("Scifi/nmats"); + scifi_n_channels_per_plane = scifi_geometry->GetConfParI("Scifi/nsipm_channels") * scifi_geometry->GetConfParI("Scifi/nmats") * scifi_geometry->GetConfParI("Scifi/nsipm_mat"); + scifi_fiber_lenght = scifi_geometry->GetConfParF("Scifi/fiber_length"); + scifi_centroid_error_x = scifi_geometry->GetConfParF("Scifi/channel_width"); + scifi_centroid_error_y = scifi_geometry->GetConfParF("Scifi/channel_width"); + scifi_centroid_error_z = scifi_geometry->GetConfParF("Scifi/channel_height"); + + veto_n_stations = muon_filter_geometry->GetConfParI("MuFilter/NVetoPlanes"); + veto_bar_per_station = muon_filter_geometry->GetConfParI("MuFilter/NVetoBars"); + + us_n_stations = muon_filter_geometry->GetConfParI("MuFilter/NUpstreamPlanes"); + us_bar_per_station = muon_filter_geometry->GetConfParI("MuFilter/NUpstreamBars"); + us_n_sipm_per_bar = muon_filter_geometry->GetConfParI("MuFilter/UpstreamnSiPMs") * muon_filter_geometry->GetConfParI("MuFilter/UpstreamnSides"); + us_n_channels_per_station = us_bar_per_station * us_n_sipm_per_bar; + us_bar_length = muon_filter_geometry->GetConfParF("MuFilter/UpstreamBarX"); + us_signal_speed = muon_filter_geometry->GetConfParF("MuFilter/VandUpPropSpeed"); + us_centroid_error_x = 3.;//estimate based on PMU data, no time calibration + us_centroid_error_y = muon_filter_geometry->GetConfParF("MuFilter/UpstreamBarY") / std::sqrt(12); + us_centroid_error_z = muon_filter_geometry->GetConfParF("MuFilter/UpstreamBarZ"); + + ds_n_stations = muon_filter_geometry->GetConfParI("MuFilter/NDownstreamPlanes"); + ds_bar_per_station = muon_filter_geometry->GetConfParI("MuFilter/NDownstreamBars"); + ds_hor_spatial_resolution_x = muon_filter_geometry->GetConfParF("MuFilter/DownstreamBarX") / std::sqrt(12); + ds_hor_spatial_resolution_y = muon_filter_geometry->GetConfParF("MuFilter/DownstreamBarY") / std::sqrt(12); + ds_hor_spatial_resolution_z = muon_filter_geometry->GetConfParF("MuFilter/DownstreamBarZ") / std::sqrt(12); + ds_ver_spatial_resolution_x = muon_filter_geometry->GetConfParF("MuFilter/DownstreamBarX_ver") / std::sqrt(12); + ds_ver_spatial_resolution_y = muon_filter_geometry->GetConfParF("MuFilter/DownstreamBarY_ver") / std::sqrt(12); + ds_ver_spatial_resolution_z = muon_filter_geometry->GetConfParF("MuFilter/DownstreamBarZ_ver") / std::sqrt(12); + + // Common parameters not present in geometry + veto_min_hit_on_bar = 0; + scifi_min_timestamp = std::nan(""); + scifi_max_timestamp = std::nan(""); + us_min_timestamp = std::nan(""); + us_max_timestamp = std::nan(""); + us_min_n_hits_for_centroid = 15; + us_qdc_to_gev = 0.0151; + us_min_hit_on_bar = 5; + + // Ad hoc parameters not present in geometry + if (option == Option::ti18_2024_2025) + { + scifi_shower_window_width = 128; + scifi_min_hits_in_window = 10; + scifi_min_n_hits_for_centroid = 5; + scifi_qdc_to_gev = 0.14; + + scifi_x_min = -50.0; + scifi_x_max = 0.0; + scifi_y_min = 10.0; + scifi_y_max = 60.0; + scifi_z_min = 280.0; + scifi_z_max = 360.0; + + us_x_min = -80.0; + us_x_max = 0.0; + us_y_min = 0.0; + us_y_max = 80.0; + us_z_min = 370.0; + us_z_max = 480.0; + + centroid_min_valid_station = 2; + + scifi_min_timestamp = -0.5; + scifi_max_timestamp = 1.2; + } + + else if (option == Option::ti18_2022_2023) + { + scifi_shower_window_width = 128; + scifi_min_hits_in_window = 10; + scifi_min_n_hits_for_centroid = 5; + scifi_qdc_to_gev = 0.14; + scifi_x_min = -50.0; + scifi_x_max = 0.0; + scifi_y_min = 10.0; + scifi_y_max = 60.0; + scifi_z_min = 280.0; + scifi_z_max = 360.0; + + us_x_min = -80.0; + us_x_max = 0.0; + us_y_min = 0.0; + us_y_max = 80.0; + us_z_min = 370.0; + us_z_max = 480.0; + + centroid_min_valid_station = 2; + + scifi_min_timestamp = -0.5; + scifi_max_timestamp = 1.2; + } + + else if (option == Option::test_beam_2023) + { + scifi_shower_window_width = 128; + scifi_min_hits_in_window = 36; + scifi_min_n_hits_for_centroid = 0; + scifi_qdc_to_gev = 0.053; // mirrors + + scifi_x_min = -47.0; + scifi_x_max = -27.0; + scifi_y_min = 35.0; + scifi_y_max = 55.0; + scifi_z_min = 310.0; + scifi_z_max = 360.0; + scifi_min_timestamp = -0.5; + scifi_max_timestamp = 0.5; + + us_x_min = -80.0; + us_x_max = 5.0; + us_y_min = 10.0; + us_y_max = 80.0; + us_z_min = 370.0; + us_z_max = 480.0; + us_min_timestamp = -0.5; + us_max_timestamp = 3.0; + + centroid_min_valid_station = 0; + } + else if (option == Option::test_beam_2024) + { + scifi_shower_window_width = 128; + scifi_min_hits_in_window = 36; + scifi_min_n_hits_for_centroid = 0; + scifi_qdc_to_gev = 0.14; + + scifi_x_min = -47.0; + scifi_x_max = -27.0; + scifi_y_min = 35.0; + scifi_y_max = 55.0; + scifi_z_min = 310.0; + scifi_z_max = 380.0; + + us_x_min = std::nan(""); + us_x_max = std::nan(""); + us_y_min = std::nan(""); + us_y_max = std::nan(""); + us_z_min = std::nan(""); + us_z_max = std::nan(""); + + centroid_min_valid_station = 0; + } + else + { + throw std::invalid_argument("Unknown configuration option"); + } +} + +snd::Configuration::Option snd::Configuration::GetOption(int run_number) +{ + if (run_number >= 100841 && run_number <= 100985) { + std::cout << "Choosing option >>>>>>>>>>\t test_beam_2024 \t<<<<<<<<<<" <= 100238 && run_number <= 100679) { + std::cout << "Choosing option >>>>>>>>>>\t test_beam_2023 \t<<<<<<<<<<" <= 7357 && run_number <= 12792) { + std::cout << "Choosing option >>>>>>>>>>\t ti18_2024_2025 \t<<<<<<<<<<" <= 4361 && run_number <= 7356){ + std::cout << "Choosing option >>>>>>>>>>\t ti18_2022_2023 \t<<<<<<<<<<" < +#include +#include +#include + +#include "MuFilter.h" +#include "MuFilterHit.h" +#include "ShipUnit.h" +#include "FairLogger.h" + +snd::analysis_tools::DSPlane::DSPlane(std::vector snd_hits, const Configuration &configuration, MuFilter *muon_filter_geometry, int station) : configuration_(configuration), station_(station) +{ + for ( auto mu_hit : snd_hits) + { + TVector3 A, B; + int detectorID = mu_hit->GetDetectorID(); + muon_filter_geometry->GetPosition(detectorID, A, B); + const int n_sipms = muon_filter_geometry->GetnSiPMs(detectorID); + const int n_sides = muon_filter_geometry->GetnSides(detectorID); + for (int i{0}; i < n_sipms * n_sides; ++i) + { + if (mu_hit->isMasked(i) || mu_hit->GetSignal(i) < -990.) continue; + DSHit hit; + hit.bar = static_cast(detectorID % 1000); + hit.channel_index = n_sipms * n_sides * hit.bar + i; + hit.is_right = (i >= n_sipms) ? true : false; + hit.timestamp = configuration_.is_mc ? mu_hit->GetTime(i) / ShipUnit::snd_TDC2ns : mu_hit->GetTime(i); + hit.qdc = mu_hit->GetSignal(i); + // use the left and right measurements to calculate the x coordinate along the bar + if (!mu_hit->isVertical()) { + hit.is_x = false; + float tmp_x = mu_hit->GetImpactXpos(true, true, false, configuration_.is_mc); + hit.x = (tmp_x < -990.) ? std::nan("") : A.X() - tmp_x; + hit.y = A.Y(); + } + else { + hit.is_x = true; + hit.x = A.X(); + hit.y = std::nan(""); + } + hit.z = A.Z(); + hits_.push_back(hit); + } + } +} + +const double snd::analysis_tools::DSPlane::GetTotQdc() const +{ + double tot_qdc = std::accumulate(hits_.begin(), hits_.end(), 0.0, + [](double sum, const auto &b) {return sum + b.qdc;}); + return tot_qdc; +} + +const int snd::analysis_tools::DSPlane::GetBarNHits(int bar_to_compute) const +{ + int bar_hit = std::count_if(hits_.begin(), hits_.end(), + [bar_to_compute](const auto &hit) {return hit.bar == bar_to_compute;}); + return bar_hit; +} + +void snd::analysis_tools::DSPlane::TimeFilter(double min_timestamp, double max_timestamp) +{ + hits_.erase(std::remove_if(hits_.begin(), hits_.end(), + [&](auto &hit) + { return hit.timestamp < min_timestamp || hit.timestamp > max_timestamp; }), + hits_.end()); +} + +const int snd::analysis_tools::DSPlane::GetNHitBars() const +{ + int count{0}; + for (int bar{0}; bar < 2*configuration_.ds_bar_per_station; ++bar) { + if (GetBarNHits(bar) > 0) count++; + } + return count; +} + +void snd::analysis_tools::DSPlane::DSHit::Print() const +{ + LOGF(INFO, "DSHit ch_idx :%d\tposition: (%f,%f,%f)\ttime: %f\tqdc: %f\tbar: %d\tis_x: %d\tis_right: %d", channel_index, x, y, z, timestamp, qdc, bar, is_x, is_right); +} \ No newline at end of file diff --git a/analysis/tools/sndDSPlane.h b/analysis/tools/sndDSPlane.h new file mode 100644 index 0000000000..f658f577e0 --- /dev/null +++ b/analysis/tools/sndDSPlane.h @@ -0,0 +1,62 @@ +#ifndef SND_DSPLANE_H +#define SND_DSPLANE_H + +#include + +#include "MuFilter.h" +#include "MuFilterHit.h" +#include "sndConfiguration.h" + +namespace snd { + namespace analysis_tools { + class DSPlane + { + public: + + // right and left side of DS bars + template + struct rl_pair + { + T right{}; + T left{}; + }; + + // hits vector, each hit has info about timestamp, qdc and position + struct DSHit + { + int channel_index; + int bar; + + double qdc; + double timestamp; + double x; + double y; + double z; + + bool is_x; // true if vertical (measures x) + bool is_right; + + void Print() const; + }; + + DSPlane(std::vector snd_hits, const Configuration &configuration, MuFilter *muon_filter_geometry, int station); + + const int GetNHits() const { return hits_.size(); }; + const int GetStation() const { return station_; } + + const double GetTotQdc() const; + const std::vector GetHits() const { return hits_; }; + const int GetBarNHits(int bar_to_compute) const; + + void TimeFilter(double min_timestamp, double max_timestamp); + const int GetNHitBars() const; + + private: + std::vector hits_; + Configuration configuration_; + int station_; + }; + } +} + +#endif diff --git a/analysis/tools/sndGeometryGetter.cxx b/analysis/tools/sndGeometryGetter.cxx index 5764f9f98d..964e710392 100644 --- a/analysis/tools/sndGeometryGetter.cxx +++ b/analysis/tools/sndGeometryGetter.cxx @@ -3,6 +3,7 @@ #include #include #include +#include #include "Scifi.h" #include "MuFilter.h" @@ -11,8 +12,11 @@ // Get geometry full path, works with test beam too // 2022 constants are included in the 2023 geofile -std::string snd::analysis_tools::GetGeoPath(const std::string& csv_file_path, int run_number) +std::string snd::analysis_tools::GetGeoPath(int run_number, std::string csv_file_path) { + if (csv_file_path.empty()) { + csv_file_path = std::string(getenv("SNDSW_ROOT")) + "/analysis/tools/geo_paths.csv"; + } std::ifstream file(csv_file_path); if (!file.is_open()) { throw std::runtime_error("Could not open CSV file: " + csv_file_path); @@ -42,7 +46,7 @@ std::string snd::analysis_tools::GetGeoPath(const std::string& csv_file_path, in } // Get SciFi and MuFilter geometries -std::pair snd::analysis_tools::GetGeometry(std::string geometry_path) +std::pair snd::analysis_tools::GetGeometry(const std::string& geometry_path) { TPython::Exec("import SndlhcGeo"); TPython::Exec(("SndlhcGeo.GeoInterface('" + geometry_path + "')").c_str()); @@ -59,9 +63,9 @@ std::pair snd::analysis_tools::GetGeometry(std::string geom } // Get SciFi and MuFilter geometries directly from run number -std::pair snd::analysis_tools::GetGeometry(const std::string& csv_file_path, int run_number) +std::pair snd::analysis_tools::GetGeometry(int run_number, const std::string& csv_file_path) { - std::string geometry_path = GetGeoPath(csv_file_path, run_number); + std::string geometry_path = GetGeoPath(run_number, csv_file_path); return GetGeometry(geometry_path); } diff --git a/analysis/tools/sndGeometryGetter.h b/analysis/tools/sndGeometryGetter.h index 59d5e9e1b7..8b53d276a6 100644 --- a/analysis/tools/sndGeometryGetter.h +++ b/analysis/tools/sndGeometryGetter.h @@ -9,9 +9,10 @@ namespace snd { namespace analysis_tools { - std::string GetGeoPath(const std::string& csv_file_path, int run_number); - std::pair GetGeometry(std::string geometry_path); - std::pair GetGeometry(const std::string& csv_file_path, int run_number); + // With empty csv_file_path, it gets the official one from the sndsw installation + std::string GetGeoPath(int run_number, std::string csv_file_path = ""); + std::pair GetGeometry(const std::string& geometry_path); + std::pair GetGeometry(int run_number, const std::string& csv_file_path = ""); } } diff --git a/analysis/tools/sndPlaneTools.cxx b/analysis/tools/sndPlaneTools.cxx new file mode 100644 index 0000000000..0dbbe0b0ff --- /dev/null +++ b/analysis/tools/sndPlaneTools.cxx @@ -0,0 +1,112 @@ +#include "sndPlaneTools.h" + +#include +#include + +#include "TClonesArray.h" +#include "Scifi.h" +#include "MuFilter.h" +#include "sndConfiguration.h" +#include "sndVetoPlane.h" +#include "sndScifiPlane.h" +#include "sndUSPlane.h" +#include "sndDSPlane.h" +#include "sndScifiHit.h" +#include "MuFilterHit.h" + +std::vector snd::analysis_tools::FillVeto(const snd::Configuration &configuration, TClonesArray *mufi_hits, MuFilter *mufilter_geometry) +{ + + std::vector veto_planes; + const int n_mufi_hits{mufi_hits->GetEntries()}; + + const int n_station = configuration.veto_n_stations; + std::vector> plane_hits(n_station); + + for (int i{0}; i < n_mufi_hits; ++i) { + auto hit = static_cast(mufi_hits->At(i)); + if (hit->GetSystem()!=1) continue; + int station_id = hit->GetPlane(); + if (station_id > -1 && station_id < n_station) { + plane_hits[station_id].push_back(hit); + } + else throw std::runtime_error{"Invalid Veto plane"}; + } + for (int st{0}; st < n_station; ++st) { + veto_planes.emplace_back(snd::analysis_tools::VetoPlane(plane_hits[st], configuration, mufilter_geometry, st+1)); + } + return veto_planes; +} + +std::vector snd::analysis_tools::FillScifi(const snd::Configuration &configuration, TClonesArray *sf_hits, Scifi *scifi_geometry) +{ + + std::vector scifi_planes; + const int n_sf_hits{sf_hits->GetEntries()}; + + const int max_station = configuration.scifi_n_stations; + std::vector> stations_hits(max_station); + + for (int i{0}; i < n_sf_hits; ++i) { + auto hit = static_cast(sf_hits->At(i)); + int station_id = hit->GetStation()-1; + + if (station_id > -1 && station_id < max_station) { + stations_hits[station_id].push_back(hit); + } + else throw std::runtime_error{"Invalid SciFi station"}; + } + for (int st{0}; st < max_station; ++st) { + scifi_planes.emplace_back(snd::analysis_tools::ScifiPlane(stations_hits[st], configuration, scifi_geometry, st+1)); + } + return scifi_planes; +} + + +std::vector snd::analysis_tools::FillUS(const snd::Configuration &configuration, TClonesArray *mufi_hits, MuFilter *mufilter_geometry, bool use_small_sipms) +{ + + std::vector us_planes; + const int n_mufi_hits{mufi_hits->GetEntries()}; + + const int n_station = configuration.us_n_stations; + std::vector> plane_hits(n_station); + + for (int i{0}; i < n_mufi_hits; ++i) { + auto hit = static_cast(mufi_hits->At(i)); + if (hit->GetSystem()!=2) continue; + int station_id = hit->GetPlane(); + if (station_id > -1 && station_id < n_station) { + plane_hits[station_id].push_back(hit); + } + else throw std::runtime_error{"Invalid US plane"}; + } + for (int st{0}; st < n_station; ++st) { + us_planes.emplace_back(snd::analysis_tools::USPlane(plane_hits[st], configuration, mufilter_geometry, st+1, use_small_sipms)); + } + return us_planes; +} + +std::vector snd::analysis_tools::FillDS(const snd::Configuration &configuration, TClonesArray *mufi_hits, MuFilter *mufilter_geometry) +{ + + std::vector ds_planes; + const int n_mufi_hits{mufi_hits->GetEntries()}; + + const int n_station = configuration.ds_n_stations; + std::vector> plane_hits(n_station); + + for (int i{0}; i < n_mufi_hits; ++i) { + auto hit = static_cast(mufi_hits->At(i)); + if (hit->GetSystem()!=3) continue; + int station_id = hit->GetPlane(); + if (station_id > -1 && station_id < n_station) { + plane_hits[station_id].push_back(hit); + } + else throw std::runtime_error{"Invalid DS plane"}; + } + for (int st{0}; st < n_station; ++st) { + ds_planes.emplace_back(snd::analysis_tools::DSPlane(plane_hits[st], configuration, mufilter_geometry, st+1)); + } + return ds_planes; +} diff --git a/analysis/tools/sndPlaneTools.h b/analysis/tools/sndPlaneTools.h new file mode 100644 index 0000000000..25e3ac7f78 --- /dev/null +++ b/analysis/tools/sndPlaneTools.h @@ -0,0 +1,25 @@ +#ifndef SND_PLANETOOLS_H +#define SND_PLANETOOLS_H + +#include + +#include "TClonesArray.h" +#include "Scifi.h" +#include "MuFilter.h" +#include "sndConfiguration.h" +#include "sndVetoPlane.h" +#include "sndScifiPlane.h" +#include "sndUSPlane.h" +#include "sndDSPlane.h" + +namespace snd { + namespace analysis_tools { + // Produce veto, scifi, us and ds planes from data + std::vector FillVeto(const Configuration &configuration, TClonesArray *mufi_hits, MuFilter *mufilter_geometry); + std::vector FillScifi(const Configuration &configuration, TClonesArray *sf_hits, Scifi *scifi_geometry); + std::vector FillUS(const Configuration &configuration, TClonesArray *mufi_hits, MuFilter *mufilter_geometry, bool use_small_sipms=false); + std::vector FillDS(const Configuration &configuration, TClonesArray *mufi_hits, MuFilter *mufilter_geometry); + } +} + +#endif diff --git a/analysis/tools/sndSciFiTools.cxx b/analysis/tools/sndSciFiTools.cxx index 8405243c90..c8ee2a5ec6 100644 --- a/analysis/tools/sndSciFiTools.cxx +++ b/analysis/tools/sndSciFiTools.cxx @@ -9,6 +9,8 @@ #include "sndScifiHit.h" #include "ROOT/TSeq.hxx" #include "Scifi.h" +#include "TROOT.h" +#include "ROOT/RRangeCast.hxx" void snd::analysis_tools::getSciFiHitsPerStation(const TClonesArray *digiHits, std::vector &horizontal_hits, std::vector &vertical_hits) @@ -22,7 +24,7 @@ void snd::analysis_tools::getSciFiHitsPerStation(const TClonesArray *digiHits, s sndScifiHit *hit; TIter hitIterator(digiHits); - while (hit = dynamic_cast(hitIterator.Next())) { + while ((hit = dynamic_cast(hitIterator.Next()))) { if (hit->isValid()) { int station = hit->GetStation(); if (hit->isVertical()) { @@ -132,8 +134,10 @@ float snd::analysis_tools::peakScifiTiming(const TClonesArray &digiHits, int bin TH1F ScifiTiming("Timing", "Scifi Timing", bins, min_x, max_x); - int refStation = ((sndScifiHit *)digiHits.At(0))->GetStation(); - bool refOrientation = ((sndScifiHit *)digiHits.At(0))->isVertical(); + Scifi *ScifiDet = dynamic_cast (gROOT->GetListOfGlobals()->FindObject("Scifi") ); + auto* first_hit = static_cast(digiHits[0]); + int refStation = first_hit->GetStation(); + bool refOrientation = first_hit->isVertical(); float hitTime = -1.0; float timeConversion = 1.; if (!isMC) { @@ -146,6 +150,10 @@ float snd::analysis_tools::peakScifiTiming(const TClonesArray &digiHits, int bin continue; } hitTime = hit->GetTime() * timeConversion; + if (!isMC){ + int id_hit = hit ->GetDetectorID(); + hitTime = ScifiDet->GetCorrectedTime(id_hit, hitTime, 0); + } if (hitTime < min_x || hitTime > max_x) { continue; } @@ -311,7 +319,6 @@ snd::analysis_tools::filterScifiHits(const TClonesArray &digiHits, LOG(info) << "\"TI18\" setup will be used by default, please provide \"H8\" for the Testbeam setup."; } - sndScifiHit *hit; TIter hitIterator(&supportArray); if (method == 0) { @@ -330,8 +337,8 @@ snd::analysis_tools::filterScifiHits(const TClonesArray &digiHits, for (auto station : ROOT::MakeSeq(1, ScifiStations + 1)) { for (auto orientation : {false, true}) { - auto supportArray = selectScifiHits(digiHits, station, orientation, selection_parameters, true, isMC); - for (auto *p : *supportArray) { + auto supportArray_p = selectScifiHits(digiHits, station, orientation, selection_parameters, true, isMC); + for (auto *p : *supportArray_p) { auto *hit = dynamic_cast(p); if (hit->isValid()) { new ((*filteredHits)[filteredHitsIndex++]) sndScifiHit(*hit); @@ -568,11 +575,39 @@ snd::analysis_tools::findCentreOfGravityPerStation(const TClonesArray* digiHits, if (!digiHits) { LOG(ERROR) << "Error: digiHits is null in findCentreOfGravityPerStation"; } - std::vector x_positions; - std::vector y_positions; - TVector3 A, B; - for (auto* obj : *digiHits) { - auto* hit = dynamic_cast(obj); + + auto [x_positions, y_positions] = snd::analysis_tools::hitPositionVectorsPerStation(digiHits, station, ScifiDet); + + if (x_positions.empty()) { + LOG(ERROR) << "Error: No hits enter."; + } + double meanX = computeMean(x_positions); + if (y_positions.empty()) { + LOG(ERROR) << "Error: No hits enter."; + } + double meanY = computeMean(y_positions); + return {meanX, meanY}; +} + +std::pair,std::vector> snd::analysis_tools::hitPositionVectorsPerStation(const TClonesArray *digiHits, int station, Scifi* ScifiDet){ + /*This function returns the hit vector position for the digiHits in both orientations + Arguments: + digiHits: A TClonesArray containing the hits in the SciFi detector. + station: The station number for which to retrieve the hit positions. + ScifiDet: A pointer to the Scifi detector object to retrieve SiPM positions. + Returns: + A pair of vectors containing the x and y positions of the hits in the specified station. + If no hits are found, it returns empty vectors and logs an error message. + */ + + if (!digiHits) { + LOG(ERROR) << "Error: digiHits is null"; + } + std::vector x_positions; + std::vector y_positions; + + TVector3 A, B; + for (auto* hit : ROOT::RRangeCast(*digiHits)) { if (!hit || !hit->isValid()) { continue; } @@ -586,14 +621,109 @@ snd::analysis_tools::findCentreOfGravityPerStation(const TClonesArray* digiHits, else { y_positions.push_back((A.Y() + B.Y()) * 0.5); } - } - if (x_positions.empty()) { - LOG(ERROR) << "Error: No hits enter."; - } - double meanX = computeMean(x_positions); - if (y_positions.empty()) { - LOG(ERROR) << "Error: No hits enter."; - } - double meanY = computeMean(y_positions); - return {meanX, meanY}; -} \ No newline at end of file + } + if (x_positions.empty()) { + LOG(ERROR) << "Error: No hits enter."; + } + if (y_positions.empty()) { + LOG(ERROR) << "Error: No hits enter."; + } + return {x_positions, y_positions}; +} + + +double hitWeightComputation(std::vector hit_position) { + /* This function returns the summation of the hit weights + where a hit weight is the number of neighbouring hits within 1 cm position (default width). + + Arguments: + hit_position: A vector containing the positions of the hits in a specific SciFi station. + + Returns: + The sum of the weights of the hits in the vector. + Returns 0 if the vector is empty or if the sum of weights is 0. + */ + + if (hit_position.empty()) { + LOG(INFO) << "Warning: The hit position vector is empty." << std::endl; + return 0; // Return 0 if the vector is empty + } + + int n_hits = hit_position.size(); + double width = 1.0; // Width around the hit to consider as a neighbour, in cm + std::vector weights; // Vector to store the weights of each hit + + // Sorting the hits for efficient neighbour counting + std::sort(hit_position.begin(), hit_position.end()); + + // Initialize pointers for the sliding window. + // Both pointers will only move forward (or stay put) across the outer loop iterations. + int left_ptr = 0; + int right_ptr = 0; + + // Loop through each hit position to calculate its neighbour count + for (int i = 0; i < n_hits; i++) { + double specific_hit = hit_position[i]; + + // Move right_ptr forward to include all hits within (specific_hit + width). + // It will point to the first element *just outside* the right boundary of the window. + while (right_ptr < n_hits && hit_position[right_ptr] <= specific_hit + width) { + right_ptr++; + } + + // Move left_ptr forward to exclude all hits *less than* (specific_hit - width). + // It will point to the first element *just inside* or at the left boundary of the window. + while (left_ptr < n_hits && hit_position[left_ptr] < specific_hit - width) { + left_ptr++; + } + + // Calculate neighbouring hits within the window: + // The number of hits in the current window is (right_ptr - left_ptr). + // This count includes the 'specific_hit' itself. + double count_in_window = right_ptr - left_ptr; + + // Subtract 1 to exclude the 'specific_hit' itself from the neighbour count. + // We must ensure that 'specific_hit' (hit_position[i]) is actually within the window + // defined by left_ptr and right_ptr. Since the array is sorted and we are iterating + // through 'i', hit_position[i] will always be between hit_position[left_ptr] and + // hit_position[right_ptr-1] (if left_ptr <= i < right_ptr). + double neighbour_no_of_hits = count_in_window - 1; + + // Ensure the neighbour count is non-negative + if (neighbour_no_of_hits < 0) { + neighbour_no_of_hits = 0; + } + + weights.push_back(neighbour_no_of_hits); + } + + // Calculate the sum of all computed weights + double sum_weights = std::accumulate(weights.begin(), weights.end(), 0.0); + + return sum_weights; +} + +std::pair snd::analysis_tools::hitDensityPerStation(const TClonesArray *digiHits, int station, Scifi* ScifiDet){ + /* + This function returns the hit density which is described as the summation of the hit weights + where a hit weigth is the number of neighbouring hits within 1 cm postion (default width) + arguments: hit position vector : vector ontaining the positions of the hits in a specific SciFi station + returns: the sum of the weights of the hits in the vector + If the vector is empty, it returns 0. + If the sum of weights is 0, it returns 0. + */ + std::pair, std::vector> hit_position_vec = hitPositionVectorsPerStation(digiHits, station, ScifiDet); + std::vector hit_position_x = hit_position_vec.first; // Assuming we are interested in X positions + std::vector hit_position_y = hit_position_vec.second; // Assuming we are interested in Y positions + + if (hit_position_x.size()==0 && hit_position_y.size()==0) { + LOG(INFO)<< "Warning: The hit position vector is empty." << std::endl; + return {0.0,0.0}; // Check if the vector is empty + } + + double sum_weights_x = hitWeightComputation(hit_position_x); + double sum_weights_y = hitWeightComputation(hit_position_y); + + return {sum_weights_x, sum_weights_y}; +} + diff --git a/analysis/tools/sndSciFiTools.h b/analysis/tools/sndSciFiTools.h index f8967dfe33..fec3ded4e4 100644 --- a/analysis/tools/sndSciFiTools.h +++ b/analysis/tools/sndSciFiTools.h @@ -39,7 +39,7 @@ namespace snd { // selection_parameters[3] // make_selection determines whether you want to select the hits for the respective plane // or not and skip that step (in case the hits are already curated) - std::unique_ptr selectScifiHits(const TClonesArray &digiHits, int station, bool orientation, int bins_x=52, float min_x=0.0, float max_x=26.0, float time_lower_range=1E9/(2*ShipUnit::snd_freq), float time_upper_range=1.2E9/(ShipUnit::snd_freq/ShipUnit::hertz), bool make_selection=true, bool isMC=false); + std::unique_ptr selectScifiHits(const TClonesArray &digiHits, int station, bool orientation, int bins_x=52, float min_x=0.0, float max_x=26.0, float time_lower_range=1E9/(2*ShipUnit::snd_freq/ShipUnit::hertz), float time_upper_range=1.2E9/(ShipUnit::snd_freq/ShipUnit::hertz), bool make_selection=true, bool isMC=false); std::unique_ptr selectScifiHits(const TClonesArray &digiHits, int station, bool orientation, const std::map &selection_parameters, bool make_selection=true, bool isMC=false); // Function to obtain the ScifiHits that are considered to be useful from an event @@ -79,5 +79,13 @@ namespace snd { // Find the Center of Particle Showering on the SciFi plane std::pair findCentreOfGravityPerStation(const TClonesArray* digiHits, int station, Scifi* ScifiDet); + + // Function to get the hit position vectors for a specific station + // Returns a pair of vectors containing the x and y positions of the hits in the specified station + std::pair, std::vector> hitPositionVectorsPerStation(const TClonesArray* digiHits, int station, Scifi* ScifiDet); + + // Function to compute the hit density per station + // Returns a pair of doubles containing the hit density for horizontal and vertical orientations of that station + std::pair hitDensityPerStation(const TClonesArray* digiHits, int station, Scifi* ScifiDet); } } diff --git a/analysis/tools/sndScifiPlane.cxx b/analysis/tools/sndScifiPlane.cxx new file mode 100644 index 0000000000..311eb22609 --- /dev/null +++ b/analysis/tools/sndScifiPlane.cxx @@ -0,0 +1,302 @@ +#include "sndScifiPlane.h" + +#include +#include +#include +#include +#include + +#include "TVector3.h" +#include "Scifi.h" +#include "sndScifiHit.h" +#include "ShipUnit.h" +#include "Math/Point3D.h" +#include "FairLogger.h" + +snd::analysis_tools::ScifiPlane::ScifiPlane(std::vector snd_hits, const Configuration &configuration, Scifi *scifi_geometry, int station) : configuration_(configuration), centroid_(std::nan(""), std::nan(""), std::nan("")), centroid_error_(std::nan(""), std::nan(""), std::nan("")), station_(station) +{ + for ( auto snd_hit : snd_hits) + { + if (!snd_hit->isValid()) continue; + ScifiHit hit; + hit.channel_index = 512 * snd_hit->GetMat() + 128 * snd_hit->GetSiPM() + snd_hit->GetSiPMChan(); + int detectorID = snd_hit->GetDetectorID(); + if (configuration_.is_mc) { + hit.timestamp = snd_hit->GetTime(0) / ShipUnit::snd_TDC2ns; // timestamp is in clock cycles + } + else { + hit.timestamp = (scifi_geometry->GetCorrectedTime(detectorID, snd_hit->GetTime(0)*ShipUnit::snd_TDC2ns, 0) / ShipUnit::snd_TDC2ns); // timestamp is in clock cycles + } + hit.qdc = snd_hit->GetSignal(0); + hit.is_x = snd_hit->isVertical(); + + TVector3 A, B; + scifi_geometry->GetSiPMPosition(detectorID, A, B); + hit.z = A.Z(); + if (hit.is_x) + { + hit.x = A.X(); + hit.y = std::nan(""); + } + else + { + hit.x = std::nan(""); + hit.y = A.Y(); + } + hits_.push_back(hit); + } +} + +const snd::analysis_tools::ScifiPlane::xy_pair snd::analysis_tools::ScifiPlane::GetNHits() const +{ + xy_pair counts{0, 0}; + counts.x = std::count_if(hits_.begin(), hits_.end(), [](auto &hit) + { return hit.is_x; }); + counts.y = static_cast(hits_.size()) - counts.x; + + return counts; +} + +bool snd::analysis_tools::ScifiPlane::HasShower() const +{ + if (configuration_.scifi_min_hits_in_window > configuration_.scifi_shower_window_width) + { + throw std::runtime_error{"min_hits > window_width"}; + } + + if (configuration_.scifi_shower_window_width > configuration_.scifi_n_channels_per_plane) + { + throw std::runtime_error{"window_width > n_channels_per_plane"}; + } + + xy_pair> is_hit; + is_hit.x.resize(configuration_.scifi_n_channels_per_plane, 0); + is_hit.y.resize(configuration_.scifi_n_channels_per_plane, 0); + + for (auto &hit : hits_) + { + if (hit.channel_index < 0 || + hit.channel_index >= configuration_.scifi_n_channels_per_plane) + { + throw std::out_of_range{"hit.channel_index out of range"}; + } + (hit.is_x ? is_hit.x : is_hit.y)[hit.channel_index] = 1; + } + + auto density = [&](const std::vector &hit_arr) + { + int count{0}; + + // Initial count for the first window + for (int i{0}; i < configuration_.scifi_shower_window_width; ++i) + { + count += hit_arr[i]; + } + + if (count >= configuration_.scifi_min_hits_in_window) + return true; + + // Slide the window across the array + for (int i = configuration_.scifi_shower_window_width; i < configuration_.scifi_n_channels_per_plane; ++i) + { + count += hit_arr[i] - hit_arr[i - configuration_.scifi_shower_window_width]; + if (count >= configuration_.scifi_min_hits_in_window) + return true; + } + + return false; + }; + + return density(is_hit.x) && density(is_hit.y); +} + +const ROOT::Math::XYZPoint snd::analysis_tools::ScifiPlane::GetCluster(int max_gap) const +{ + std::vector pos_x(configuration_.scifi_n_channels_per_plane, std::nan("")); + std::vector pos_y(configuration_.scifi_n_channels_per_plane, std::nan("")); + + for (auto &hit : hits_) + { + if (hit.qdc > 0) + { + if (hit.is_x) + { + pos_x[hit.channel_index] = hit.x; + } + else + { + pos_y[hit.channel_index] = hit.y; + } + } + } + + auto largest_cluster = [&](const std::vector &positions) + { + int n = static_cast(positions.size()); + + int best_start = -1, best_end = -1, best_size = 0; + int start = -1, gap_count = 0, size = 0; + + // Find the largest cluster + for (int i = 0; i < n; ++i) + { + if (!std::isnan(positions[i])) + { + if (start == -1) + start = i; // Start a new cluster + size++; + gap_count = 0; // Reset consecutive gap counter + } + else + { + gap_count++; + if (gap_count > max_gap) + { // Too many consecutive gaps, finalize previous cluster + if (size > best_size) + { + best_start = start; + best_end = i - gap_count; // End before the excessive gaps + best_size = size; + } + // Reset for a new potential cluster + start = -1; + gap_count = 0; + size = 0; + } + else + { + size++; // Gaps within max_gap still count in cluster size + } + } + } + + // Check last cluster + if (size > best_size) + { + best_start = start; + best_end = n - 1; + best_size = size; + } + + // Compute the average of non-gap values in the best cluster + if (best_start == -1 || best_end == -1) + return std::nan(""); // No valid cluster found + + double sum = 0.0; + int count = 0; + for (int i = best_start; i <= best_end; ++i) + { + if (!std::isnan(positions[i])) + { + sum += positions[i]; + count++; + } + } + + return (count > 0) ? sum / count : std::nan(""); + }; + + double cluster_x = largest_cluster(pos_x); + double cluster_y = largest_cluster(pos_y); + if (!(std::isnan(cluster_x) || std::isnan(cluster_y))) { + return ROOT::Math::XYZPoint(cluster_x, cluster_y, hits_[0].z); + } + return ROOT::Math::XYZPoint(std::nan(""), std::nan(""), std::nan("")); +} + +void snd::analysis_tools::ScifiPlane::TimeFilter(double min_timestamp, double max_timestamp) +{ + hits_.erase(std::remove_if(hits_.begin(), hits_.end(), + [&](auto &hit) + { return hit.timestamp < min_timestamp || hit.timestamp > max_timestamp; }), + hits_.end()); +} + +snd::analysis_tools::ScifiPlane::xy_pair snd::analysis_tools::ScifiPlane::GetPointQdc(const ROOT::Math::XYZPoint &point, double radius) const +{ + xy_pair qdc{0.0, 0.0}; + for (const auto &hit : hits_) { + if (hit.is_x && std::abs(hit.x - point.X()) <= radius) { + qdc.x += hit.qdc; + } + else if (!hit.is_x && std::abs(hit.y - point.Y()) <= radius) { + qdc.y += hit.qdc; + } + } + return qdc; +} + +void snd::analysis_tools::ScifiPlane::FindCentroid() +{ + centroid_.SetXYZ(0, 0, 0); + double tot_qdc_x{0}, sum_qdc2_x{0.0}; + double tot_qdc_y{0}, sum_qdc2_y{0.0}; + std::vector cleaned_hits = hits_; + + cleaned_hits.erase(std::remove_if(cleaned_hits.begin(), cleaned_hits.end(), + [&](auto &hit) + { return hit.qdc <= 0; }), + cleaned_hits.end()); + int counts_x = std::count_if(cleaned_hits.begin(), cleaned_hits.end(), [](auto &hit) + { return hit.is_x; }); + int counts_y = static_cast(cleaned_hits.size())-counts_x; + if (counts_x < configuration_.scifi_min_n_hits_for_centroid && counts_y < configuration_.scifi_min_n_hits_for_centroid ) { + centroid_ = ROOT::Math::XYZPoint(std::nan(""), std::nan(""), std::nan("")); + return; + } + + for (auto &hit : cleaned_hits) + { + if (hit.qdc > 0) + { + if (hit.is_x) + { + centroid_.SetX(centroid_.X() + hit.x * hit.qdc); + tot_qdc_x += hit.qdc; + sum_qdc2_x += hit.qdc*hit.qdc; + } + else + { + centroid_.SetY(centroid_.Y() + hit.y * hit.qdc); + tot_qdc_y += hit.qdc; + sum_qdc2_y += hit.qdc*hit.qdc; + } + + centroid_.SetZ(centroid_.Z() + hit.z * hit.qdc); + } + } + centroid_ = ROOT::Math::XYZPoint((tot_qdc_x > 0) ? centroid_.X() / tot_qdc_x : std::nan(""), (tot_qdc_y > 0) ? centroid_.Y() / tot_qdc_y : std::nan(""), (tot_qdc_x+tot_qdc_y > 0) ? centroid_.Z() / (tot_qdc_x+tot_qdc_y) : std::nan("")); + auto qdc_x_error_scaler = sqrt(sum_qdc2_x)/tot_qdc_x; + auto qdc_y_error_scaler = sqrt(sum_qdc2_y)/tot_qdc_y; + auto qdc_z_error_scaler = sqrt(sum_qdc2_x+sum_qdc2_y)/(tot_qdc_x+tot_qdc_y); + centroid_error_ = ROOT::Math::XYZPoint(configuration_.scifi_centroid_error_x*qdc_x_error_scaler, + configuration_.scifi_centroid_error_y*qdc_y_error_scaler, + configuration_.scifi_centroid_error_z*qdc_z_error_scaler); +} + +const snd::analysis_tools::ScifiPlane::xy_pair snd::analysis_tools::ScifiPlane::GetTotQdc(bool only_positive) const +{ + xy_pair qdc_sum{0.0, 0.0}; + qdc_sum.x = std::accumulate(hits_.begin(), hits_.end(), 0.0, [&](double current_sum, auto &hit) + { return (hit.is_x && (hit.qdc > 0 || !only_positive)) ? (current_sum + hit.qdc) : current_sum; }); + qdc_sum.y = std::accumulate(hits_.begin(), hits_.end(), 0.0, [&](double current_sum, auto &hit) + { return (!hit.is_x && (hit.qdc > 0 || !only_positive)) ? (current_sum + hit.qdc) : current_sum; }); + + return qdc_sum; +} + +const snd::analysis_tools::ScifiPlane::xy_pair snd::analysis_tools::ScifiPlane::GetTotEnergy(bool only_positive) const +{ + xy_pair energy{0.0, 0.0}; + auto qdc = GetTotQdc(only_positive); + + energy.x = qdc.x * configuration_.scifi_qdc_to_gev; + energy.y = qdc.y * configuration_.scifi_qdc_to_gev; + + return energy; +} + +void snd::analysis_tools::ScifiPlane::ScifiHit::Print() const +{ + LOGF(INFO, "ScifiHit ch_idx :%d\tposition: (%f,%f,%f)\ttime: %f\tqdc: %f\tis_x: %d", channel_index, x, y, z, timestamp, qdc, is_x); +} \ No newline at end of file diff --git a/analysis/tools/sndScifiPlane.h b/analysis/tools/sndScifiPlane.h new file mode 100644 index 0000000000..2867344eee --- /dev/null +++ b/analysis/tools/sndScifiPlane.h @@ -0,0 +1,66 @@ +#ifndef SND_SCIFIPLANE_H +#define SND_SCIFIPLANE_H + +#include + +#include "Scifi.h" +#include "sndConfiguration.h" +#include "sndScifiHit.h" +#include "Math/Point3D.h" + +namespace snd { + namespace analysis_tools { + class ScifiPlane + { + public: + // x and y planes + template + struct xy_pair + { + T x{}; + T y{}; + }; + + struct ScifiHit + { + double qdc{}; + double timestamp{}; // timestamp is in clock cycles + double x{}; + double y{}; + double z{}; + int channel_index{}; + bool is_x{}; + + void Print() const; + }; + + ScifiPlane(std::vector snd_hits, const Configuration &configuration, Scifi *scifi_geometry, int station); + + const int GetStation() const { return station_; } + const std::vector GetHits() const { return hits_; } + const ROOT::Math::XYZPoint GetCentroid() const { return centroid_; } + const ROOT::Math::XYZPoint GetCentroidError() const { return centroid_error_; } + const xy_pair GetTotQdc(bool only_positive = false) const; + const xy_pair GetTotEnergy(bool only_positive = false) const; + const xy_pair GetNHits() const; + // Position of larger cluster of consecutive hits, allowing at most max_gap channels with no hit + const ROOT::Math::XYZPoint GetCluster(int max_gap) const; + // The centroid is the qdc-weighted mean of hit positions, considering only hits with positive qdc + void FindCentroid(); + bool HasShower() const; + void TimeFilter(double min_timestamp, double max_timestamp); + // qdc from hits within a given point and radius (square, not circle) + xy_pair GetPointQdc(const ROOT::Math::XYZPoint &point, double radius) const; + + private: + std::vector hits_; + Configuration configuration_; + ROOT::Math::XYZPoint centroid_; + ROOT::Math::XYZPoint centroid_error_; + + int station_; + }; + } +} + +#endif \ No newline at end of file diff --git a/analysis/tools/sndShowerTools.cxx b/analysis/tools/sndShowerTools.cxx new file mode 100644 index 0000000000..9b064b76f1 --- /dev/null +++ b/analysis/tools/sndShowerTools.cxx @@ -0,0 +1,279 @@ +#include "sndShowerTools.h" + +#include +#include +#include + +#include "sndConfiguration.h" +#include "sndScifiPlane.h" +#include "sndUSPlane.h" +#include "TMatrixD.h" +#include "TMatrixDSym.h" +#include "TMatrixDSymEigen.h" +#include "TVectorD.h" +#include "TDecompChol.h" +#include "Math/Vector3D.h" +#include "Math/Point3D.h" + +int snd::analysis_tools::GetScifiShowerStart(const std::vector &scifi_planes) +{ + // Assuming planes are ordered + for (const auto &p : scifi_planes) + { + if (p.HasShower()) + { + return p.GetStation(); + } + } + return -1; +} + +int snd::analysis_tools::GetScifiShowerEnd(const std::vector &scifi_planes) +{ + // Assuming planes are ordered + for (auto p = scifi_planes.rbegin(); p != scifi_planes.rend(); p++) { + if (p->HasShower()) + { + return p->GetStation(); + } + } + return -1; +} + +int snd::analysis_tools::GetUSShowerStart(const std::vector &us_planes) +{ + // Assuming planes are ordered + for (const auto &p : us_planes) + { + if (p.HasShower()) + { + return p.GetStation(); + } + } + return -1; +} + +int snd::analysis_tools::GetUSShowerEnd(const std::vector &us_planes) +{ + // Assuming planes are ordered + for (auto p = us_planes.rbegin(); p != us_planes.rend(); p++) { + if (p->HasShower()) + { + return p->GetStation(); + } + } + return -1; +} + +std::pair snd::analysis_tools::GetShowerInterceptAndDirection(const snd::Configuration &configuration, const std::vector &scifi_planes, const std::vector &us_planes) +{ + std::vector positions_x; + std::vector positions_y; + std::vector positions_z; + std::vector err_positions_x; + std::vector err_positions_y; + std::vector err_positions_z; + + auto collect_centroids = [&](const auto &planes) { + for (const auto &p : planes) { + auto centroid = p.GetCentroid(); + auto centroid_error = p.GetCentroidError(); + // Check that centroid is valid + if (!(std::isnan(centroid.X()) || std::isnan(centroid.Y()) || std::isnan(centroid.Z()))) { + positions_x.push_back(centroid.X()); + positions_y.push_back(centroid.Y()); + positions_z.push_back(centroid.Z()); + err_positions_x.push_back(centroid_error.X()); + err_positions_y.push_back(centroid_error.Y()); + err_positions_z.push_back(centroid_error.Z()); + } + } + }; + + collect_centroids(scifi_planes); + collect_centroids(us_planes); + + if (static_cast(positions_z.size()) < configuration.centroid_min_valid_station) { + return std::make_pair(ROOT::Math::XYZPoint(std::nan(""), std::nan(""), std::nan("")), ROOT::Math::XYZVector(std::nan(""), std::nan(""), std::nan(""))); + } + + auto weighted_linear_fit = [](const std::vector &z, const std::vector &val, const std::vector &err) + { + const int n_points = z.size(); + const int n_variables = 2; // intercept + slope + + // Wrap existing memory into TVectorD + TVectorD v_z; + v_z.Use(n_points, const_cast(z.data())); + TVectorD v_y; + v_y.Use(n_points, const_cast(val.data())); + TVectorD v_e; + v_e.Use(n_points, const_cast(err.data())); + + TMatrixD matrix(n_points, n_variables); + TMatrixDColumn(matrix,0) = 1.0; + TMatrixDColumn(matrix,1) = v_z; + + // Solve normal equations (weighted least squares) + TVectorD coeffs = NormalEqn(matrix, v_y, v_e); + + return std::make_pair(coeffs[0], coeffs[1]); // (intercept, slope) + }; + + auto [intercept_x, slope_x] = weighted_linear_fit(positions_z, positions_x, err_positions_x); + auto [intercept_y, slope_y] = weighted_linear_fit(positions_z, positions_y, err_positions_y); + + ROOT::Math::XYZPoint shower_intercept(intercept_x, intercept_y, 0.0); + ROOT::Math::XYZVector shower_direction(slope_x, slope_y, 1.0); + + return std::make_pair(shower_intercept, shower_direction.Unit()); +} + +std::pair, std::vector> snd::analysis_tools::GetShoweringPlanes(const std::vector &scifi_planes, const std::vector &us_planes) +{ + std::vector sh_scifi_planes; + std::vector sh_us_planes; + + // Filter showering Scifi planes + std::copy_if(scifi_planes.begin(), scifi_planes.end(), std::back_inserter(sh_scifi_planes), [](const snd::analysis_tools::ScifiPlane& p) { return p.HasShower(); }); + // Filter showering US planes + std::copy_if(us_planes.begin(), us_planes.end(), std::back_inserter(sh_us_planes), [](const snd::analysis_tools::USPlane& p) { return p.HasShower(); }); + + for (auto &p: sh_scifi_planes) { + p.FindCentroid(); + } + for (auto &p: sh_us_planes) { + p.FindCentroid(); + } + + return std::make_pair(sh_scifi_planes, sh_us_planes); +} + +std::pairPCA(const std::vector& u, const std::vector& z) +{ + // Need at least 3 measurements + assert(u.size() == z.size() && u.size() >= 2); + const size_t n_samples = u.size(); + const int n_variables = 2; + + // Find the mean values + double mean_u = 0.0, mean_z = 0.0; + for (size_t i = 0; i < n_samples; ++i) + { + mean_u += u[i]; + mean_z += z[i]; + } + mean_u /= n_samples; + mean_z /= n_samples; + + // Find the covariance components + double cov_uu = 0.0, cov_zz = 0.0, cov_uz = 0.0; + for (size_t i = 0; i < n_samples; ++i) + { + double du = u[i] - mean_u; + double dz = z[i] - mean_z; + cov_uu += du * du; + cov_zz += dz * dz; + cov_uz += du * dz; + } + + const double denom = (n_samples > 1) ? (n_samples - 1.0) : 1.0; + cov_uu /= denom; + cov_zz /= denom; + cov_uz /= denom; + + // Fill the covariance matrix + TMatrixDSym cov_matrix(n_variables); + cov_matrix[0][0] = cov_uu; + cov_matrix[0][1] = cov_uz; + cov_matrix[1][0] = cov_uz; + cov_matrix[1][1] = cov_zz; + + // Find the eigenvalues and eigenvectors. For now we only need the former. + // Perhaps the eigenvectors could be used in the future, so left that part in comments. + TMatrixDSymEigen eig_decomp(cov_matrix); + TVectorD evals = eig_decomp.GetEigenValues(); + //TMatrixD evecs = eig_decomp.GetEigenVectors(); + + // Make sure the order is descending, swap if needed + if (evals[0] < evals[1]) + { + std::swap(evals[0], evals[1]); + //evecs[0] = eig_decomp.GetEigenVectors()[1]; + //evecs[1] = eig_decomp.GetEigenVectors()[0]; + } + + return {evals[0], evals[1]}; +} + +std::pair snd::analysis_tools::GetSciFiSpatialAnisotropy(const std::vector &scifi_planes, bool use_all_centroids) +{ + std::vector x, y, zx, zy; + for (auto &p : scifi_planes) { + // Add measurements per projection and in the respective coordinate vectors + if (use_all_centroids == false) { + for (auto &hit : p.GetHits()) { + if (hit.is_x) { + x.push_back(hit.x); + zx.push_back(hit.z); + } else { + y.push_back(hit.y); + zy.push_back(hit.z); + } + } + } else { + auto centroid = p.GetCentroid(); + // Check that centroid is valid + if (!std::isnan(centroid.Z())) { + if (!std::isnan(centroid.X())) { + x.push_back(centroid.X()); + zx.push_back(centroid.Z()); + } + if (!std::isnan(centroid.Y())) { + y.push_back(centroid.Y()); + zy.push_back(centroid.Z()); + } + } + } + } + + // Need at least 3 measurements + if (x.size() < 2 || y.size() < 2 || zx.size() < 2 || zy.size() < 2) { + return {1., 1.}; + } + double lambda1, lambda2; + std::tie(lambda1, lambda2) = PCA(x, zx); + double anisotropy_xz = (lambda1 > 0) ? 1.0 - lambda2 / lambda1 : 0.0; + std::tie(lambda1, lambda2) = PCA(y, zy); + double anisotropy_yz = (lambda1 > 0) ? 1.0 - lambda2 / lambda1 : 0.0; + return {anisotropy_xz, anisotropy_yz}; +} + +double snd::analysis_tools::GetUSSpatialAnisotropy(const std::vector &us_planes, bool use_all_centroids) +{ + std::vector y, zy; + for (auto &p : us_planes) { + // Add measurements per projection and in the respective coordinate vectors + if (use_all_centroids == false) { + for (auto &hit : p.GetHits()) { + y.push_back(hit.y); + zy.push_back(hit.z); + } + } else { + auto centroid = p.GetCentroid(); + // Check that centroid is valid + if (!(std::isnan(centroid.Z()) || std::isnan(centroid.Y()))) { + y.push_back(centroid.Y()); + zy.push_back(centroid.Z()); + } + } + } + + // Need at least 3 measurements + if (y.size() < 2 || zy.size() < 2) { + return 1.; + } + double lambda1, lambda2; + std::tie(lambda1, lambda2) = PCA(y, zy); + return (lambda1 > 0) ? 1.0 - lambda2 / lambda1 : 0.0; +} diff --git a/analysis/tools/sndShowerTools.h b/analysis/tools/sndShowerTools.h new file mode 100644 index 0000000000..79faaeb3bc --- /dev/null +++ b/analysis/tools/sndShowerTools.h @@ -0,0 +1,33 @@ +#ifndef SND_SHOWERTOOLS_H +#define SND_SHOWERTOOLS_H + +#include + +#include "sndConfiguration.h" +#include "sndScifiPlane.h" +#include "sndUSPlane.h" +#include "Math/Vector3D.h" +#include "Math/Point3D.h" + +namespace snd { + namespace analysis_tools { + // Returns first SciFi station with shower. If no shower is found in SciFi, returns -1 + int GetScifiShowerStart(const std::vector &scifi_planes); + // Returns last SciFi station with shower. If no shower is found in SciFi, returns -1 + int GetScifiShowerEnd(const std::vector &scifi_planes); + // Returns first US station with shower. If no shower is found in US, returns -1 + int GetUSShowerStart(const std::vector &us_planes); + // Returns last US station with shower. If no shower is found in US, returns -1 + int GetUSShowerEnd(const std::vector &us_planes); + // Returns a 3D reference point and the normalized direction of the shower fitting centroids in SciFi and US + std::pair GetShowerInterceptAndDirection(const Configuration &configuration, const std::vector &scifi_planes, const std::vector &us_planes); + // Filters out non showering planes + std::pair, std::vector> GetShoweringPlanes(const std::vector &scifi_planes, const std::vector &us_planes); + // Returns the SciFi spatial anisotropy(track-likeness) per projection + std::pair GetSciFiSpatialAnisotropy(const std::vector &scifi_planes, bool use_all_centroids = false); + // Returns the US spatial anisotropy + double GetUSSpatialAnisotropy(const std::vector &us_planes, bool use_all_centroids = false); + } +} + +#endif diff --git a/analysis/tools/sndTchainGetter.cxx b/analysis/tools/sndTchainGetter.cxx index 6e65267a61..fd64e88c1a 100644 --- a/analysis/tools/sndTchainGetter.cxx +++ b/analysis/tools/sndTchainGetter.cxx @@ -5,11 +5,15 @@ #include #include #include +#include #include "TChain.h" #include -std::string snd::analysis_tools::GetDataBasePath(const std::string& csv_file_path, int run_number) { +std::string snd::analysis_tools::GetDataBasePath(int run_number, std::string csv_file_path) { + if (csv_file_path.empty()) { + csv_file_path = std::string(getenv("SNDSW_ROOT")) + "/analysis/tools/data_paths.csv"; + } std::ifstream file(csv_file_path); if (!file.is_open()) { throw std::runtime_error("Could not open CSV file: " + csv_file_path); @@ -38,8 +42,8 @@ std::string snd::analysis_tools::GetDataBasePath(const std::string& csv_file_pat throw std::runtime_error("Run number not found in CSV mapping."); } -std::unique_ptr snd::analysis_tools::GetTChain(const std::string& csv_file_path, int run_number, int n_files){ - std::string base_folder = GetDataBasePath(csv_file_path, run_number); +std::unique_ptr snd::analysis_tools::GetTChain(int run_number, int n_files, const std::string& csv_file_path){ + std::string base_folder = GetDataBasePath(run_number, csv_file_path); auto tchain = std::make_unique("rawConv"); if (n_files == -1) { tchain->Add(Form("%srun_%06d/sndsw_raw-*", base_folder.c_str(), run_number)); @@ -52,7 +56,7 @@ std::unique_ptr snd::analysis_tools::GetTChain(const std::string& csv_fi return tchain; }; -std::unique_ptr snd::analysis_tools::GetTChain(std::string file_name){ +std::unique_ptr snd::analysis_tools::GetTChain(const std::string& file_name){ auto tchain = std::make_unique("rawConv"); tchain->Add(file_name.c_str()); return tchain; diff --git a/analysis/tools/sndTchainGetter.h b/analysis/tools/sndTchainGetter.h index ca4a0423d2..4013689725 100644 --- a/analysis/tools/sndTchainGetter.h +++ b/analysis/tools/sndTchainGetter.h @@ -8,9 +8,11 @@ namespace snd { namespace analysis_tools { - std::unique_ptr GetTChain(const std::string& csv_file_path, int run_number, int n_files = -1); - std::unique_ptr GetTChain(std::string file_name); - std::string GetDataBasePath(const std::string& csv_file_path, int run_number); + // With n_files == -1 all run partitions are added to the chain. + // With empty csv_file_path, it gets the official one from the sndsw installation + std::unique_ptr GetTChain(int run_number, int n_files = -1, const std::string& csv_file_path = ""); + std::unique_ptr GetTChain(const std::string& file_name); + std::string GetDataBasePath(int run_number, std::string csv_file_path = ""); } } diff --git a/analysis/tools/sndUSPlane.cxx b/analysis/tools/sndUSPlane.cxx new file mode 100644 index 0000000000..b3e29a091a --- /dev/null +++ b/analysis/tools/sndUSPlane.cxx @@ -0,0 +1,202 @@ +#include "sndUSPlane.h" + +#include +#include +#include +#include + +#include "TVector3.h" +#include "MuFilter.h" +#include "MuFilterHit.h" +#include "ShipUnit.h" + +snd::analysis_tools::USPlane::USPlane(std::vector snd_hits, const Configuration &configuration, MuFilter *muon_filter_geometry, int station, bool use_small_sipms) : configuration_(configuration), centroid_(std::nan(""), std::nan(""), std::nan("")), centroid_error_(std::nan(""), std::nan(""),std::nan("")), station_(station) +{ + for ( auto mu_hit : snd_hits) + { + TVector3 A, B; + int detectorID = mu_hit->GetDetectorID(); + muon_filter_geometry->GetPosition(detectorID, A, B); + for (int i{0}; i < 16; ++i) + { + if (mu_hit->isMasked(i) || mu_hit->GetSignal(i) < -990.) continue; + USHit hit; + hit.bar = static_cast(detectorID % 1000); + hit.channel_index = 16 * hit.bar + i; + hit.is_large = !mu_hit->isShort(i); + hit.is_right = i > 7 ? true : false; + + if (!hit.is_large && !use_small_sipms) + { + hit.timestamp = std::nan(""); + hit.qdc = std::nan(""); + hit.x = std::nan(""); + hit.y = std::nan(""); + hit.z = std::nan(""); + } + else + { + hit.timestamp = configuration_.is_mc ? mu_hit->GetTime(i) / ShipUnit::snd_TDC2ns : mu_hit->GetTime(i); + hit.qdc = mu_hit->GetSignal(i); + // use the left and right measurements to calculate the x coordinate along the bar + float tmp_x = mu_hit->GetImpactXpos(true, true, false, configuration_.is_mc); + hit.x = (tmp_x < -990.) ? std::nan("") : A.X() - tmp_x; + hit.y = A.Y(); + hit.z = A.Z(); + } + hits_.push_back(hit); + } + } +} + +void snd::analysis_tools::USPlane::FindCentroid() +{ + // select large SiPM channels with positive qdc + std::vector cleaned_hits = hits_; + + cleaned_hits.erase(std::remove_if(cleaned_hits.begin(), cleaned_hits.end(), + [&](auto &hit) + { return (hit.qdc <= 0 || hit.is_large == false || std::isnan(hit.x)); }), + cleaned_hits.end()); + + // min number of hit in the plane to attempt to find a centroid + if (static_cast(cleaned_hits.size()) < configuration_.us_min_n_hits_for_centroid) + { + return; + } + + // weigthed sum calculated per plane + double weighted_sum_x{0.0}, weighted_sum_y{0.0}, weighted_sum_z{0.0}; + double total_qdc_positive{0.0}, sum_qdc2_positive{0.0}; + // loop over hits in the plane + for (const auto &hit : cleaned_hits) + { + weighted_sum_x += hit.x * hit.qdc; + weighted_sum_y += hit.y * hit.qdc; + weighted_sum_z += hit.z * hit.qdc; + total_qdc_positive += hit.qdc; + sum_qdc2_positive += hit.qdc*hit.qdc; + } + weighted_sum_x /= total_qdc_positive; + weighted_sum_y /= total_qdc_positive; + weighted_sum_z /= total_qdc_positive; + centroid_.SetXYZ(weighted_sum_x, weighted_sum_y, weighted_sum_z); + auto qdc_error_scaler = sqrt(sum_qdc2_positive)/total_qdc_positive; + centroid_error_.SetXYZ(configuration_.us_centroid_error_x*qdc_error_scaler, + configuration_.us_centroid_error_y*qdc_error_scaler, + configuration_.us_centroid_error_z*qdc_error_scaler); +} + +const snd::analysis_tools::USPlane::sl_pair snd::analysis_tools::USPlane::GetTotQdc() const +{ + sl_pair totQdc{0.0, 0.0}; + for (const auto &hit : hits_) + { + if (hit.is_large) + totQdc.large += hit.qdc; + else + totQdc.small += hit.qdc; + } + return totQdc; +} + +const snd::analysis_tools::USPlane::sl_pair snd::analysis_tools::USPlane::GetTotEnergy() const +{ + sl_pair tot_energy{0.0, 0.0}; + sl_pair tot_qdc = GetTotQdc(); + + tot_energy.large = tot_qdc.large *configuration_.us_qdc_to_gev; + tot_energy.small = tot_qdc.small *configuration_.us_qdc_to_gev; + + return tot_energy; +} + + +const snd::analysis_tools::USPlane::rl_pair snd::analysis_tools::USPlane::GetSideQdc() const +{ + rl_pair side_qdc{0.0, 0.0}; + for (const auto &hit : hits_) + { + if (hit.is_large) + { + if (hit.is_right) + side_qdc.right += hit.qdc; + else + side_qdc.left += hit.qdc; + } + } + return side_qdc; +} + +const snd::analysis_tools::USPlane::rl_pair snd::analysis_tools::USPlane::GetBarQdc(int bar_to_compute) const +{ + rl_pair bar_qdc{0.0, 0.0}; + for (const auto &hit : hits_) + { + if (hit.bar != bar_to_compute) + continue; + else + { + if (hit.is_large) + { + if (hit.is_right) + bar_qdc.right += hit.qdc; + else + bar_qdc.left += hit.qdc; + } + } + } + return bar_qdc; +} + +const snd::analysis_tools::USPlane::sl_pair snd::analysis_tools::USPlane::GetBarNHits(int bar_to_compute) const +{ + sl_pair bar_hit{0, 0}; + for (const auto &hit : hits_) + { + if (hit.bar != bar_to_compute) + continue; + else + { + if (hit.is_large) + + bar_hit.large++; + else + bar_hit.small++; + + } + } + return bar_hit; +} + +void snd::analysis_tools::USPlane::TimeFilter(double min_timestamp, double max_timestamp) +{ + hits_.erase(std::remove_if(hits_.begin(), hits_.end(), + [&](auto &hit) + { return hit.timestamp < min_timestamp || hit.timestamp > max_timestamp; }), + hits_.end()); +} + +const snd::analysis_tools::USPlane::sl_pair snd::analysis_tools::USPlane::GetNHits() const +{ + sl_pair counts{0, 0}; + counts.large = std::count_if(hits_.begin(), hits_.end(), [](auto &hit) + { return hit.is_large; }); + counts.small = static_cast(hits_.size()) - counts.large; + + return counts; +} + +const int snd::analysis_tools::USPlane::GetNHitBars() const +{ + int count{0}; + for (int bar{0}; bar < configuration_.us_bar_per_station; ++bar) { + if (GetBarNHits(bar).large > configuration_.us_min_hit_on_bar) count++; + } + return count; +} + +void snd::analysis_tools::USPlane::USHit::Print() const +{ + LOGF(INFO, "USHit ch_idx :%d\tposition: (%f,%f,%f)\ttime: %f\tqdc: %f\tbar: %d\tis_right: %d\tis_large: %d", channel_index, x, y, z, timestamp, qdc, bar, is_right, is_large); +} \ No newline at end of file diff --git a/analysis/tools/sndUSPlane.h b/analysis/tools/sndUSPlane.h new file mode 100644 index 0000000000..12261666e3 --- /dev/null +++ b/analysis/tools/sndUSPlane.h @@ -0,0 +1,81 @@ +#ifndef SND_USPLANE_H +#define SND_USPLANE_H + +#include + +#include "Math/Point3D.h" +#include "MuFilter.h" +#include "MuFilterHit.h" +#include "sndConfiguration.h" + +namespace snd { + namespace analysis_tools { + class USPlane + { + public: + // small and large SiPMs + template + struct sl_pair + { + T small{}; + T large{}; + }; + + // right and left side of US bars + template + struct rl_pair + { + T right{}; + T left{}; + }; + + // hits vector, each hit has info about timestamp, qdc and position + struct USHit + { + int channel_index; + int bar; + + double qdc; + double timestamp; + double x; + double y; + double z; + + bool is_large; + bool is_right; + + void Print() const; + }; + + USPlane(std::vector snd_hits, const Configuration &configuration, MuFilter *muon_filter_geometry, int station, bool use_small_sipms=false); + + const sl_pair GetNHits() const; + const int GetStation() const { return station_; } + + const sl_pair GetTotQdc() const; + const sl_pair GetTotEnergy() const; + const rl_pair GetSideQdc() const; + const rl_pair GetBarQdc(int bar_to_compute) const; + const sl_pair GetBarNHits(int bar_to_compute) const; + const std::vector GetHits() const { return hits_; }; + double HasShower() const { return GetNHits().large >= configuration_.us_min_n_hits_for_centroid; } + // The centroid is the qdc-weighted mean of hit positions, considering only hits with positive qdc + void FindCentroid(); + ROOT::Math::XYZPoint GetCentroid() const { return centroid_; } + ROOT::Math::XYZPoint GetCentroidError() const { return centroid_error_; } + const int GetNHitBars() const; + + + void TimeFilter(double min_timestamp, double max_timestamp); + + private: + std::vector hits_; + Configuration configuration_; + ROOT::Math::XYZPoint centroid_; + ROOT::Math::XYZPoint centroid_error_; + int station_; + }; + } +} + +#endif diff --git a/analysis/tools/sndVetoPlane.cxx b/analysis/tools/sndVetoPlane.cxx new file mode 100644 index 0000000000..ba0509e705 --- /dev/null +++ b/analysis/tools/sndVetoPlane.cxx @@ -0,0 +1,90 @@ +#include "sndVetoPlane.h" + +#include +#include +#include +#include + +#include "MuFilter.h" +#include "MuFilterHit.h" +#include "ShipUnit.h" +#include "FairLogger.h" + +snd::analysis_tools::VetoPlane::VetoPlane(std::vector snd_hits, const Configuration &configuration, MuFilter *muon_filter_geometry, int station) : configuration_(configuration), station_(station) +{ + for ( auto mu_hit : snd_hits) + { + TVector3 A, B; + int detectorID = mu_hit->GetDetectorID(); + muon_filter_geometry->GetPosition(detectorID, A, B); + const int n_sipms = muon_filter_geometry->GetnSiPMs(detectorID); + const int n_sides = muon_filter_geometry->GetnSides(detectorID); + for (int i{0}; i < n_sipms * n_sides; ++i) + { + if (mu_hit->isMasked(i) || mu_hit->GetSignal(i) < -990.) continue; + VetoHit hit; + hit.bar = static_cast(detectorID % 1000); + hit.channel_index = n_sipms * n_sides * hit.bar + i; + hit.is_right = (i >= n_sipms) ? true : false; + hit.timestamp = configuration_.is_mc ? mu_hit->GetTime(i) / ShipUnit::snd_TDC2ns : mu_hit->GetTime(i); + hit.qdc = mu_hit->GetSignal(i); + // use the left and right measurements to calculate the x coordinate along the bar + if (!mu_hit->isVertical()) { + hit.is_x = false; + float tmp_x = mu_hit->GetImpactXpos(true, true, false, configuration_.is_mc); + hit.x = (tmp_x < -990.) ? std::nan("") : A.X() - tmp_x; + hit.y = A.Y(); + } + else { + hit.is_x = true; + hit.x = A.X(); + hit.y = std::nan(""); + } + hit.z = A.Z(); + hits_.push_back(hit); + } + } +} + +const double snd::analysis_tools::VetoPlane::GetTotQdc() const +{ + double tot_qdc = std::accumulate(hits_.begin(), hits_.end(), 0.0, + [](double sum, const auto &b) {return sum + b.qdc;}); + return tot_qdc; +} + +const double snd::analysis_tools::VetoPlane::GetBarQdc(int bar_to_compute) const +{ + double bar_qdc = std::accumulate(hits_.begin(), hits_.end(), 0.0, + [bar_to_compute](double sum, const auto &b) {return (b.bar == bar_to_compute) ? sum + b.qdc : sum;}); + return bar_qdc; +} + +const int snd::analysis_tools::VetoPlane::GetBarNHits(int bar_to_compute) const +{ + int bar_hit = std::count_if(hits_.begin(), hits_.end(), + [bar_to_compute](const auto &hit) {return hit.bar == bar_to_compute;}); + return bar_hit; +} + +void snd::analysis_tools::VetoPlane::TimeFilter(double min_timestamp, double max_timestamp) +{ + hits_.erase(std::remove_if(hits_.begin(), hits_.end(), + [&](auto &hit) + { return hit.timestamp < min_timestamp || hit.timestamp > max_timestamp; }), + hits_.end()); +} + +const int snd::analysis_tools::VetoPlane::GetNHitBars() const +{ + int count{0}; + for (int bar{0}; bar < configuration_.veto_bar_per_station; ++bar) { + if (GetBarNHits(bar) > configuration_.veto_min_hit_on_bar) count++; + } + return count; +} + +void snd::analysis_tools::VetoPlane::VetoHit::Print() const +{ + LOGF(INFO, "VetoHit ch_idx :%d\tposition: (%f,%f,%f)\ttime: %f\tqdc: %f\tbar: %d\tis_x: %d\tis_right: %d", channel_index, x, y, z, timestamp, qdc, bar, is_x, is_right); +} \ No newline at end of file diff --git a/analysis/tools/sndVetoPlane.h b/analysis/tools/sndVetoPlane.h new file mode 100644 index 0000000000..eaa765dec7 --- /dev/null +++ b/analysis/tools/sndVetoPlane.h @@ -0,0 +1,63 @@ +#ifndef SND_VETOPLANE_H +#define SND_VETOPLANE_H + +#include + +#include "MuFilter.h" +#include "MuFilterHit.h" +#include "sndConfiguration.h" + +namespace snd { + namespace analysis_tools { + class VetoPlane + { + public: + + // right and left side of Veto bars + template + struct rl_pair + { + T right{}; + T left{}; + }; + + // hits vector, each hit has info about timestamp, qdc and position + struct VetoHit + { + int channel_index; + int bar; + + double qdc; + double timestamp; + double x; + double y; + double z; + + bool is_x; // true if vertical (measures x) + bool is_right; + + void Print() const; + }; + + VetoPlane(std::vector snd_hits, const Configuration &configuration, MuFilter *muon_filter_geometry, int station); + + const int GetNHits() const { return hits_.size(); }; + const int GetStation() const { return station_; } + + const double GetTotQdc() const; + const double GetBarQdc(int bar_to_compute) const; + const int GetBarNHits(int bar_to_compute) const; + const std::vector GetHits() const { return hits_; }; + const int GetNHitBars() const; + + void TimeFilter(double min_timestamp, double max_timestamp); + + private: + std::vector hits_; + Configuration configuration_; + int station_; + }; + } +} + +#endif diff --git a/gconfig/DecayConfigPy8.C b/gconfig/DecayConfigPy8.C index 3eba51e693..634b506b72 100644 --- a/gconfig/DecayConfigPy8.C +++ b/gconfig/DecayConfigPy8.C @@ -1,40 +1,59 @@ -void DecayConfig() { - - // This script uses the external decayer TPythia8Decayer in place of the - // concrete Monte Carlo native decay mechanisms only for the - // specific types of decays defined below. - - // Access the external decayer singleton and initialize it - TVirtualMCDecayer* decayer = new TPythia8Decayer(); - decayer->Init(); - - // Tell the concrete monte carlo to use the external decayer. The - // external decayer will be used for: - // i)particle decays not defined in concrete monte carlo, or - //ii)particles for which the concrete monte carlo is told - // to use the external decayer for its type via: - // gMC->SetUserDecay(pdgId); - // If this is invoked, the external decayer will be used for particles - // of type pdgId even if the concrete monte carlo has a decay mode - // already defined for that particle type. - gMC->SetExternalDecayer(decayer); - // to get the rare muon decays, HOWEVER does not work with Geant4 logic - //gMC->SetUserDecay(221); // eta - //gMC->SetUserDecay(223); // omega - //gMC->SetUserDecay(113); // rho0 - //gMC->SetUserDecay(331); // eta_prime - //gMC->SetUserDecay(333); // phi - gMC->SetUserDecay(411); - gMC->SetUserDecay(-411); - gMC->SetUserDecay(421); - gMC->SetUserDecay(-421); - gMC->SetUserDecay(4122); - gMC->SetUserDecay(-4122); - gMC->SetUserDecay(431); - gMC->SetUserDecay(-431); - gMC->SetUserDecay(15); - gMC->SetUserDecay(-15); - cout<< "External decayer DecayConfigPy8 initialized"<Init(); // Currently does nothing, but leave it in case it gets implemented. + + TVirtualMC::GetMC()->SetExternalDecayer(decayer); + + // Set list of particles to decay with external decayer via geant4 command + std::string g4cmd = "/mcPhysics/setExtDecayerSelection"; + std::cout << "DecayConfigPy8 macro: setting heavy flavour decays via Pythia8" << std::endl; + // First loop through all particles in database and add heavy flavour + for (TObject * obj : *(db->ParticleList())) { + TParticlePDG * p = dynamic_cast(obj); + + if (!p) { + std::cout << "DecayConfigPy8 WARNING: got object other than TParticlePDG from particle list" << std::endl; + continue; + } + + bool isCharm = (strcmp("CharmedMeson", p->ParticleClass()) == 0) or + (strcmp("CharmedBaryon", p->ParticleClass()) == 0) or + p->Charm(); // p->Charm() doesn't work. Leave here for future compat. + bool isBeauty = (strcmp("B-Meson", p->ParticleClass()) == 0) or + (strcmp("B-Baryon", p->ParticleClass()) == 0) or + p->Beauty(); // p->Beauty() doesn't work. Leave here for future compat. + + if (isCharm or isBeauty){ + g4cmd += " "; + g4cmd += p->GetName(); + } + } + + std::cout << "DecayConfigPy8 macro: setting other decays (short-lived and tau) via Pythia8" << std::endl; + std::vector pdgs_to_decay_with_pythia8 = { + // Short-lived resonances, for rare decays. + 221, // eta + // 223, // omega CRASHES + // 113, // rho CRASHES + 331, // eta_prime + // 333, // phi CRASHES + // Tau + 15, // tau + -15 // tau bar + }; + + for (int pdg : pdgs_to_decay_with_pythia8){ + g4cmd += " "; + g4cmd += db->GetParticle(pdg)->GetName(); + } + dynamic_cast(TVirtualMC::GetMC())->ProcessGeantCommand(g4cmd.c_str()); +} +void DecayConfig() { DecayConfigPythia8(); } diff --git a/gconfig/SetCuts.C b/gconfig/SetCuts.C index d59b6ebc24..8b97e2ccc8 100755 --- a/gconfig/SetCuts.C +++ b/gconfig/SetCuts.C @@ -1,6 +1,8 @@ /** Configuration macro for setting common cuts and processes for G3, G4 and Fluka (M. Al-Turany 27.03.2008) specific cuts and processes to g3 or g4 should be set in the g3Config.C, g4Config.C or flConfig.C + Comment out all cuts defined in this macro (S. Ilieva 3.03.2026) + essentially removing the 1 MeV particle production cut-off for SND@LHC G4 simulations */ @@ -33,20 +35,22 @@ void SetCuts() gMC->SetProcess("LOSS",1); /**energy loss*/ gMC->SetProcess("MULS",1); /**multiple scattering*/ - Double_t cut1 = 1.0E-3; // GeV --> 1 MeV - Double_t cutb = 1.0E4; // GeV --> 10 TeV - Double_t tofmax = 1.E10; // seconds - cout << "SetCuts Macro: Setting cuts.." <SetCut("CUTGAM",cut1); /** gammas (GeV)*/ - gMC->SetCut("CUTELE",cut1); /** electrons (GeV)*/ - gMC->SetCut("CUTNEU",cut1); /** neutral hadrons (GeV)*/ - gMC->SetCut("CUTHAD",cut1); /** charged hadrons (GeV)*/ - gMC->SetCut("CUTMUO",cut1); /** muons (GeV)*/ - gMC->SetCut("BCUTE",cut1); /** electron bremsstrahlung (GeV)*/ - gMC->SetCut("BCUTM",cut1); /** muon and hadron bremsstrahlung(GeV)*/ - gMC->SetCut("DCUTE",cut1); /** delta-rays by electrons (GeV)*/ - gMC->SetCut("DCUTM",cut1); /** delta-rays by muons (GeV)*/ - gMC->SetCut("PPCUTM",cut1); /** direct pair production by muons (GeV)*/ - gMC->SetCut("TOFMAX",tofmax); /**time of flight cut in seconds*/ + cout << "SetCuts Macro: No cuts to be set in here." < 1 MeV + //Double_t cutb = 1.0E4; // GeV --> 10 TeV + //Double_t tofmax = 1.E10; // seconds + //cout << "SetCuts Macro: Setting cuts.." <SetCut("CUTGAM",cut1); /** gammas (GeV)*/ + //gMC->SetCut("CUTELE",cut1); /** electrons (GeV)*/ + //gMC->SetCut("CUTNEU",cut1); /** neutral hadrons (GeV)*/ + //gMC->SetCut("CUTHAD",cut1); /** charged hadrons (GeV)*/ + //gMC->SetCut("CUTMUO",cut1); /** muons (GeV)*/ + //gMC->SetCut("BCUTE",cut1); /** electron bremsstrahlung (GeV)*/ + //gMC->SetCut("BCUTM",cut1); /** muon and hadron bremsstrahlung(GeV)*/ + //gMC->SetCut("DCUTE",cut1); /** delta-rays by electrons (GeV)*/ + //gMC->SetCut("DCUTM",cut1); /** delta-rays by muons (GeV)*/ + //gMC->SetCut("PPCUTM",cut1); /** direct pair production by muons (GeV)*/ + //gMC->SetCut("TOFMAX",tofmax); /**time of flight cut in seconds*/ } diff --git a/gconfig/g4Config.C b/gconfig/g4Config.C index d698dd305f..912404603d 100755 --- a/gconfig/g4Config.C +++ b/gconfig/g4Config.C @@ -25,7 +25,7 @@ void Config() /// When more than one options are selected, they should be separated with '+' /// character: eg. stepLimit+specialCuts. TG4RunConfiguration* runConfiguration - = new TG4RunConfiguration("geomRoot", "QGSP_BERT_HP_PEN", "stepLimiter+specialCuts+specialControls"); + = new TG4RunConfiguration("geomRoot", "FTFP_BERT_HP_EMZ", "stepLimiter+specialCuts+specialControls"); /// Create the G4 VMC TGeant4* geant4 = new TGeant4("TGeant4", "The Geant4 Monte Carlo", runConfiguration); diff --git a/gconfig/g4config.in b/gconfig/g4config.in index de79f5aeea..5f918a17eb 100644 --- a/gconfig/g4config.in +++ b/gconfig/g4config.in @@ -4,8 +4,13 @@ /mcVerbose/all 0 /mcTracking/loopVerbose 0 # suggested by Ivana Hrivnacova to switch off *** Particle reached max step number .... -/mcPhysics/g4NeutronHPVerbose 0 # - suppress the warnings about missing neutron data -/mcPhysics/g4HadronicProcessStoreVerbose 0 + +# Geant4 range cut for all geo regions and all supported particle types: e+,e-, gamma, proton +/mcPhysics/rangeCuts 2 mm +# range cut for all geo regions and selected particle e.g. +#/mcPhysics/rangeCutForGamma 2 mm +# print the Geant4 and VMC cuts per material +#/mcRegions/print true # switch on other sources of di muon /physics_lists/em/PositronToMuons true diff --git a/geometry/sndLHC_H4geom_config.py b/geometry/sndLHC_H4geom_config.py index e7349b5a6a..ebc877d4ca 100644 --- a/geometry/sndLHC_H4geom_config.py +++ b/geometry/sndLHC_H4geom_config.py @@ -259,16 +259,17 @@ c.MuFilter.DS4ShiftX = 1.39 * u.cm #digitization parameters - c.MuFilter.DsAttenuationLength = 350 * u.cm # values between 300 cm and 400cm observed for H6 testbeam - c.MuFilter.DsTAttenuationLength = 700 * u.cm # top readout with mirror on bottom - c.MuFilter.VandUpAttenuationLength = 999 * u.cm # no significante attenuation observed for H6 testbeam - c.MuFilter.VandUpSiPMcalibrationL = 25.*1000. # 1.65 MeV = 41 qcd - c.MuFilter.VandUpSiPMcalibrationS = 25.*1000. + c.MuFilter.DsAttenuationLength = 230*u.cm # values between 130cm and 330cm are observed for TI18 in years 2022-2025, but between 300 cm and 400cm for H6 testbeam + c.MuFilter.DsTAttenuationLength = 700*u.cm # top readout with mirror on bottom, TI18 and H6 observables agree + c.MuFilter.VandUpAttenuationLength = 210*u.cm # while no significant attenuation observed for H6 testbeam + c.MuFilter.VTAttenuationLength = 999*u.cm # Veto 3, no significant attenuation observed in 2024-2025 data + c.MuFilter.VandUpSiPMcalibrationL = 50.*1000. # 1.65 MeV = 41 qcd over 6 Large SiPMs(one side) + c.MuFilter.VandUpSiPMcalibrationS = 0. # no MIP signal for small SiPMs, delayed and compromised response in general c.MuFilter.DsSiPMcalibration = 25.*1000. c.MuFilter.timeResol = 150.*u.picosecond - c.MuFilter.VandUpPropSpeed = 12.5*u.cm/u.nanosecond - c.MuFilter.DsPropSpeed = 14.3*u.cm/u.nanosecond - + c.MuFilter.VandUpPropSpeed = 13.6*u.cm/u.nanosecond + c.MuFilter.DsPropSpeed = 15.1*u.cm/u.nanosecond + c.Floor = AttrDict(z=48000.*u.cm) # to place tunnel in SND_@LHC coordinate system c.Floor.DX = 1.0*u.cm c.Floor.DY = -4.5*u.cm # subtract 4.5cm to avoid overlaps diff --git a/geometry/sndLHC_H6geom_config.py b/geometry/sndLHC_H6geom_config.py index 0951ef5a75..17560ef099 100644 --- a/geometry/sndLHC_H6geom_config.py +++ b/geometry/sndLHC_H6geom_config.py @@ -245,15 +245,16 @@ c.MuFilter.VETOBoxY2 = c.MuFilter.VETOLocY + c.MuFilter.VetoBarZ/2 + c.MuFilter.SupportBoxD #digitization parameters - c.MuFilter.DsAttenuationLength = 350 * u.cm # values between 300 cm and 400cm observed for H6 testbeam - c.MuFilter.DsTAttenuationLength = 700 * u.cm # top readout with mirror on bottom - c.MuFilter.VandUpAttenuationLength = 999 * u.cm # no significante attenuation observed for H6 testbeam - c.MuFilter.VandUpSiPMcalibrationL = 25.*1000. # 1.65 MeV = 41 qcd - c.MuFilter.VandUpSiPMcalibrationS = 25.*1000. + c.MuFilter.DsAttenuationLength = 230*u.cm # values between 130cm and 330cm are observed for TI18 in years 2022-2025, but between 300 cm and 400cm for H6 testbeam + c.MuFilter.DsTAttenuationLength = 700*u.cm # top readout with mirror on bottom, TI18 and H6 observables agree + c.MuFilter.VandUpAttenuationLength = 210*u.cm # while no significant attenuation observed for H6 testbeam + c.MuFilter.VTAttenuationLength = 999*u.cm # Veto 3, no significant attenuation observed in 2024-2025 data + c.MuFilter.VandUpSiPMcalibrationL = 50.*1000. # 1.65 MeV = 41 qcd over 6 Large SiPMs(one side) + c.MuFilter.VandUpSiPMcalibrationS = 0. # no MIP signal for small SiPMs, delayed and compromised response in general c.MuFilter.DsSiPMcalibration = 25.*1000. c.MuFilter.timeResol = 150.*u.picosecond - c.MuFilter.VandUpPropSpeed = 12.5*u.cm/u.nanosecond - c.MuFilter.DsPropSpeed = 14.3*u.cm/u.nanosecond + c.MuFilter.VandUpPropSpeed = 13.6*u.cm/u.nanosecond + c.MuFilter.DsPropSpeed = 15.1*u.cm/u.nanosecond c.Floor = AttrDict(z=48000.*u.cm) # to place tunnel in SND_@LHC coordinate system c.Floor.DX = 1.0*u.cm diff --git a/geometry/sndLHC_HXgeom_config.py b/geometry/sndLHC_HXgeom_config.py index be48a76012..3b57be8b3d 100644 --- a/geometry/sndLHC_HXgeom_config.py +++ b/geometry/sndLHC_HXgeom_config.py @@ -256,15 +256,16 @@ c.MuFilter.DS4ShiftX = 1.39 * u.cm #digitization parameters - c.MuFilter.DsAttenuationLength = 350 * u.cm # values between 300 cm and 400cm observed for H6 testbeam - c.MuFilter.DsTAttenuationLength = 700 * u.cm # top readout with mirror on bottom - c.MuFilter.VandUpAttenuationLength = 999 * u.cm # no significante attenuation observed for H6 testbeam - c.MuFilter.VandUpSiPMcalibrationL = 25.*1000. # 1.65 MeV = 41 qcd - c.MuFilter.VandUpSiPMcalibrationS = 25.*1000. + c.MuFilter.DsAttenuationLength = 230*u.cm # values between 130cm and 330cm are observed for TI18 in years 2022-2025, but between 300 cm and 400cm for H6 testbeam + c.MuFilter.DsTAttenuationLength = 700*u.cm # top readout with mirror on bottom, TI18 and H6 observables agree + c.MuFilter.VandUpAttenuationLength = 210*u.cm # while no significant attenuation observed for H6 testbeam + c.MuFilter.VTAttenuationLength = 999*u.cm # Veto 3, no significant attenuation observed in 2024-2025 data + c.MuFilter.VandUpSiPMcalibrationL = 50.*1000. # 1.65 MeV = 41 qcd over 6 Large SiPMs(one side) + c.MuFilter.VandUpSiPMcalibrationS = 0. # no MIP signal for small SiPMs, delayed and compromised response in general c.MuFilter.DsSiPMcalibration = 25.*1000. c.MuFilter.timeResol = 150.*u.picosecond - c.MuFilter.VandUpPropSpeed = 12.5*u.cm/u.nanosecond - c.MuFilter.DsPropSpeed = 14.3*u.cm/u.nanosecond + c.MuFilter.VandUpPropSpeed = 13.6*u.cm/u.nanosecond + c.MuFilter.DsPropSpeed = 15.1*u.cm/u.nanosecond c.Floor = AttrDict(z=48000.*u.cm) # to place tunnel in SND_@LHC coordinate system c.Floor.DX = 1.0*u.cm diff --git a/geometry/sndLHC_TI18geom_config.py b/geometry/sndLHC_TI18geom_config.py index 48db4b0e68..27136b94b5 100644 --- a/geometry/sndLHC_TI18geom_config.py +++ b/geometry/sndLHC_TI18geom_config.py @@ -317,15 +317,16 @@ c.MuFilter.DS4ShiftX = 0.0 #digitization parameters - c.MuFilter.DsAttenuationLength = 350 * u.cm # values between 300 cm and 400cm observed for H6 testbeam - c.MuFilter.DsTAttenuationLength = 700 * u.cm # top readout with mirror on bottom - c.MuFilter.VandUpAttenuationLength = 999 * u.cm # no significante attenuation observed for H6 testbeam - c.MuFilter.VandUpSiPMcalibrationL = 25.*1000. # 1.65 MeV = 41 qcd - c.MuFilter.VandUpSiPMcalibrationS = 25.*1000. + c.MuFilter.DsAttenuationLength = 230*u.cm # values between 130cm and 330cm are observed for TI18 in years 2022-2025, but between 300 cm and 400cm for H6 testbeam + c.MuFilter.DsTAttenuationLength = 700*u.cm # top readout with mirror on bottom, TI18 and H6 observables agree + c.MuFilter.VandUpAttenuationLength = 210*u.cm # while no significant attenuation observed for H6 testbeam + c.MuFilter.VTAttenuationLength = 999*u.cm # Veto 3, no significant attenuation observed in 2024-2025 data + c.MuFilter.VandUpSiPMcalibrationL = 50.*1000. # 1.65 MeV = 41 qcd over 6 Large SiPMs(one side) + c.MuFilter.VandUpSiPMcalibrationS = 0. # no MIP signal for small SiPMs, delayed and compromised response in general c.MuFilter.DsSiPMcalibration = 25.*1000. c.MuFilter.timeResol = 150.*u.picosecond - c.MuFilter.VandUpPropSpeed = 12.5*u.cm/u.nanosecond - c.MuFilter.DsPropSpeed = 14.9*u.cm/u.nanosecond + c.MuFilter.VandUpPropSpeed = 13.6*u.cm/u.nanosecond + c.MuFilter.DsPropSpeed = 15.1*u.cm/u.nanosecond c.Floor = AttrDict(z=48000.*u.cm) # to place tunnel in SND_@LHC coordinate system c.Floor.DX = 1.0*u.cm diff --git a/macro/MufluxReco.py b/macro/MufluxReco.py index eea4ab39ff..18ea26f370 100644 --- a/macro/MufluxReco.py +++ b/macro/MufluxReco.py @@ -97,7 +97,7 @@ def mem_monitor(): # check if simulation or raw data f=ROOT.TFile.Open(outFile) -if f.cbmsim.FindBranch('MCTrack'): simulation = True +if f.Get("cbmsim").FindBranch('MCTrack'): simulation = True else: simulation = False f.Close() diff --git a/macro/run_simScript.py b/macro/run_simScript.py index 3b871a8859..79cf1f82d7 100755 --- a/macro/run_simScript.py +++ b/macro/run_simScript.py @@ -673,6 +673,7 @@ # checking for overlaps if checking4overlaps: + ROOT.gROOT.SetWebDisplay("off") # Workaround for https://github.com/root-project/root/issues/18881 fGeo = ROOT.gGeoManager fGeo.SetNmeshPoints(10000) fGeo.CheckOverlaps(0.1) # 1 micron takes 5minutes @@ -707,7 +708,7 @@ nm = ff.GetName().split('/') if nm[len(nm)-1] == check: fin = ff if not fin: fin = ROOT.TFile.Open(outFile) - t = fin.cbmsim + t = fin.Get("cbmsim") fout = ROOT.TFile(tmpFile,'recreate') sTree = t.CloneTree(0) nEvents = 0 diff --git a/python/SndlhcDigi.py b/python/SndlhcDigi.py index 55d83ac8cf..5900982843 100644 --- a/python/SndlhcDigi.py +++ b/python/SndlhcDigi.py @@ -16,7 +16,7 @@ def __init__(self,fout, makeClusterScifi): outdir=os.getcwd() outfile=outdir+"/"+fout self.fn = ROOT.TFile(fout,'update') - self.sTree = self.fn.cbmsim + self.sTree = self.fn.Get("cbmsim") # event header self.header = ROOT.SNDLHCEventHeader() @@ -103,11 +103,13 @@ def digitizeScifi(self): allWeights.push_back(p[1]) aHit = ROOT.sndScifiHit(detID,allPoints,allWeights) if self.digiScifi.GetSize() == index: self.digiScifi.Expand(index+100) - self.digiScifi[index]=aHit + aHit_TCA = self.digiScifi.ConstructedAt(index) + ROOT.std.swap(aHit, aHit_TCA) index+=1 for k in mcPoints[detID]: mcLinks.Add(detID,k, mcPoints[detID][k]/norm[detID]) - self.digiScifi2MCPoints[0]=mcLinks + mcLinks_TCA = self.digiScifi2MCPoints.ConstructedAt(0) + ROOT.std.swap(mcLinks, mcLinks_TCA) def digitizeMuFilter(self): hitContainer = {} @@ -134,11 +136,13 @@ def digitizeMuFilter(self): aHit = ROOT.MuFilterHit(detID,allPoints) if self.digiMuFilter.GetSize() == index: self.digiMuFilter.Expand(index+100) - self.digiMuFilter[index]=aHit + aHit_TCA = self.digiMuFilter.ConstructedAt(index) + ROOT.std.swap(aHit, aHit_TCA) index+=1 for k in mcPoints[detID]: mcLinks.Add(detID,k, mcPoints[detID][k]/norm[detID]) - self.digiMuFilter2MCPoints[0]=mcLinks + mcLinks_TCA = self.digiMuFilter2MCPoints.ConstructedAt(0) + ROOT.std.swap(mcLinks, mcLinks_TCA) def clusterScifi(self): index = 0 @@ -167,7 +171,8 @@ def clusterScifi(self): for aHit in tmp: hitlist.push_back( self.sTree.Digi_ScifiHits[hitDict[aHit]],) aCluster = ROOT.sndCluster(first,N,hitlist,self.scifiDet) if self.clusScifi.GetSize() == index: self.clusScifi.Expand(index+10) - self.clusScifi[index]=aCluster + aCluster_TCA = self.clusScifi.ConstructedAt(index) + ROOT.std.swap(aCluster, aCluster_TCA) index+=1 if c!=hitList[last]: ncl+=1 diff --git a/python/SndlhcMuonReco.py b/python/SndlhcMuonReco.py index c12de6491e..d38a8ca06d 100644 --- a/python/SndlhcMuonReco.py +++ b/python/SndlhcMuonReco.py @@ -946,7 +946,8 @@ def Exec(self, opt) : this_track.setRawMeasTimes(pointTimes) this_track.setTrackType(self.track_type) # Save the track in sndRecoTrack format - self.kalman_tracks[i_muon] = this_track + this_track_TCA = self.kalman_tracks.ConstructedAt(i_muon) + ROOT.std.swap(this_track, this_track_TCA) # Delete the Kalman track object theTrack.Delete() diff --git a/python/SndlhcTracking.py b/python/SndlhcTracking.py index ff92e8cb00..6476b00166 100644 --- a/python/SndlhcTracking.py +++ b/python/SndlhcTracking.py @@ -132,7 +132,8 @@ def ExecuteTask(self,option='ScifiDS'): this_track.setRawMeasTimes(pointTimes) this_track.setTrackType(rc.GetUniqueID()) # Store the track in sndRecoTrack format - self.fittedTracks[i_muon] = this_track + this_track_TCA = self.fittedTracks.ConstructedAt(i_muon) + ROOT.std.swap(this_track, this_track_TCA) def DStrack(self): diff --git a/python/rootUtils.py b/python/rootUtils.py index 2b962fee96..05829a0040 100644 --- a/python/rootUtils.py +++ b/python/rootUtils.py @@ -146,7 +146,7 @@ def container_sizes(sTree,perEvent=False): def stripOffBranches(fout): f = TFile(fout) - sTree = f.cbmsim + sTree = f.Get("cbmsim") nEvents = sTree.GetEntries() strip = False oldTargetClass = False @@ -178,7 +178,7 @@ def stripOffBranches(fout): recf.Close() # should do some sanity checks before deleting old file f = TFile(sFile) - sTree = f.cbmsim + sTree = f.Get("cbmsim") if nEvents == sTree.GetEntries(): print("looks ok, could be deleted",os.path.abspath('.')) else: print("stripping failed, keep old file",os.path.abspath('.')) # os.system('mv '+sFile +' '+fout) diff --git a/shipLHC/EmulsionDet.cxx b/shipLHC/EmulsionDet.cxx index 33a5734ceb..ec7ecaed0f 100644 --- a/shipLHC/EmulsionDet.cxx +++ b/shipLHC/EmulsionDet.cxx @@ -195,6 +195,9 @@ void EmulsionDet::ConstructGeometry() int NTungstenPlatesTB24 = conf_ints["EmulsionDet/n_tungsten_plates_tb24"]; bool testbeam_2024_setup = false; + if (NTungstenPlatesTB24 > 0) { + testbeam_2024_setup = true; + } TGeoVolume *top=gGeoManager->FindVolumeFast("Detector"); if(!top) LOG(ERROR) << "no Detector volume found " ; diff --git a/shipLHC/Large__SND_Logo_white_cut.png b/shipLHC/Large__SND_Logo_white_cut.png new file mode 100644 index 0000000000..a4f8b1b614 Binary files /dev/null and b/shipLHC/Large__SND_Logo_white_cut.png differ diff --git a/shipLHC/MuFilter.cxx b/shipLHC/MuFilter.cxx index 9ab91b3f98..e8ef291fbb 100644 --- a/shipLHC/MuFilter.cxx +++ b/shipLHC/MuFilter.cxx @@ -714,16 +714,18 @@ void MuFilter::GetPosition(Int_t fDetectorID, TVector3& vLeft, TVector3& vRight) } Int_t MuFilter::GetnSides(Int_t detID){ int subsystem = floor(detID/10000)-1; - if (subsystem==0){ - // vertical Veto 3 has the readout on the top only - if (detID>=12000) return conf_ints["MuFilter/VetonSides"]-1; - else {return conf_ints["MuFilter/VetonSides"];} - } - if (subsystem==1){return conf_ints["MuFilter/UpstreamnSides"];} - if (subsystem==2){ - if (detID%1000>59) return conf_ints["MuFilter/DownstreamnSides"]-1; - else {return conf_ints["MuFilter/DownstreamnSides"];} - } + switch (subsystem) { + case 0: + // vertical Veto 3 has the readout on the top only + return (detID>=12000) ? conf_ints["MuFilter/VetonSides"]-1 : conf_ints["MuFilter/VetonSides"]; + case 1: + return conf_ints["MuFilter/UpstreamnSides"]; + case 2: + return (detID%1000>59) ? conf_ints["MuFilter/DownstreamnSides"]-1 : conf_ints["MuFilter/DownstreamnSides"]; + default: + LOG(FATAL) << "Unknown subsystem"; + return -1; + } } ClassImp(MuFilter) diff --git a/shipLHC/MuFilterHit.cxx b/shipLHC/MuFilterHit.cxx index a871785c7b..ed2c0a8fe6 100644 --- a/shipLHC/MuFilterHit.cxx +++ b/shipLHC/MuFilterHit.cxx @@ -1,5 +1,6 @@ #include "MuFilterHit.h" #include "MuFilter.h" +#include "ShipUnit.h" #include "TROOT.h" #include "FairRunSim.h" #include "TGeoNavigator.h" @@ -8,7 +9,6 @@ #include #include -Double_t speedOfLight = TMath::C() *100./1000000000.0 ; // from m/sec to cm/ns // ----- Default constructor ------------------------------------------- MuFilterHit::MuFilterHit() : SndlhcHit() @@ -60,10 +60,10 @@ MuFilterHit::MuFilterHit(Int_t detID, std::vector V) else { if (floor(detID/10000)==1 && nSides==1){ // top readout with mirror on bottom - attLength = 2*MuFilterDet->GetConfParF("MuFilter/VandUpAttenuationLength"); + attLength = 2*MuFilterDet->GetConfParF("MuFilter/VTAttenuationLength"); } else {attLength = MuFilterDet->GetConfParF("MuFilter/VandUpAttenuationLength");} - siPMcalibration = MuFilterDet->GetConfParF("MuFilter/VandUpSiPMcalibration"); + siPMcalibration = MuFilterDet->GetConfParF("MuFilter/VandUpSiPMcalibrationL"); siPMcalibrationS = MuFilterDet->GetConfParF("MuFilter/VandUpSiPMcalibrationS"); propspeed = MuFilterDet->GetConfParF("MuFilter/VandUpPropSpeed"); } @@ -104,7 +104,7 @@ MuFilterHit::MuFilterHit(Int_t detID, std::vector V) // In the SndlhcHit class the 'signals' array starts from 0. for (unsigned int j=0; jGaus(earliestToAL, timeResol); }else{ signals[j] = signalLeft/float(nSiPMs) * siPMcalibration; // most simplest model, divide signal individually. @@ -125,11 +125,12 @@ MuFilterHit::~MuFilterHit() { } // ------------------------------------------------------------------------- // ----- Public method GetEnergy ------------------------------------------- -Float_t MuFilterHit::GetEnergy() +Float_t MuFilterHit::GetEnergy(Bool_t use_small_sipms) { // to be calculated from digis and calibration constants, missing! Float_t E = 0; for (unsigned int j=0; j1){ E+=signals[j+nSiPMs];} } @@ -137,20 +138,21 @@ Float_t MuFilterHit::GetEnergy() } bool MuFilterHit::isVertical(){ - if ( (floor(fDetectorID/10000)==3&&fDetectorID%1000>59) || - (floor(fDetectorID/10000)==1&&int(fDetectorID/1000)%10==2) ) { + if ( (GetSystem()==3&&fDetectorID%1000>59) || + (GetSystem()==1&&GetPlane()==2) ) { return kTRUE; } else{return kFALSE;} } bool MuFilterHit::isShort(Int_t i){ - if (i%8==2 || i%8==5) {return kTRUE;} + // only US has short SiPMs + if (GetSystem()==2 && (i%8==2 || i%8==5)) {return kTRUE;} else{return kFALSE;} } // ----- Public method Get List of signals ------------------------------------------- -std::map MuFilterHit::GetAllSignals(Bool_t mask,Bool_t positive) +std::map MuFilterHit::GetAllSignals(Bool_t mask,Bool_t positive,Bool_t use_small_sipms) { std::map allSignals; for (unsigned int s=0; s MuFilterHit::GetAllSignals(Bool_t mask,Bool_t positive) if (signals[channel]<-900){continue;} if (signals[channel]> 0 || !positive){ if (!fMasked[channel] || !mask){ + if (!isShort(channel) || use_small_sipms){ allSignals[channel] = signals[channel]; } + } } } } @@ -168,16 +172,18 @@ std::map MuFilterHit::GetAllSignals(Bool_t mask,Bool_t positive) } // ----- Public method Get List of time measurements ------------------------------------------- -std::map MuFilterHit::GetAllTimes(Bool_t mask) +std::map MuFilterHit::GetAllTimes(Bool_t mask,Bool_t positive,Bool_t use_small_sipms) { std::map allTimes; for (unsigned int s=0; s 0){ + if (signals[channel]> 0 || !positive){ if (!fMasked[channel] || !mask){ + if (!isShort(channel) || use_small_sipms){ allTimes[channel] = times[channel]; } + } } } } @@ -185,7 +191,7 @@ std::map MuFilterHit::GetAllTimes(Bool_t mask) } // ----- Public method Get time difference mean Left - mean Right ----------------- -Float_t MuFilterHit::GetDeltaT(Bool_t mask) +Float_t MuFilterHit::GetDeltaT(Bool_t mask,Bool_t positive,Bool_t use_small_sipms) // based on mean TDC measured on Left and Right { Float_t mean[] = {0,0}; @@ -194,11 +200,13 @@ Float_t MuFilterHit::GetDeltaT(Bool_t mask) for (unsigned int s=0; s 0){ + if (signals[channel]> 0 || !positive){ if (!fMasked[channel] || !mask){ + if (!isShort(channel) || use_small_sipms){ mean[s] += times[channel]; count[s] += 1; } + } } } } @@ -207,7 +215,7 @@ Float_t MuFilterHit::GetDeltaT(Bool_t mask) } return dT; } -Float_t MuFilterHit::GetFastDeltaT(Bool_t mask) +Float_t MuFilterHit::GetFastDeltaT(Bool_t mask,Bool_t positive,Bool_t use_small_sipms) // based on fastest (earliest) TDC measured on Left and Right { Float_t first[] = {1E20,1E20}; @@ -215,10 +223,12 @@ Float_t MuFilterHit::GetFastDeltaT(Bool_t mask) for (unsigned int s=0; s 0){ + if (signals[channel]> 0 || !positive){ if (!fMasked[channel] || !mask){ - if (times[channel] (gROOT->GetListOfGlobals()->FindObject("MuFilter")); - if (floor(fDetectorID/10000==3)) { - dL = MuFilterDet->GetConfParF("MuFilterDet/DownstreamBarX") / MuFilterDet->GetConfParF("MuFilter/DsPropSpeed");} - else if (floor(fDetectorID/10000==2)) { - dL = MuFilterDet->GetConfParF("MuFilterDet/UpstreamBarX") / MuFilterDet->GetConfParF("MuFilter/VandUpPropSpeed");} + if (GetSystem()==3) { + dL = MuFilterDet->GetConfParF("MuFilter/DownstreamBarX") / MuFilterDet->GetConfParF("MuFilter/DsPropSpeed");} + else if (GetSystem()==2) { + dL = MuFilterDet->GetConfParF("MuFilter/UpstreamBarX") / MuFilterDet->GetConfParF("MuFilter/VandUpPropSpeed");} else { - dL = MuFilterDet->GetConfParF("MuFilterDet/VetoBarX") / MuFilterDet->GetConfParF("MuFilter/VandUpPropSpeed");} + dL = MuFilterDet->GetConfParF("MuFilter/VetoBarX") / MuFilterDet->GetConfParF("MuFilter/VandUpPropSpeed");} for (unsigned int s=0; s 0){ + if (signals[channel]> 0 || !positive){ if (!fMasked[channel] || !mask){ - mean[s] += times[channel]; - count[s] += 1; - } + if (!isShort(channel) || use_small_sipms){ + mean[s] += times[channel]; + count[s] += 1; + } + } } } } if (count[0]>0 && count[1]>0) { - dT = (mean[0]/count[0] + mean[1]/count[1])/2.*6.25 - dL/2.; // TDC to ns = 6.25 + dT = (mean[0]/count[0] + mean[1]/count[1])/2.*ShipUnit::snd_TDC2ns - dL/2.; // TDC to ns = 6.25 } return dT; } +// ----- Public method Get position of impact point along the bar ----------------- +Float_t MuFilterHit::GetImpactXpos(Bool_t mask,Bool_t positive,Bool_t use_small_sipms,Bool_t isMC) +{ + if ( nSides!=2 ){ + return -999.; + } + Float_t dT = GetDeltaT(mask,positive,use_small_sipms); + if (dT==-999.){ + return -999.; + } + + MuFilter* MuFilterDet = dynamic_cast (gROOT->GetListOfGlobals()->FindObject("MuFilter")); + Float_t bar_length = MuFilterDet->GetConfParF("MuFilter/UpstreamBarX"); + Float_t signal_speed = MuFilterDet->GetConfParF("MuFilter/VandUpPropSpeed"); + if (GetSystem()==3) { + signal_speed = MuFilterDet->GetConfParF("MuFilter/DsPropSpeed"); + bar_length = MuFilterDet->GetConfParF("MuFilter/DownstreamBarX"); + } + else if (GetSystem()==2) { + bar_length = MuFilterDet->GetConfParF("MuFilter/UpstreamBarX"); + } + else { + bar_length = MuFilterDet->GetConfParF("MuFilter/VetoBarX"); + } + double timeConversion = 1.; + if (!isMC) timeConversion = ShipUnit::snd_TDC2ns; + return 0.5*(bar_length + dT*timeConversion*signal_speed); +} + std::map MuFilterHit::SumOfSignals(Bool_t mask) { /* use cases, for Veto and DS small/large ignored @@ -277,7 +318,7 @@ std::map MuFilterHit::SumOfSignals(Bool_t mask) for (unsigned int s=0; s 0){ + if (signals[channel]> 0){ // makes sense to sum up positive signals only if (!fMasked[channel] || !mask){ if (s==0 and !isShort(j)){theSumL+= signals[channel];} if (s==0 and isShort(j)){theSumLS+= signals[channel];} diff --git a/shipLHC/MuFilterHit.h b/shipLHC/MuFilterHit.h index c448a7f015..4c5e1e0363 100644 --- a/shipLHC/MuFilterHit.h +++ b/shipLHC/MuFilterHit.h @@ -27,21 +27,21 @@ class MuFilterHit : public SndlhcHit /** Output to screen **/ void Print() const; - Float_t GetEnergy(); - Float_t SumOfSignals(char* opt,Bool_t mask=kTRUE); + Float_t GetEnergy(Bool_t use_small_sipms=kFALSE); std::map SumOfSignals(Bool_t mask=kTRUE); - std::map GetAllSignals(Bool_t mask=kTRUE,Bool_t positive=kTRUE); - std::map GetAllTimes(Bool_t mask=kTRUE); - Float_t GetDeltaT(Bool_t mask=kTRUE); - Float_t GetFastDeltaT(Bool_t mask=kTRUE); - Float_t GetImpactT(Bool_t mask=kTRUE); + std::map GetAllSignals(Bool_t mask=kTRUE,Bool_t positive=kTRUE,Bool_t use_small_sipms=kFALSE); + std::map GetAllTimes(Bool_t mask=kTRUE,Bool_t positive=kTRUE,Bool_t use_small_sipms=kFALSE); + Float_t GetDeltaT(Bool_t mask=kTRUE,Bool_t positive=kTRUE,Bool_t use_small_sipms=kFALSE); + Float_t GetFastDeltaT(Bool_t mask=kTRUE,Bool_t positive=kTRUE,Bool_t use_small_sipms=kFALSE); + Float_t GetImpactT(Bool_t mask=kTRUE,Bool_t positive=kTRUE,Bool_t use_small_sipms=kFALSE); + Float_t GetImpactXpos(Bool_t mask=kTRUE,Bool_t positive=kTRUE,Bool_t use_small_sipms=kFALSE,Bool_t isMC=kFALSE); bool isValid() const {return flag;} bool isMasked(Int_t i) const {return fMasked[i];} void SetMasked(Int_t i) {fMasked[i]=kTRUE;} int GetSystem(){return floor(fDetectorID/10000);} int GetPlane(){return int(fDetectorID/1000)%10;} bool isVertical(); - bool isShort(Int_t); + bool isShort(Int_t);// short==small sipm private: Float_t flag; ///< flag diff --git a/shipLHC/MuFilterPoint.h b/shipLHC/MuFilterPoint.h index c6784c78f3..df33249e72 100644 --- a/shipLHC/MuFilterPoint.h +++ b/shipLHC/MuFilterPoint.h @@ -38,17 +38,17 @@ class MuFilterPoint : public FairMCPoint Int_t PdgCode() const {return fPdgCode;} + + /** Copy constructor **/ + MuFilterPoint(const MuFilterPoint& point); + + MuFilterPoint operator=(const MuFilterPoint& point); private: Int_t fPdgCode; - /** Copy constructor **/ - - MuFilterPoint(const MuFilterPoint& point); - MuFilterPoint operator=(const MuFilterPoint& point); - ClassDef(MuFilterPoint,1) }; diff --git a/shipLHC/ana_thermalNeutrons.py b/shipLHC/ana_thermalNeutrons.py index f48fd7203a..181dd2a181 100644 --- a/shipLHC/ana_thermalNeutrons.py +++ b/shipLHC/ana_thermalNeutrons.py @@ -30,7 +30,7 @@ def count(hFile): ut.bookHist(h,'xyz','', 100,-0.1,0.1,100,-0.1,0.1,200,-1.,1.) ut.bookHist(h,'dxyz','',100,-0.1,0.1,100,-0.1,0.1,200,-1.,1.) Npassed = 0 - for sTree in f.cbmsim: + for sTree in f.Get("cbmsim"): N = sTree.MCTrack[0] Ekin = N.GetP()**2/(2*N.GetMass())*1000. logEkin = ROOT.TMath.Log10(Ekin) @@ -152,7 +152,7 @@ def absorptionLength(plotOnly=True): myPrint(h['T'+X],'fracEveWith'+X) # fntuple = ROOT.TFile.Open('neutronsTI18.root') - nt=fntuple.nt + nt=fntuple.Get("nt") ROOT.gROOT.cd() tcanv = 'TFig12' if tcanv in h: h.pop(tcanv) @@ -201,7 +201,7 @@ def absorptionLength(plotOnly=True): n = 0 RateIntegrated = 0 RateIntegratedW = 0 - for nt in fntuple.nt: + for nt in fntuple.Get("nt"): E = (nt.Eleft+nt.Eright)/2. dE = nt.Eright - nt.Eleft h['Fig12'].SetPoint(n,ROOT.TMath.Log10(E),ROOT.TMath.Log10(nt.N*E)) @@ -312,7 +312,7 @@ def reactions(hFile,distance=1E10): ut.bookHist(h,'E',';log10(Ekin [MeV])',1000,-12.,4.) stats = {} n=-1 - for sTree in f.cbmsim: + for sTree in f.Get("cbmsim"): n+=1 daughters = [] for m in sTree.MCTrack: @@ -397,7 +397,7 @@ def flukaRateIntegrated(save=False): h['neutronRate'] = ROOT.TGraph() h['N'] = ROOT.TGraph() n = 0 - for nt in fntuple.nt: + for nt in fntuple.Get("nt"): # nt.N = [cm/MeV] E = (nt.Eleft+nt.Eright)/2. dE = nt.Eright - nt.Eleft @@ -516,8 +516,8 @@ def coldBox(plotOnly=True,pas=''): ut.bookHist(h,'multN','mult neutrons',100,-0.5,99.5) flukaRateIntegrated() - Nsim = f.cbmsim.GetEntries() - for sTree in f.cbmsim: + Nsim = f.Get("cbmsim").GetEntries() + for sTree in f.Get("cbmsim"): rc=h['multN'].Fill(sTree.MCTrack.GetEntries()) neutron = sTree.MCTrack[0] start = ROOT.TVector3(neutron.GetStartX(),neutron.GetStartY(),neutron.GetStartZ()) diff --git a/shipLHC/modifyGeoFileDict.py b/shipLHC/modifyGeoFileDict.py index 9bc7a66858..bf4586b6b9 100644 --- a/shipLHC/modifyGeoFileDict.py +++ b/shipLHC/modifyGeoFileDict.py @@ -9,17 +9,18 @@ theClient = client.FileSystem('root://eospublic.cern.ch') commonPath = "/eos/experiment/sndlhc/convertedData/physics/" -'''supportedGeoFiles = ["geofile_sndlhc_TI18_V1_06July2022.root","geofile_sndlhc_TI18_V2_12July2022.root","geofile_sndlhc_TI18_V3_08August2022.root", - "geofile_sndlhc_TI18_V4_10August2022.root","geofile_sndlhc_TI18_V5_14August2022.root","geofile_sndlhc_TI18_V6_08October2022.root", - "geofile_sndlhc_TI18_V7_22November2022.root"] -''' -supportedGeoFiles = {"geofile_sndlhc_TI18_V7_22November2022.root":commonPath+"2022/", - "geofile_sndlhc_TI18_V0_2024.root":commonPath+"2024/"} +supportedGeoFiles = {} +supported_years = [2022, 2023, 2024, 2025, 2026] +for year in supported_years: + supportedGeoFiles["geofile_sndlhc_TI18_V0_"+str(year)+".root"] = commonPath+str(year)+"/" def modifyDicts(year=2024): + if year not in supported_years: + print("Geometry for the required year "+str(year)+"is not supported.\n\ + The list of supported years is:", supported_years) + return for geoFileName in supportedGeoFiles: - if (year == 2024 and geoFileName.find(str(year)) > 0) or \ - (year != 2024 and geoFileName.find('2022') > 0): + if geoFileName.find(str(year)) > 0: # override locally existing file with the same name status = theClient.copy(supportedGeoFiles[geoFileName]+geoFileName, os.getcwd()+'/'+geoFileName, force=True) @@ -36,7 +37,9 @@ def modifyDicts(year=2024): mufi_spatial_aligment_consts = { 2022: [ 0.11,-0.04, 0.00, 0.10, 0.26, 0.24, 0.31, 0.34, 0.43, 1.13, 0.53, 1.31, 0.61, 1.35, 1.39], 2023: [ 0.15, 0.03, 0.00,-0.01, 0.09,-0.01, 0.06, 0.06, 0.49, 0.86, 0.20, 0.98, 0.24, 1.01, 1.11], - 2024: [-0.06, 0.04, 0.65,-0.15,-0.12,-0.26,-0.29,-0.36, 0.00, 0.35,-0.37, 0.38,-0.41, 0.30, 0.35] + 2024: [-0.06, 0.04, 0.65,-0.15,-0.12,-0.26,-0.29,-0.36, 0.00, 0.35,-0.37, 0.38,-0.41, 0.30, 0.35], + 2025: [-0.04, 0.05, 0.66,-0.16,-0.15,-0.31,-0.34,-0.45,-0.10, 0.39,-0.48, 0.43,-0.91,-0.36,-2.02], + 2026: [-0.06, 0.06, 0.62,-0.17,-0.17,-0.30,-0.34,-0.43,-0.08, 0.44,-0.46, 0.48,-0.89,-0.34,-2.01] } mufi_spatial_aligment_keys = ['Veto1ShiftY','Veto2ShiftY','Veto3ShiftX', 'US1ShiftY','US2ShiftY','US3ShiftY','US4ShiftY','US5ShiftY', @@ -69,13 +72,28 @@ def modifyDicts(year=2024): constants['t_9692'] = [-6.58,-6.67,-6.76,-6.95,-6.47,-6.53,-6.07,-6.20,-6.22,-6.39,-6.01,-6.18,-8.23,-8.26,-8.43,-8.65,-8.19,-8.27,-8.20,-8.35] constants['t_9882'] = [-6.22,-6.13,-6.37,-6.25,-6.02,-6.18,-5.65,-5.58,-5.76,-5.79,-5.53,-5.75,-7.81,-7.69,-8.02,-8.07,-7.74,-7.85,-7.81,-7.98] constants['t_10012']= [-6.70,-6.84,-6.73,-6.86,-6.50,-6.64,-6.04,-6.19,-6.21,-6.41,-6.05,-6.24,-8.23,-8.38,-8.40,-8.64,-8.22,-8.38,-8.31,-8.52] - slopes_dict_2022 = {"t_0":0.000, "t_4361":0.082, "t_5117":0.085, "t_5478":0.082, "t_6208":0.086, "t_6443":0.082, "t_6677":0.084} - # new constants for 2024 and later on - slopes_dict_2024 = {"t_0":0.000, "t_7649":0.081, "t_8318":0.082, "t_8583":0.081, "t_8942":0.080, - "t_9156":0.083, "t_9286":0.082, "t_9379":0.083, "t_9462":0.083, "t_9613":0.082, - "t_9692":0.078, "t_9882":0.084, "t_10012":0.085} - if year == 2024: slopes_dict = slopes_dict_2024 - else: slopes_dict = slopes_dict_2022 + # 2025 + constants['t_10423']= [-6.61,-6.69,-6.68,-6.77,-6.44,-6.58,-5.98,-6.06,-6.11,-6.29,-5.98,-6.16,-8.16,-8.17,-8.24,-8.43,-8.13,-8.22,-8.26,-8.41] + constants['t_11158']= [-6.60,-6.68,-6.79,-6.86,-6.51,-6.60,-6.08,-6.13,-6.20,-6.32,-6.04,-6.20,-8.19,-8.26,-8.34,-8.49,-8.19,-8.29,-8.31,-8.44] + constants['t_11576']= [-6.41,-6.48,-6.72,-6.78,-6.43,-6.53,-5.94,-6.06,-6.19,-6.36,-6.00,-6.17,-8.17,-8.20,-8.30,-8.49,-8.08,-8.16,-8.11,-8.21] + constants['t_11676']= [-6.42,-6.47,-6.74,-6.79,-6.43,-6.53,-5.95,-6.06,-6.21,-6.35,-6.01,-6.17,-8.16,-8.22,-8.32,-8.50,-8.09,-8.17,-8.13,-8.22] + constants['t_11795']= [-6.43,-6.45,-6.71,-6.79,-6.43,-6.52,-5.99,-6.09,-6.15,-6.32,-5.98,-6.13,-8.17,-8.25,-8.29,-8.48,-8.08,-8.14,-8.16,-8.24] + constants['t_12021']= [-6.46,-6.53,-6.79,-6.89,-6.49,-6.62,-6.02,-6.13,-6.24,-6.40,-6.05,-6.22,-8.21,-8.28,-8.38,-8.57,-8.15,-8.25,-8.22,-8.32] + constants['t_12123']= [-6.60,-6.65,-6.74,-6.84,-6.46,-6.57,-6.06,-6.13,-6.15,-6.33,-6.00,-6.16,-8.20,-8.26,-8.30,-8.49,-8.13,-8.23,-8.27,-8.40] + #2026 + constants['t_12794']= [-6.52,-6.59,-6.70,-6.77,-6.42,-6.59,-5.89,-5.98,-6.19,-6.34,-5.97,-6.18,-8.15,-8.19,-8.25,-8.45,-8.08,-8.19,-8.08,-8.26] + constants['t_13118']= [-6.50,-6.56,-6.66,-6.74,-6.39,-6.51,-5.88,-5.99,-6.16,-6.31,-5.95,-6.13,-8.13,-8.16,-8.20,-8.39,-8.04,-8.20,-8.05,-8.23] + ds_time_aligment_consts = { + 2022: {"t_0":0.000, "t_4361":0.082, "t_5117":0.085}, + 2023: {"t_0":0.000, "t_5478":0.082, "t_6208":0.086, "t_6443":0.082, "t_6677":0.084}, + 2024: {"t_0":0.000, "t_7649":0.081, "t_8318":0.082, "t_8583":0.081, "t_8942":0.080, + "t_9156":0.083, "t_9286":0.082, "t_9379":0.083, "t_9462":0.083, "t_9613":0.082, + "t_9692":0.078, "t_9882":0.084, "t_10012":0.085}, + 2025: {"t_0": 0.000, "t_10423":0.083, "t_11158":0.082, "t_11576":0.079, "t_11676":0.080, + "t_11795": 0.079, "t_12021":0.077, "t_12123":0.081}, + 2026: {"t_0": 0.000, "t_12794":0.082, "t_13118":0.081} + } + slopes_dict = ds_time_aligment_consts[year] #time delay corrections first order, only for DS at the moment for p in slopes_dict.keys(): setattr(sGeo.MuFilter,'DSTcorslope'+p,slopes_dict[p]) @@ -182,12 +200,63 @@ def modifyDicts(year=2024): -1.119*u.ns, 0.000*u.ns, -0.005*u.ns, 0.582*u.ns, 0.875*u.ns, -0.390*u.ns, 1.059*u.ns, -0.406*u.ns, 0.000*u.ns, -1.319*u.ns, 0.337*u.ns, -2.012*u.ns, -1.059*u.ns, -0.618*u.ns, -0.825*u.ns, 0.000*u.ns, -1.205*u.ns, -0.952*u.ns, 0.241*u.ns, -1.363*u.ns, 1.123*u.ns ] +#2025 + constants['t_10423']=[ 0.000*u.ns, 0.000*u.ns, -0.337*u.ns, -0.509*u.ns, 0.125*u.ns, -0.442*u.ns, -0.049*u.ns, + -0.936*u.ns, 0.000*u.ns, 0.238*u.ns, -0.780*u.ns, -0.484*u.ns, 0.172*u.ns, -0.326*u.ns, + -1.097*u.ns, 0.000*u.ns, -0.028*u.ns, 0.543*u.ns, 0.764*u.ns, -0.568*u.ns, 0.962*u.ns, + -0.425*u.ns, 0.000*u.ns, -1.339*u.ns, 0.324*u.ns, -2.087*u.ns, -1.170*u.ns, -0.765*u.ns, + -0.768*u.ns, 0.000*u.ns, -1.312*u.ns, -1.038*u.ns, 0.204*u.ns, -1.518*u.ns, 0.960*u.ns] + constants['t_11158']=[ 0.000*u.ns, 0.000*u.ns, -0.377*u.ns, -0.542*u.ns, 0.115*u.ns, -0.456*u.ns, -0.184*u.ns, + -0.982*u.ns, 0.000*u.ns, 0.233*u.ns, -0.827*u.ns, -1.809*u.ns, 0.228*u.ns, -0.281*u.ns, + -1.140*u.ns, 0.000*u.ns, -0.020*u.ns, 0.606*u.ns, 0.813*u.ns, -0.462*u.ns, 1.016*u.ns, + -1.687*u.ns, 0.000*u.ns, -0.101*u.ns, 1.552*u.ns, -0.874*u.ns, 0.028*u.ns, 0.478*u.ns, + -0.791*u.ns, 0.000*u.ns, -1.310*u.ns, -1.055*u.ns, 0.445*u.ns, -1.393*u.ns, 0.794*u.ns] + constants['t_11576']=[ 0.000*u.ns, 0.000*u.ns, -0.373*u.ns, -0.541*u.ns, 0.143*u.ns, -0.430*u.ns, -0.117*u.ns, + -0.939*u.ns, 0.000*u.ns, 0.213*u.ns, -0.861*u.ns, -0.396*u.ns, 0.231*u.ns, -0.286*u.ns, + -1.118*u.ns, 0.000*u.ns, -0.005*u.ns, 0.569*u.ns, 0.823*u.ns, -0.506*u.ns, 1.050*u.ns, + -0.455*u.ns, 0.000*u.ns, -1.324*u.ns, 0.337*u.ns, -2.054*u.ns, -1.458*u.ns, -0.756*u.ns, + -0.790*u.ns, 0.000*u.ns, -1.279*u.ns, -1.031*u.ns, -0.108*u.ns, -1.363*u.ns, 1.033*u.ns] + constants['t_11676']=[ 0.000*u.ns, 0.000*u.ns, -0.374*u.ns, -0.527*u.ns, 0.144*u.ns, -0.406*u.ns, -0.090*u.ns, + -1.324*u.ns, 0.000*u.ns, 0.618*u.ns, -0.929*u.ns, -0.019*u.ns, 0.613*u.ns, 0.056*u.ns, + -1.098*u.ns, 0.000*u.ns, -0.045*u.ns, 0.584*u.ns, 0.915*u.ns, -0.543*u.ns, 1.006*u.ns, + -0.433*u.ns, 0.000*u.ns, -1.322*u.ns, 0.322*u.ns, -2.106*u.ns, -1.526*u.ns, -0.797*u.ns, + -0.816*u.ns, 0.000*u.ns, -1.254*u.ns, -1.008*u.ns, 0.299*u.ns, -1.388*u.ns, 1.048*u.ns] + constants['t_11795']=[ 0.000*u.ns, 0.000*u.ns, -0.359*u.ns, -0.524*u.ns, 0.124*u.ns, -0.417*u.ns, -0.084*u.ns, + -1.327*u.ns, 0.000*u.ns, 0.628*u.ns, -0.722*u.ns, -0.167*u.ns, 0.608*u.ns, 0.140*u.ns, + -1.078*u.ns, 0.000*u.ns, -0.014*u.ns, 0.569*u.ns, 0.930*u.ns, -0.501*u.ns, 1.000*u.ns, + -0.417*u.ns, 0.000*u.ns, -1.360*u.ns, 0.336*u.ns, -2.079*u.ns, -1.514*u.ns, -0.775*u.ns, + -0.736*u.ns, 0.000*u.ns, -1.317*u.ns, -1.076*u.ns, 0.296*u.ns, -1.402*u.ns, 0.989*u.ns] + constants['t_12021']=[ 0.000*u.ns, 0.000*u.ns, -0.362*u.ns, -0.538*u.ns, 0.199*u.ns, -0.368*u.ns, 0.065*u.ns, + -1.309*u.ns, 0.000*u.ns, 0.621*u.ns, -0.677*u.ns, 0.018*u.ns, 0.642*u.ns, 0.122*u.ns, + -1.107*u.ns, 0.000*u.ns, -0.012*u.ns, 0.655*u.ns, 0.945*u.ns, -0.479*u.ns, 1.043*u.ns, + -0.409*u.ns, 0.000*u.ns, -1.327*u.ns, 0.336*u.ns, -2.151*u.ns, -1.517*u.ns, -0.686*u.ns, + -1.292*u.ns, 0.000*u.ns, -0.719*u.ns, 0.005*u.ns, -2.059*u.ns, -0.669*u.ns, 1.800*u.ns] + constants['t_12123']=[ 0.000*u.ns, 0.000*u.ns, -0.381*u.ns, -0.537*u.ns, 0.162*u.ns, -0.414*u.ns, 0.002*u.ns, + -1.353*u.ns, 0.000*u.ns, 0.614*u.ns, -0.849*u.ns, 0.014*u.ns, 0.634*u.ns, 0.060*u.ns, + -1.092*u.ns, 0.000*u.ns, -0.038*u.ns, 0.568*u.ns, 0.909*u.ns, -0.522*u.ns, 0.968*u.ns, + -0.445*u.ns, 0.000*u.ns, -1.315*u.ns, 0.328*u.ns, -2.095*u.ns, -1.545*u.ns, -0.743*u.ns, + -1.662*u.ns, 0.000*u.ns, -0.410*u.ns, -0.160*u.ns, 1.179*u.ns, -0.530*u.ns, 1.911*u.ns] +#2026 + constants['t_12794']=[ 0.000*u.ns, 0.000*u.ns, -0.325*u.ns, -0.529*u.ns, 0.148*u.ns, -0.449*u.ns, -0.082*u.ns, + -1.478*u.ns, 0.000*u.ns, 0.755*u.ns, -0.195*u.ns, 0.165*u.ns, 0.760*u.ns, 0.255*u.ns, + -1.026*u.ns, 0.000*u.ns, -0.713*u.ns, 0.431*u.ns, 0.806*u.ns, -0.560*u.ns, 0.880*u.ns, + -0.461*u.ns, 0.000*u.ns, -1.388*u.ns, 0.304*u.ns, -2.082*u.ns, -1.489*u.ns, -0.679*u.ns, + -1.676*u.ns, 0.000*u.ns, -0.464*u.ns, -0.170*u.ns, 1.052*u.ns, -0.551*u.ns, 1.918*u.ns] + constants['t_13118']=[ 0.000*u.ns, 0.000*u.ns, -0.330*u.ns, -0.506*u.ns, 0.170*u.ns, -0.408*u.ns, -0.052*u.ns, + -1.409*u.ns, 0.000*u.ns, 0.714*u.ns, -0.277*u.ns, 0.116*u.ns, 0.725*u.ns, 0.222*u.ns, + -0.997*u.ns, 0.000*u.ns, -1.387*u.ns, 0.493*u.ns, 0.817*u.ns, -0.552*u.ns, 0.854*u.ns, + -0.431*u.ns, 0.000*u.ns, -1.428*u.ns, 0.305*u.ns, -2.108*u.ns, -1.494*u.ns, -0.690*u.ns, + -1.652*u.ns, 0.000*u.ns, -0.456*u.ns, -0.151*u.ns, 1.117*u.ns, -0.559*u.ns, 1.926*u.ns] # - scifi_time_tags_2022 = ['t_0', 't_4361','t_5117', 't_5478', 't_6208', 't_6443', 't_6677'] - scifi_time_tags_2024 = ['t_0', 't_7649', 't_8318', 't_8583', 't_8942', 't_9156', 't_9286', - 't_9379', 't_9462', 't_9613', 't_9692', 't_9882', 't_10012'] - if year == 2024: scifi_time_tags = scifi_time_tags_2024 - else: scifi_time_tags = scifi_time_tags_2022 + scifi_time_aligment_consts = { + 2022: ['t_0', 't_4361','t_5117'], + 2023: ['t_0', 't_5478', 't_6208', 't_6443', 't_6677'], + 2024: ['t_0', 't_7649', 't_8318', 't_8583', 't_8942', 't_9156', 't_9286', + 't_9379', 't_9462', 't_9613', 't_9692', 't_9882', 't_10012'], + 2025: ['t_0', 't_10423', 't_11158', 't_11576', 't_11676', 't_11795', 't_12021', 't_12123'], + 2026: ['t_0', 't_12794', 't_13118'] + } + scifi_time_tags = scifi_time_aligment_consts[year] for c in scifi_time_tags: k=0 for s in range(1,6): @@ -521,12 +590,160 @@ def modifyDicts(year=2024): 0*u.mrad, 1.65*u.mrad, 0*u.mrad, 0*u.mrad, -0.79*u.mrad, 0*u.mrad, 0*u.mrad, 2.93*u.mrad, 0*u.mrad] + alignment['t_10423']=[ #2025 emulsion run 19 run_251 + 908.9748*u.um, 573.9855*u.um, 812.3022*u.um, + 86.541*u.um, 123.4473*u.um, 88.648*u.um, + 357.9727*u.um, 241.2571*u.um, 309.628*u.um, + 18.2807*u.um, 87.3165*u.um, 94.6991*u.um, + -1.8476*u.um, 78.9774*u.um, 13.9775*u.um, + 100.7577*u.um, -59.748*u.um, 174.542*u.um, + -100*u.um, 236.436*u.um, 196.9426*u.um, + -11.044*u.um, -185.1257*u.um, 17.2473*u.um, + -593.3074*u.um, -118.6736*u.um, -339.1444*u.um, + -101.3628*u.um, -25.2405*u.um, 184.0245*u.um, + 0*u.mrad, -0.82*u.mrad, 0*u.mrad, + 0*u.mrad, -0.73*u.mrad, 0*u.mrad, + 0*u.mrad, 0*u.mrad, 0*u.mrad, + 0*u.mrad, -0.43*u.mrad, 0*u.mrad, + 0*u.mrad, 0.12*u.mrad, 0*u.mrad] + alignment['t_11158']=[ #2025 emulsion run 20 run_252 + 758.176*u.um, 84.0408*u.um, 221.625*u.um, + 126.9696*u.um, -100.0*u.um, -542.1658*u.um, + 56.0*u.um, -234.5533*u.um, -264.4273*u.um, + 210.2911*u.um, 130.0*u.um, 9.9216*u.um, + -150.0*u.um, -102.0*u.um, -250.0*u.um, + 100.7577*u.um, -59.748*u.um, 304.0*u.um, + -200.0*u.um, 256.0*u.um, 96.9426*u.um, + -550.0*u.um, -600.0*u.um, 50*u.um, + -892.3995*u.um, -108.0*u.um, -489.6239*u.um, + -300.0*u.um, 129.0*u.um, 911.9577*u.um, + 0*u.mrad, -3.16*u.mrad, 0*u.mrad, + 0*u.mrad, -1.13*u.mrad, 0*u.mrad, + 0*u.mrad, 0*u.mrad, 0*u.mrad, + 0*u.mrad, -2.00*u.mrad, 0*u.mrad, + 0*u.mrad, 1.86*u.mrad, 0*u.mrad] + alignment['t_11576']=[ #2025 emulsion run 21 run_253 + 208.20*u.um, 125.00*u.um, 195.00*u.um, + -8.53*u.um, 44.46*u.um, 52.20*u.um, + -8.71*u.um, -30.53*u.um, -63.62*u.um, + -98.60*u.um, 10.00*u.um, -45.00*u.um, + -30.00*u.um, -8.00*u.um, -51.29*u.um, + 180.00*u.um, 98.05*u.um, 228.20*u.um, + 25.19*u.um, 67.28*u.um, 127.54*u.um, + 7.54*u.um, -141.72*u.um, -70.68*u.um, + -5.00*u.um, 0.00*u.um, -50.29*u.um, + -64.05*u.um, 86.94*u.um, 4.89*u.um, + 0.00*u.mrad, 0.01*u.mrad, -0.01*u.mrad, + -0.47*u.mrad, -0.01*u.mrad, 0.08*u.mrad, + 0.90*u.mrad, 0.05*u.mrad, 0.00*u.mrad, + -0.40*u.mrad, -0.01*u.mrad, 0.08*u.mrad, + -0.23*u.mrad, 0.01*u.mrad, 0.04*u.mrad] + alignment['t_11676']=[ #2025 emulsion run 22 run_254 + 233.2526*u.um, 203.3113*u.um, 261.1019*u.um, + -248.7328*u.um, -168.8227*u.um, -98.6278*u.um, + -238.4141*u.um, -238.2084*u.um, -260.8086*u.um, + 72.2430*u.um, 152.4410*u.um, 162.6312*u.um, + 0.0000*u.um, -37.1101*u.um, -101.2731*u.um, + 130.0000*u.um, 13.4088*u.um, 150.0951*u.um, + 223.6887*u.um, 223.8778*u.um, 254.0594*u.um, + -215.8022*u.um, -413.1773*u.um, -344.0464*u.um, + -177.8921*u.um, -214.8997*u.um, -295.3263*u.um, + 125.5983*u.um, 255.9729*u.um, 212.1852*u.um, + -0.4488*u.mrad, -0.3531*u.mrad, 0.3427*u.mrad, + 0.5384*u.mrad, 0.0341*u.mrad, 0.1195*u.mrad, + 0.7200*u.mrad, 0.2000*u.mrad, 0.0000*u.mrad, + -1.0244*u.mrad, -0.5787*u.mrad, 0.0171*u.mrad, + 0.9128*u.mrad, 0.4085*u.mrad, 0.0171*u.mrad] + alignment['t_11795']=[ #2025 emulsion run 23 run_255 + 327.80*u.um, 296.15*u.um, 357.10*u.um, + 49.60*u.um, 102.00*u.um, 25.95*u.um, + -81.10*u.um, -82.15*u.um, -114.40*u.um, + 21.55*u.um, 90.50*u.um, -4.55*u.um, + -45.00*u.um, -36.60*u.um, -83.76*u.um, + 65.00*u.um, 2.20*u.um, 96.00*u.um, + -126.80*u.um, -67.85*u.um, -94.90*u.um, + 29.05*u.um, -20.00*u.um, -29.05*u.um, + 400.60*u.um, 430.85*u.um, 331.80*u.um, + -146.00*u.um, 119.80*u.um, 17.65*u.um, + -0.24*u.mrad, -0.05*u.mrad, -0.30*u.mrad, + -0.06*u.mrad, -0.18*u.mrad, -0.15*u.mrad, + 0.48*u.mrad, 0.00*u.mrad, 0.00*u.mrad, + 0.32*u.mrad, -0.28*u.mrad, 0.23*u.mrad, + -0.20*u.mrad, 0.40*u.mrad, 1.32*u.mrad] + alignment['t_12021']=[ #2025 emulsion run 24 run_256 scifi 4v mat 0 is out + 385.78*u.um, 414.57*u.um, 345.23*u.um, + -188.43*u.um, -125.87*u.um, -75.15*u.um, + -184.34*u.um, -108.25*u.um, -208.02*u.um, + 152.45*u.um, 244.56*u.um, 228.59*u.um, + 1.00*u.um, 48.14*u.um, -56.59*u.um, + 159.00*u.um, 131.02*u.um, 191.70*u.um, + 269.80*u.um, 291.84*u.um, 330.59*u.um, + -433.21*u.um, -465.53*u.um, -428.92*u.um, + -175.32*u.um, -225.90*u.um, -249.68*u.um, + 0.00*u.um, 325.05*u.um, 191.00*u.um, + 0.03*u.mrad, -0.35*u.mrad, -0.13*u.mrad, + 1.22*u.mrad, 0.00*u.mrad, -0.24*u.mrad, + 1.50*u.mrad, 0.10*u.mrad, 0.00*u.mrad, + -1.52*u.mrad, -0.21*u.mrad, 0.21*u.mrad, + 1.83*u.mrad, -0.24*u.mrad, 0.50*u.mrad] + alignment['t_12123']=[ #2025 dummy target run_257 and ion run_258 + 230.00*u.um, 213.75*u.um, 268.00*u.um, + -70.75*u.um, 9.00*u.um, -36.25*u.um, + -101.50*u.um, -91.75*u.um, -118.20*u.um, + -12.25*u.um, 52.50*u.um, -18.60*u.um, + -33.00*u.um, -26.00*u.um, -63.25*u.um, + 58.00*u.um, -7.50*u.um, 56.25*u.um, + -104.00*u.um, -104.25*u.um, -26.50*u.um, + 105.25*u.um, 2.00*u.um, 31.75*u.um, + 238.00*u.um, 216.25*u.um, 214.00*u.um, + -239.25*u.um, 2.00*u.um, -119.75*u.um, + -0.32*u.mrad, -0.01*u.mrad, 0.23*u.mrad, + -0.18*u.mrad, -0.01*u.mrad, 0.23*u.mrad, + 0.45*u.mrad, 0.05*u.mrad, 0.00*u.mrad, + 0.39*u.mrad, -0.09*u.mrad, 0.01*u.mrad, + -0.59*u.mrad, 0.01*u.mrad, 0.00*u.mrad] + alignment['t_12794']=[ #2026 emulsion run 26 run_261 + 148.00*u.um, 180.00*u.um, 176.00*u.um, + -160.00*u.um, -12.00*u.um, -72.00*u.um, + -278.00*u.um, -198.00*u.um, -285.00*u.um, + 75.00*u.um, 225.00*u.um, 117.00*u.um, + 20.00*u.um, 60.00*u.um, -34.00*u.um, + 105.00*u.um, 90.00*u.um, 190.00*u.um, + 188.00*u.um, 240.00*u.um, 290.00*u.um, + -200.00*u.um, -322.00*u.um, -236.00*u.um, + -110.00*u.um, -130.00*u.um, -160.00*u.um, + 210.00*u.um, 400.00*u.um, 345.00*u.um, + -0.38*u.mrad, -0.10*u.mrad, 0.00*u.mrad, + 0.30*u.mrad, 0.10*u.mrad, 0.00*u.mrad, + 0.35*u.mrad, 0.30*u.mrad, 0.00*u.mrad, + -1.00*u.mrad, -0.40*u.mrad, 0.00*u.mrad, + 0.90*u.mrad, 0.30*u.mrad, 0.00*u.mrad] + alignment['t_13118']=[ #2026 emulsion run 27 run_262 + 110.00*u.um, 170.00*u.um, 138.00*u.um, + -88.00*u.um, 64.00*u.um, 0.00*u.um, + -154.00*u.um, -81.00*u.um, -140.00*u.um, + -7.00*u.um, 126.20*u.um, 41.00*u.um, + -5.00*u.um, 55.00*u.um, -47.00*u.um, + 7.00*u.um, -31.00*u.um, 80.20*u.um, + -109.80*u.um, -33.00*u.um, -10.00*u.um, + 152.00*u.um, 58.00*u.um, 96.00*u.um, + 165.80*u.um, 153.00*u.um, 95.00*u.um, + -92.00*u.um, 126.00*u.um, 71.00*u.um, + -0.42*u.mrad, -0.08*u.mrad, 0.20*u.mrad, + -0.08*u.mrad, 0.00*u.mrad, 0.00*u.mrad, + 0.15*u.mrad, 0.00*u.mrad, 0.00*u.mrad, + 0.20*u.mrad, 0.00*u.mrad, 0.00*u.mrad, + -0.20*u.mrad, 0.00*u.mrad, 0.00*u.mrad] - scifi_spatial_tags_2022 = ['t_0', 't_4361','t_4575','t_4855','t_5172','t_5431', 't_6443', 't_6677'] - scifi_spatial_tags_2024 = ['t_0', 't_7649', 't_8318', 't_8583', 't_8942', 't_9156', 't_9286', - 't_9379', 't_9462', 't_9613', 't_9692', 't_9882', 't_10012'] - if year == 2024: scifi_spatial_tags = scifi_spatial_tags_2024 - else: scifi_spatial_tags = scifi_spatial_tags_2022 + scifi_spatial_aligment_consts = { + 2022: ['t_0', 't_4361','t_4575','t_4855','t_5172'], + 2023: ['t_0', 't_5431', 't_6443', 't_6677'], + 2024: ['t_0', 't_7649', 't_8318', 't_8583', 't_8942', 't_9156', 't_9286', + 't_9379', 't_9462', 't_9613', 't_9692', 't_9882', 't_10012'], + 2025: ['t_0', 't_10423', 't_11158', 't_11576', 't_11676', 't_11795', 't_12021', 't_12123'], + 2026: ['t_0', 't_12794', 't_13118'] + } + scifi_spatial_tags = scifi_spatial_aligment_consts[year] for c in scifi_spatial_tags: k=0 for s in range(1,6): diff --git a/shipLHC/muonDis.py b/shipLHC/muonDis.py index 25401882fe..14e55b4580 100755 --- a/shipLHC/muonDis.py +++ b/shipLHC/muonDis.py @@ -181,7 +181,7 @@ def convertAscii2Root(fname,version=2): def muonPreTransport(): f = ROOT.TFile(options.muonIn) - nt = f.nt + nt = f.Get("nt") foutName = options.muonIn.replace('.root','_z'+str(options.z)+'.root') fout = ROOT.TFile(foutName,'recreate') variables = "" @@ -461,17 +461,23 @@ def makeMuDISEvents(withElossFunction=False,nucleon='p+'): myPythia.SetMRPY(1,R) mutype = {-13:'gamma/mu+',13:'gamma/mu-'} # DIS event +# run number +# event number # incoming muon, id:px:py:pz:x:y:z:w # outgoing particles, id:px:py:pz + run = array('i', [0]) + event = array('i', [0]) fout = ROOT.TFile('muonDis_'+str(options.run)+'.root','recreate') dTree = ROOT.TTree('DIS','muon DIS') + dTree.Branch("run", run, "run/I") + dTree.Branch("event", event, "event/I") iMuon = ROOT.TClonesArray("TParticle") iMuonBranch = dTree.Branch("InMuon",iMuon,32000,-1) dPart = ROOT.TClonesArray("TParticle") dPartBranch = dTree.Branch("Particles",dPart,32000,-1) # read file with muons hitting concrete wall fin = ROOT.TFile(options.muonIn) # id:px:py:pz:x:y:z:w - sTree = fin.nt + sTree = fin.Get("nt") nTOT = sTree.GetEntries() nEnd = min(nTOT,options.nStart + options.nEvents) # stop pythia printout during loop @@ -480,6 +486,8 @@ def makeMuDISEvents(withElossFunction=False,nucleon='p+'): nMade = 0 for k in range(options.nStart,nEnd): rc = sTree.GetEvent(k) + run[0] = int(sTree.run) + event[0] = int(sTree.event) # make n events / muon px,py,pz = sTree.px,sTree.py,sTree.pz x,y,z = sTree.x,sTree.y,sTree.z-SND_Z @@ -504,7 +512,9 @@ def makeMuDISEvents(withElossFunction=False,nucleon='p+'): for n in range(options.nMult): dPart.Clear() iMuon.Clear() - iMuon[0] = muPart + muPart_replica = ROOT.TParticle(muPart) + muPart_TCA = iMuon.ConstructedAt(0) + ROOT.std.swap(muPart_replica, muPart_TCA) myPythia.GenerateEvent() # remove all unnecessary stuff myPythia.Pyedit(2) @@ -518,7 +528,8 @@ def makeMuDISEvents(withElossFunction=False,nucleon='p+'): # copy to branch nPart = dPart.GetEntries() if dPart.GetSize() == nPart: dPart.Expand(nPart+10) - dPart[nPart] = part + part_TCA = dPart.ConstructedAt(nPart) + ROOT.std.swap(part, part_TCA) nMade+=1 if nMade%options.heartbeat==0: print('made so far ',options.run,' :',nMade,time.ctime()) dTree.Fill() @@ -832,7 +843,7 @@ def muonRateAtSND(withFaser=False,withEff=False,version=1): ut.bookHist(h,'tanThetaXYMS_'+str(mu), 'tan theta X/Y',200,-0.01,0.01,200,-0.01,0.01) ut.bookHist(h,'tanThetaXYMSfaser_'+str(mu), 'tan theta X/Y',200,-0.01,0.01,200,-0.01,0.01) ut.bookHist(h,'theta', 'mult scattering angle',200,-1.,1.) - for sTree in fin.nt: + for sTree in fin.Get("nt"): px,py,pz = sTree.px,sTree.py,sTree.pz m2cm = 1. if version == 0: m2cm = 100 @@ -1024,7 +1035,7 @@ def flukaMuons(version=1,Plimit=False,withFaser=True): ut.bookHist(h,'W_'+str(mu), 'w',100,0.,0.15) for mu in fnames: fin = ROOT.TFile(fnames[mu]) - for sTree in fin.nt: + for sTree in fin.Get("nt"): m2cm = 1. if version == 0: m2cm = 100. P = ROOT.TVector3(sTree.px,sTree.py,sTree.pz) @@ -1146,7 +1157,7 @@ def muInterGeant4(version=2,njobs=100): f = ROOT.TFile(fname) ROOT.gROOT.cd() nEv = -1 - for sTree in f.cbmsim: + for sTree in f.Get("cbmsim"): nEv+=1 muon = sTree.MCTrack[0] w = norm*muon.GetWeight() @@ -1214,7 +1225,7 @@ def muondEdX(version=2,njobs=100,path='',withFaser=False, plotOnly=True): h[fname] = ROOT.TFile(fname) ROOT.gROOT.cd() nEv = -1 - for sTree in h[fname].cbmsim: + for sTree in h[fname].Get("cbmsim"): nEv+=1 muon = sTree.MCTrack[0] w = norm*muon.GetWeight() @@ -1644,7 +1655,7 @@ def muonDISfull(cycle = 0, sMin=0,sMax=200,rMin=1,rMax=11,path = '/eos/experimen if path.find('eos')<0: h['f'] = ROOT.TFile.Open(path+prod+datafile) else: h['f'] = ROOT.TFile.Open(os.environ['EOSSHIP']+path+prod+datafile) nEv = -1 - for sTree in h['f'].cbmsim: + for sTree in h['f'].Get("cbmsim"): nEv+=1 if nEv%1000 == 0: print('N ',nEv,prod) muon = sTree.MCTrack[0] @@ -1787,7 +1798,7 @@ def thermNeutron(): ROOT.gROOT.cd() for pid in ['13','-13']: h['g_'+pid] = fin.Get('g_'+pid).Clone('g_'+pid) - for sTree in f.cbmsim: + for sTree in f.Get("cbmsim"): muon = sTree.MCTrack[0] muonEnergy = muon.GetEnergy() mupid = muon.GetPdgCode() diff --git a/shipLHC/rawData/ConvRawData.py b/shipLHC/rawData/ConvRawData.py index b01682e3bb..6c8f8ddf17 100644 --- a/shipLHC/rawData/ConvRawData.py +++ b/shipLHC/rawData/ConvRawData.py @@ -46,6 +46,7 @@ def Init(self,options): if options.path.find('2022')!=-1: fpath = "/eos/experiment/sndlhc/convertedData/physics/2022/" elif options.path.find('2023')!=-1: fpath = "/eos/experiment/sndlhc/convertedData/physics/2023/" elif options.path.find('2024')!=-1: fpath = "/eos/experiment/sndlhc/convertedData/physics/2024/" + elif options.path.find('2025')!=-1: fpath = "/eos/experiment/sndlhc/convertedData/physics/2025/" else: fpath = "/eos/experiment/sndlhc/convertedData/commissioning/TI18/" fg = ROOT.TFile.Open(options.server+fpath+"/FSdict.root") pkl = Unpickler(fg) @@ -93,11 +94,11 @@ def Init(self,options): if self.fiN.Get('event'): self.newFormat = False # old format if not self.monitoring: if self.newFormat: - if options.nEvents<0: self.nEvents = self.fiN.data.GetEntries() - else: self.nEvents = min(options.nEvents,self.fiN.data.GetEntries()) + if options.nEvents<0: self.nEvents = self.fiN.Get("data").GetEntries() + else: self.nEvents = min(options.nEvents,self.fiN.Get("data").GetEntries()) else: - if options.nEvents<0: self.nEvents = self.fiN.event.GetEntries() - else: self.nEvents = min(options.nEvents,self.fiN.event.GetEntries()) + if options.nEvents<0: self.nEvents = self.fiN.Get("event").GetEntries() + else: self.nEvents = min(options.nEvents,self.fiN.Get("event").GetEntries()) print('converting ',self.nEvents,' events ',' of run',options.runNumber) # Pass input parameters to the task - runN, paths, etc. ioman.RegisterInputObject('runN', ROOT.TObjString(str(options.runNumber))) @@ -429,7 +430,7 @@ def executeEvent1(self,eventNumber): if eventNumber%self.options.heartBeat==0 or self.debug: print('run ',self.options.runNumber, ' event',eventNumber," ",time.ctime()) - event = self.fiN.data + event = self.fiN.Get("data") event.GetEvent(eventNumber) self.header.SetEventTime(event.evt_timestamp) self.header.SetUTCtimestamp(int(event.evt_timestamp*6.23768*1e-9 + self.run_startUTC)) @@ -504,6 +505,8 @@ def executeEvent1(self,eventNumber): system = self.MufiSystem[board_id][tofpet_id] key = (tofpet_id%2)*1000 + tofpet_channel tmp = self.boardMaps['MuFilter'][board][self.slots[tofpet_id]] + # mini DTs are labelled DS 5V, board is 32 + if (tmp == "DS_5Vert"): continue if self.options.debug or not tmp.find('not')<0: print('debug',tmp,system,key,board,tofpet_id,tofpet_id%2,tofpet_channel) sipmChannel = 99 if not key in self.TofpetMap[system]: @@ -551,11 +554,15 @@ def executeEvent1(self,eventNumber): # copy hits to detector branches for sipmID in digiSciFiStore: if self.digiSciFi.GetSize() == indexSciFi: self.digiSciFi.Expand(indexSciFi+100) - self.digiSciFi[indexSciFi]=digiSciFiStore[sipmID] + aHit = digiSciFiStore[sipmID] + aHit_TCA = self.digiSciFi.ConstructedAt(indexSciFi) + ROOT.std.swap(aHit, aHit_TCA) indexSciFi+=1 for detID in digiMuFilterStore: if self.digiMuFilter.GetSize() == indexMuFilter: self.digiMuFilter.Expand(indexMuFilter+100) - self.digiMuFilter[indexMuFilter]=digiMuFilterStore[detID] + aHit = digiMuFilterStore[detID] + aHit_TCA = self.digiMuFilter.ConstructedAt(indexMuFilter) + ROOT.std.swap(aHit, aHit_TCA) indexMuFilter+=1 def executeEvent0(self,eventNumber): if self.options.FairTask_convRaw: @@ -568,7 +575,7 @@ def executeEvent0(self,eventNumber): return if eventNumber%self.options.heartBeat==0 or self.debug: print('run ',self.options.runNumber, ' event',eventNumber," ",time.ctime()) - event = self.fiN.event + event = self.fiN.Get("event") rc = event.GetEvent(eventNumber) self.header.SetEventTime(event.timestamp) self.header.SetRunId( self.options.runNumber ) @@ -635,6 +642,8 @@ def executeEvent0(self,eventNumber): system = self.MufiSystem[board_id][tofpet_id] key = (tofpet_id%2)*1000 + tofpet_channel tmp = self.boardMaps['MuFilter'][board][self.slots[tofpet_id]] + # mini DTs are labelled DS 5V, board is 32 + if (tmp == "DS_5Vert"): continue if self.options.debug or not tmp.find('not')<0: print('debug',tmp,system,key,board,tofpet_id,tofpet_id%2,tofpet_channel) sipmChannel = 99 if not key in self.TofpetMap[system]: @@ -682,11 +691,15 @@ def executeEvent0(self,eventNumber): # copy hits to detector branches for sipmID in digiSciFiStore: if self.digiSciFi.GetSize() == indexSciFi: self.digiSciFi.Expand(indexSciFi+100) - self.digiSciFi[indexSciFi]=digiSciFiStore[sipmID] + aHit = digiSciFiStore[sipmID] + aHit_TCA = self.digiSciFi.ConstructedAt(indexSciFi) + ROOT.std.swap(aHit, aHit_TCA) indexSciFi+=1 for detID in digiMuFilterStore: if self.digiMuFilter.GetSize() == indexMuFilter: self.digiMuFilter.Expand(indexMuFilter+100) - self.digiMuFilter[indexMuFilter]=digiMuFilterStore[detID] + aHit = digiMuFilterStore[detID] + aHit_TCA = self.digiMuFilter.ConstructedAt(indexMuFilter) + ROOT.std.swap(aHit, aHit_TCA) indexMuFilter+=1 @@ -701,7 +714,7 @@ def Finalize(self): self.outfile.Close() else: if self.options.debug: - print('number of events processed',self.sTree.GetEntries(),self.fiN.event.GetEntries()) + print('number of events processed',self.sTree.GetEntries(),self.fiN.Get("event").GetEntries()) self.sTree.Write() self.fiN.Close() self.fSink.Close() diff --git a/shipLHC/rawData/convertRawData_muTestbeam.py b/shipLHC/rawData/convertRawData_muTestbeam.py index 857c72c083..9e6a20cd7f 100644 --- a/shipLHC/rawData/convertRawData_muTestbeam.py +++ b/shipLHC/rawData/convertRawData_muTestbeam.py @@ -269,8 +269,8 @@ def channel(tofpet_id,tofpet_channel,position): if not (options.partition < 0): part = '_'+str(options.partition).zfill(4) dataName = 'data'+part+'.root' f0=ROOT.TFile.Open(X+path+dataName) -if options.nEvents<0: nEvent = f0.event.GetEntries() -else: nEvent = min(options.nEvents,f0.event.GetEntries()) +if options.nEvents<0: nEvent = f0.Get("event").GetEntries() +else: nEvent = min(options.nEvents,f0.Get("event").GetEntries()) print('converting ',nEvent,' events ',' of run',options.runNumber) boards = {} @@ -303,7 +303,7 @@ def channel(tofpet_id,tofpet_channel,position): import time def run(nEvent): - event = f0.event + event = f0.Get("event") for eventNumber in range(options.nStart,nEvent): ncreated = 0 rc = event.GetEvent(eventNumber) @@ -471,11 +471,15 @@ def run(nEvent): for sipmID in digiSciFiStore: if digiSciFi.GetSize() == indexSciFi: digiSciFi.Expand(indexSciFi+100) - digiSciFi[indexSciFi]=digiSciFiStore[sipmID] + aHit = digiSciFiStore[sipmID] + aHit_TCA = digiSciFi.ConstructedAt(indexSciFi) + ROOT.std.swap(aHit, aHit_TCA) indexSciFi+=1 for detID in digiMuFilterStore: if digiMuFilter.GetSize() == indexMuFilter: digiMuFilter.Expand(indexMuFilter+100) - digiMuFilter[indexMuFilter]=digiMuFilterStore[detID] + aHit = digiMuFilterStore[detID] + aHit_TCA = digiMuFilter.ConstructedAt(indexMuFilter) + ROOT.std.swap(aHit, aHit_TCA) indexMuFilter+=1 # make simple clustering for scifi, only works with geometry file. Don't think it is a good idea at the moment if withGeoFile: @@ -509,7 +513,8 @@ def run(nEvent): aCluster = ROOT.sndCluster(first,N,hitlist,scifiDet) print("cluster created") if clusScifi.GetSize() == index: clusScifi.Expand(index+10) - clusScifi[index]=aCluster + aCluster_TCA = clusScifi.ConstructedAt(index) + ROOT.std.swap(aCluster, aCluster_TCA) index+=1 if c!=hitList[last]: ncl+=1 @@ -520,7 +525,7 @@ def run(nEvent): if options.debug: print('====> number of created hits',ncreated) if options.debug: - print('number of events processed',sTree.GetEntries(),f0.event.GetEntries()) + print('number of events processed',sTree.GetEntries(),f0.Get("event").GetEntries()) sTree.Write() # https://gitlab.cern.ch/snd-scifi/software/-/wikis/Raw-data-format @@ -533,7 +538,7 @@ def run(nEvent): theMap = {} def getMapEvent2Time(): k = 0 - for event in f0.event: + for event in f0.Get("event"): theMap[event.timestamp]=k k+=1 if k%10000000==0: print('key time map: now at ',k) diff --git a/shipLHC/rawData/mufiHitMaps.py b/shipLHC/rawData/mufiHitMaps.py index 52df158965..7b3700ed73 100644 --- a/shipLHC/rawData/mufiHitMaps.py +++ b/shipLHC/rawData/mufiHitMaps.py @@ -45,10 +45,10 @@ def smallSiPMchannel(i): if options.runNumber>0: f=ROOT.TFile.Open(path+'sndsw_raw_'+str(options.runNumber).zfill(6)+'.root') - eventTree = f.rawConv + eventTree = f.Get("rawConv") else: f=ROOT.TFile.Open(options.fname) - eventTree = f.cbmsim + eventTree = f.Get("cbmsim") # backward compatbility for early converted events eventTree.GetEvent(0) diff --git a/shipLHC/rawData/runProd.py b/shipLHC/rawData/runProd.py index d826235a1e..aeb037fdc2 100644 --- a/shipLHC/rawData/runProd.py +++ b/shipLHC/rawData/runProd.py @@ -209,15 +209,15 @@ def checkFile(self,path,r,p): print('checkfile',path,r,p) inFile = self.options.server+path+'run_'+ r+'/data_'+p+'.root' fI = ROOT.TFile.Open(inFile) - if fI.Get('event'): Nraw = fI.event.GetEntries() - else: Nraw = fI.data.GetEntries() + if fI.Get('event'): Nraw = fI.Get("event").GetEntries() + else: Nraw = fI.Get("data").GetEntries() outFile = 'sndsw_raw_'+r+'-'+p+self.Mext+'.root' fC = ROOT.TFile(outFile) test = fC.Get('rawConv') if not test: print('Error:',path,r,p,' rawConv not found') return -2 - Nconv = fC.rawConv.GetEntries() + Nconv = fC.Get("rawConv").GetEntries() if Nraw != Nconv: print('Error:',path,r,p,':',Nraw,Nconv) return -1 @@ -331,7 +331,7 @@ def getConvStats(self,runList): for run in runList: try: f=ROOT.TFile.Open("sndsw_raw_"+str(run).zfill(6)+'.root') - print(run,':',f.rawConv.GetEntries()) + print(run,':',f.Get("rawConv").GetEntries()) except: print('problem:',run) @@ -339,8 +339,8 @@ def rawStats(self,runList): for run in runList: runNr = str(run).zfill(6) r = ROOT.TFile.Open(os.environ['EOSSHIP']+path+"/run_"+runNr+"/data.root") - if fI.Get('event'): raw = r.event.GetEntries() - else: raw = r.data.GetEntries() + if fI.Get('event'): raw = r.Get("event").GetEntries() + else: raw = r.Get("data").GetEntries() print(run,':',raw) def makeHistos(self,runList): diff --git a/shipLHC/rawData/scifiHitMaps.py b/shipLHC/rawData/scifiHitMaps.py index 61a875ca44..fbb510af2a 100644 --- a/shipLHC/rawData/scifiHitMaps.py +++ b/shipLHC/rawData/scifiHitMaps.py @@ -32,10 +32,10 @@ def xPos(detID): if options.runNumber>0: f=ROOT.TFile.Open(options.path+'sndsw_raw_'+str(options.runNumber).zfill(6)+'.root') - eventTree = f.rawConv + eventTree = f.Get("rawConv") else: f=ROOT.TFile.Open(options.fname) - eventTree = f.cbmsim + eventTree = f.Get("cbmsim") def slopes(Nev=-1): A,B = ROOT.TVector3(),ROOT.TVector3() diff --git a/shipLHC/run_MCEventBuilder.py b/shipLHC/run_MCEventBuilder.py new file mode 100644 index 0000000000..3c7d67efaf --- /dev/null +++ b/shipLHC/run_MCEventBuilder.py @@ -0,0 +1,62 @@ +import os +import ROOT +from ROOT import TObjString +import time +from argparse import ArgumentParser +import SndlhcGeo +import shipunit as u + +# ------------ Argument parser ---------------- +parser = ArgumentParser() +parser.add_argument("-f", "--inputFile", help="Input file", required=True) +parser.add_argument("-g", "--geoFile", help="geo file", required=True) +parser.add_argument("-o", "--outputFile", help="Output file", required=True) +parser.add_argument("-i", "--firstEvent", type=int, default=0, help="First event to process") +parser.add_argument("-n", "--nEvents", type=int, default=0, help="Number of events to process (0 = all)") +parser.add_argument("--saveFirst25nsOnly", action="store_true", help="Only store the first 25-ns chunk of each event") +options = parser.parse_args() + +# ------------ Geo setup ---------------- +geo = SndlhcGeo.GeoInterface(options.geoFile) +lsOfGlobals = ROOT.gROOT.GetListOfGlobals() +lsOfGlobals.Add(geo.modules['Scifi']) +scifiDet = lsOfGlobals.FindObject('Scifi') +scifiDet.SetConfPar("Scifi/signalSpeed", 15*u.cm/u.nanosecond) +lsOfGlobals.Add(geo.modules['MuFilter']) + + +#-----------Executioner-------------- +start = time.time() +inRootTFile = ROOT.TFile(options.inputFile) +print(f"Input file: {options.inputFile}") + +# Use FairRoot framework to arrange the workflow +# A FairRun is a wrapper of a collection of tasks +run = ROOT.FairRunAna() + +# Input/output manager +ioman = ROOT.FairRootManager.Instance() +source = ROOT.FairFileSource(inRootTFile) +run.SetSource(source) +outFile = ROOT.TMemFile('dummy','CREATE') #IGNORE +sink = ROOT.FairRootFileSink(outFile) +run.SetSink(sink) +ioman.InitSink() +run.SetEventHeaderPersistence(False) + +#Avoiding some error messages +xrdb = ROOT.FairRuntimeDb.instance() +xrdb.getContainer("FairBaseParSet").setStatic() +xrdb.getContainer("FairGeoParSet").setStatic() + +# Add tasks +eventBuilder = ROOT.MCEventBuilder(options.outputFile, options.saveFirst25nsOnly) +run.AddTask(eventBuilder) + +# Initialize and run +run.Init() +run.Run(options.firstEvent, options.firstEvent+options.nEvents) + +end = time.time() +elapsed = end - start +print(f"Elapsed time: {elapsed:.2f} seconds") diff --git a/shipLHC/run_digiSND.py b/shipLHC/run_digiSND.py index 91efc660e3..04a9e26e64 100644 --- a/shipLHC/run_digiSND.py +++ b/shipLHC/run_digiSND.py @@ -7,7 +7,6 @@ def pyExit(): # This is needed to bypass seg violation with exiting cpp digitization # Most likely related to file ownership. os.system('kill '+str(os.getpid())) -atexit.register(pyExit) import resource def mem_monitor(): @@ -72,19 +71,48 @@ def mem_monitor(): lsOfGlobals = ROOT.gROOT.GetListOfGlobals() scifiDet = lsOfGlobals.FindObject('Scifi') mufiDet = lsOfGlobals.FindObject('MuFilter') -mufiDet.SetConfPar("MuFilter/DsAttenuationLength",350 * u.cm) # values between 300 cm and 400cm observed for H6 testbeam -mufiDet.SetConfPar("MuFilter/DsTAttenuationLength",700 * u.cm) # top readout with mirror on bottom -mufiDet.SetConfPar("MuFilter/VandUpAttenuationLength",999 * u.cm) # no significante attenuation observed for H6 testbeam -mufiDet.SetConfPar("MuFilter/DsSiPMcalibrationS",25.*1000.) # in MC: 1.65 keV are about 41.2 qdc -mufiDet.SetConfPar("MuFilter/VandUpSiPMcalibration",25.*1000.); -mufiDet.SetConfPar("MuFilter/VandUpSiPMcalibrationS",25.*1000.); -mufiDet.SetConfPar("MuFilter/VandUpPropSpeed",12.5*u.cm/u.nanosecond); -mufiDet.SetConfPar("MuFilter/DsPropSpeed",14.3*u.cm/u.nanosecond); scifiDet.SetConfPar("Scifi/nphe_min",options.ts) # threshold scifiDet.SetConfPar("Scifi/nphe_max",options.ss) # saturation -scifiDet.SetConfPar("Scifi/timeResol",150.*u.picosecond) # time resolution in ps -scifiDet.SetConfPar("MuFilter/timeResol",150.*u.picosecond) # time resolution in ps, first guess +#### +# The lines below aim to reproduce the original digitization case, but urge the user to regenerate +# the sample shall it be outdated from before the removal of the 1MeV production cut, which +# coincides with MuFi digi const update. +# +# in MC productions generated before July 2022 Scifi signal speed is missing from the geofile +better_update = False +if scifiDet.GetConfParF("Scifi/signalSpeed")==0: + scifiDet.SetConfPar("Scifi/signalSpeed", 15*u.cm/u.nanosecond) + better_update = True +# geofiles before March 2026 don't have the Veto 3 atten.length +if mufiDet.GetConfParF("MuFilter/VTAttenuationLength")==0: + mufiDet.SetConfPar("MuFilter/VTAttenuationLength",999*u.cm) + better_update = True +# old digi constants from before the MuFi response update +if mufiDet.GetConfParF("MuFilter/VandUpAttenuationLength")==999*u.cm: + better_update = True +# in very ancient MC productions it is possible some digitization params are missing +# set them here values updated in March 2026. +if mufiDet.GetConfParF("MuFilter/DsAttenuationLength")==0 or\ + mufiDet.GetConfParF("MuFilter/VandUpPropSpeed")==0 : + mufiDet.SetConfPar("MuFilter/DsAttenuationLength",230*u.cm) + mufiDet.SetConfPar("MuFilter/DsTAttenuationLength",700*u.cm) + mufiDet.SetConfPar("MuFilter/VandUpAttenuationLength",210*u.cm) + mufiDet.SetConfPar("MuFilter/VTAttenuationLength",999*u.cm) + mufiDet.SetConfPar("MuFilter/DsSiPMcalibration",25.*1000.) + # 1.65 MeV = 41 qcd over 6 Large SiPMs(one side) + mufiDet.SetConfPar("MuFilter/VandUpSiPMcalibrationL",50.*1000.) + # no MIP signal for small SiPMs, delayed and compromised response in general + mufiDet.SetConfPar("MuFilter/VandUpSiPMcalibrationS",0.) + mufiDet.SetConfPar("MuFilter/VandUpPropSpeed",13.6*u.cm/u.nanosecond); + mufiDet.SetConfPar("MuFilter/DsPropSpeed",15.1*u.cm/u.nanosecond); + scifiDet.SetConfPar("Scifi/timeResol",150.*u.picosecond) + mufiDet.SetConfPar("MuFilter/timeResol",150.*u.picosecond) # time resolution in ps, first guess + better_update = True + +if better_update: + print("WARNING: Simulation file preceding the production cut change! Consider regenerating from scratch!") +#### # Fair digitization task if options.FairTask_digi: diff --git a/shipLHC/run_muonRecoSND.py b/shipLHC/run_muonRecoSND.py index 9692868a60..c3a990f265 100644 --- a/shipLHC/run_muonRecoSND.py +++ b/shipLHC/run_muonRecoSND.py @@ -75,7 +75,7 @@ treename = None for test_treename in ["rawConv", "cbmsim"] : - if hasattr(F, test_treename) : + if F.Get(test_treename): treename = test_treename break diff --git a/shipLHC/run_simSND.py b/shipLHC/run_simSND.py index 85d5e19778..824c12322a 100644 --- a/shipLHC/run_simSND.py +++ b/shipLHC/run_simSND.py @@ -32,7 +32,9 @@ parser.add_argument("--pID", dest="pID", help="id of particle used by the gun (default=22)", required=False, default=22, type=int) parser.add_argument("--Estart", dest="Estart", help="start of energy range of particle gun for muflux detector (default=10 GeV)", required=False, default=10, type=float) parser.add_argument("--Eend", dest="Eend", help="end of energy range of particle gun for muflux detector (default=10 GeV)", required=False, default=10, type=float) +parser.add_argument('--eMin_store', type=float, help="kinetic energy cut for !particle storage! in MeV", dest='emin_store', default=100.) +parser.add_argument("--PGrunID", dest="PGrunID",help="PG run ID", required=False, type=int) parser.add_argument("--multiplePGSources", help="Multiple particle guns in a x-y plane at a fixed z or in a 3D volume", action="store_true") parser.add_argument("--EVx", dest="EVx", help="particle gun start xpos", required=False, default=0, type=float) parser.add_argument("--EVy", dest="EVy", help="particle gun start ypos", required=False, default=0, type=float) @@ -173,6 +175,9 @@ def Exec(self,opt): # -----Particle Gun----------------------- if simEngine == "PG": + if not options.PGrunID: + print("Missing option '--PGrunID', which provides PG run ID. Set it and run again!") + exit() myPgun = ROOT.FairBoxGenerator(options.pID,1) myPgun.SetPRange(options.Estart,options.Eend) myPgun.SetPhiRange(0, 360) # // Azimuth angle range [degree] @@ -197,6 +202,7 @@ def Exec(self,opt): f'({options.EVx},{options.EVy},{options.EVz})[cm × cm × cm] \n' f'with a uniform x-y spread of (Dx,Dy)=({options.Dx},{options.Dy})[cm × cm]' f' and {options.nZSlices} z slices in steps of {options.zSliceStep}[cm].') + run.SetPythiaDecayer('DecayConfigPy8.C') ROOT.FairLogger.GetLogger().SetLogScreenLevel("WARNING") # otherwise stupid printout for each event # -----muon DIS Background------------------------ if simEngine == "muonDIS": @@ -209,6 +215,7 @@ def Exec(self,opt): primGen.AddGenerator(DISgen) options.nEvents = min(options.nEvents,DISgen.GetNevents()) inactivateMuonProcesses = True # avoid unwanted hadronic events of "incoming" muon flying backward + run.SetPythiaDecayer('DecayConfigPy8.C') print('MuDIS position info input=',mu_start, mu_end) print('Generate ',options.nEvents,' with DIS input', ' first event',options.firstEvent) @@ -257,6 +264,7 @@ def Exec(self,opt): Ntuplegen.Init(inputFile,options.firstEvent) primGen.AddGenerator(Ntuplegen) options.nEvents = min(options.nEvents,Ntuplegen.GetNevents()) + run.SetPythiaDecayer('DecayConfigPy8.C') if simEngine == "MuonBack": # reading muon tracks from FLUKA @@ -305,6 +313,12 @@ def Exec(self,opt): # -----Initialize simulation run------------------------------------ run.Init() + +if simEngine == "PG": + # set the runID for + theHeader = run.GetMCEventHeader() + theHeader.SetRunID(options.PGrunID) + gMC = ROOT.TVirtualMC.GetMC() fStack = gMC.GetStack() if MCTracksWithHitsOnly: @@ -312,10 +326,10 @@ def Exec(self,opt): fStack.SetEnergyCut(-100.*u.MeV) elif MCTracksWithEnergyCutOnly: fStack.SetMinPoints(-1) - fStack.SetEnergyCut(100.*u.MeV) + fStack.SetEnergyCut(options.emin_store*u.MeV) elif MCTracksWithHitsOrEnergyCut: fStack.SetMinPoints(1) - fStack.SetEnergyCut(100.*u.MeV) + fStack.SetEnergyCut(options.emin_store*u.MeV) elif options.deepCopy: fStack.SetMinPoints(0) fStack.SetEnergyCut(0.*u.MeV) @@ -404,7 +418,7 @@ def Exec(self,opt): if options.genie == 4 : f_input = ROOT.TFile(inputFile) - gst = f_input.gst + gst = f_input.Get("gst") selection_string = "(Entry$ >= "+str(options.firstEvent)+")" if (options.firstEvent + options.nEvents) < gst.GetEntries() : @@ -433,6 +447,7 @@ def Exec(self,opt): # ------------------------------------------------------------------------ def checkOverlaps(removeScifi=False): + ROOT.gROOT.SetWebDisplay("off") # Workaround for https://github.com/root-project/root/issues/18881 sGeo = ROOT.gGeoManager if removeScifi: for n in range(1,6): diff --git a/shipLHC/scripts/2dEventDisplay.py b/shipLHC/scripts/2dEventDisplay.py index 204804064f..300aad33b9 100644 --- a/shipLHC/scripts/2dEventDisplay.py +++ b/shipLHC/scripts/2dEventDisplay.py @@ -14,7 +14,12 @@ from datetime import datetime from pathlib import Path +import logging +logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.WARNING) +logging.basicConfig(format='%(levelname)s: %(message)s', level=logging.FATAL) + ROOT.gStyle.SetPalette(ROOT.kViridis) +ROOT.gInterpreter.ProcessLine('#include "'+os.environ['SNDSW_ROOT']+'/analysis/tools/sndSciFiTools.h"') def pyExit(): "unfortunately need as bypassing an issue related to use xrootd" @@ -35,12 +40,26 @@ def pyExit(): parser.add_argument("-g", "--geoFile", dest="geoFile", help="geofile", default=os.environ["EOSSHIP"]+"/eos/experiment/sndlhc/convertedData/physics/2022/geofile_sndlhc_TI18_V0_2022.root") parser.add_argument("-P", "--partition", dest="partition", help="partition of data", type=int,required=False,default=-1) parser.add_argument("--server", dest="server", help="xrootd server",default=os.environ["EOSSHIP"]) -parser.add_argument("-X", dest="extraInfo", help="print extra event info",default=True) +parser.add_argument("--no-extraInfo", dest="extraInfo", action="store_false", help="Do not print extra information on the event.") + +parser.add_argument("--extension", help="Extension of file to save. E.g. png, pdf, root, etc.", default="png") +parser.add_argument("--rootbatch", help="Run ROOT in batch mode.", action="store_true") + +parser.add_argument("--collision_axis", dest="drawCollAxis", help="Draw collision axis", action="store_true") parser.add_argument("-par", "--parFile", dest="parFile", help="parameter file", default=os.environ['SNDSW_ROOT']+"/python/TrackingParams.xml") parser.add_argument("-hf", "--HoughSpaceFormat", dest="HspaceFormat", help="Hough space representation. Should match the 'Hough_space_format' name in parFile, use quotes", default='linearSlopeIntercept') +parser.add_argument("--shower_dir", dest="drawShowerDir", help="Draw shower direction if a shower is found", action="store_true") + options = parser.parse_args() + +resolution_factor = 1 +if options.rootbatch: + ROOT.gROOT.SetBatch() + # Produce figures larger than the screen resolution. E.g., for printing. + resolution_factor = 2 + options.storePic = '' trans2local = False runInfo = False @@ -89,11 +108,11 @@ def pyExit(): f=ROOT.TFile.Open(options.path+options.inputFile) if f.FindKey('cbmsim'): - eventTree = f.cbmsim + eventTree = f.Get("cbmsim") runId = 'sim' if eventTree.GetBranch('ScifiPoint'): mc = True else: - eventTree = f.rawConv + eventTree = f.Get("rawConv") ioman.SetTreeName('rawConv') outFile = ROOT.TMemFile('dummy','CREATE') @@ -242,7 +261,7 @@ def getSciFiHitDensity(g, x_range=0.5): ret.append(density) return ret -def drawLegend(max_density, max_QDC, n_legend_points): +def drawLegend(max_density, max_QDC, n_legend_points, darkMode=False): """Draws legend for hit colour""" h['simpleDisplay'].cd(1) n_legend_points = 5 @@ -254,6 +273,7 @@ def drawLegend(max_density, max_QDC, n_legend_points): text_scifi_legend.SetTextAlign(11) text_scifi_legend.SetTextFont(42) text_scifi_legend.SetTextSize(.15) + if darkMode: text_scifi_legend.SetTextColor(ROOT.kWhite) for i in range(n_legend_points) : if i < (n_legend_points - 1) : text_scifi_legend.DrawLatex((i+0.3)*(1./(n_legend_points+2)), 0.2, "{:d}".format(int(i*max_density/(n_legend_points-1)))) @@ -263,11 +283,17 @@ def drawLegend(max_density, max_QDC, n_legend_points): text_scifi_legend.DrawLatex((i+0.3)*(1./(n_legend_points+2)), 0., "{:.0f} QDC units".format(int(i*max_QDC/(n_legend_points-1)))) h["markerCollection"].append(ROOT.TEllipse((i+0.15)*(1./(n_legend_points+2)), 0.26, 0.05/4, 0.05)) - h["markerCollection"][-1].SetFillColor(ROOT.TColor.GetPalette()[int(float(i*max_density/(n_legend_points-1))/max_density*(len(ROOT.TColor.GetPalette())-1))]) + if not darkMode: + h["markerCollection"][-1].SetFillColor(ROOT.TColor.GetPalette()[int(float(i*max_density/(n_legend_points-1))/max_density*(len(ROOT.TColor.GetPalette())-1))]) + else: + h["markerCollection"][-1].SetFillColor(ROOT.TColor.GetPalette()[int((1/4+3/4*float(i*max_density/(n_legend_points-1))/max_density)*(len(ROOT.TColor.GetPalette())-1))]) h["markerCollection"][-1].Draw("SAME") h["markerCollection"].append(ROOT.TBox((i+0.15)*(1./(n_legend_points+2))-0.05/4 , 0.06 - 0.05, (i+0.15)*(1./(n_legend_points+2))+0.05/4, 0.06 + 0.05)) - h["markerCollection"][-1].SetFillColor(ROOT.TColor.GetPalette()[int(float(i*max_QDC/(n_legend_points-1))/max_QDC*(len(ROOT.TColor.GetPalette())-1))]) + if not darkMode: + h["markerCollection"][-1].SetFillColor(ROOT.TColor.GetPalette()[int(float(i*max_QDC/(n_legend_points-1))/max_QDC*(len(ROOT.TColor.GetPalette())-1))]) + else: + h["markerCollection"][-1].SetFillColor(ROOT.TColor.GetPalette()[int((1/4+3/4*float(i*max_QDC/(n_legend_points-1))/max_QDC)*(len(ROOT.TColor.GetPalette())-1))]) h["markerCollection"][-1].Draw("SAME") def drawSciFiHits(g, colour): @@ -296,12 +322,38 @@ def loopEvents( minSipmMult=1, withTiming=False, option=None, - Setup='', + Setup='TI18', verbose=0, auto=False, - hitColour=None + hitColour=None, + FilterScifiHits=None, + darkMode=False ): - if 'simpleDisplay' not in h: ut.bookCanvas(h,key='simpleDisplay',title='simple event display',nx=1200,ny=1600,cx=1,cy=2) + if darkMode: + ROOT.gStyle.SetCanvasColor(ROOT.kBlack) + + # check the format of FilterScifiHits if set + if FilterScifiHits: + important_keys = {"bins_x", "min_x", "max_x", "time_lower_range", "time_upper_range"} + all_keys = important_keys.copy() + all_keys.add("method") + filter_parameters = {"bins_x":52., "min_x":0., "max_x":26., + "time_lower_range":1E9/(2*u.snd_freq/u.hertz), + "time_upper_range":2E9/(u.snd_freq/u.hertz), + "method":0} + if FilterScifiHits!="default" and not important_keys.issubset(FilterScifiHits): + logging.fatal("Invalid FilterScifiHits format. Two options are supported:\n" + "#1 FilterScifiHits = 'default'\nwhich sets the default parameters:\n"+ + str(filter_parameters)+" or\n" + "#2 FilterScifiHits = filter_dictionary \nwhere filter_dictionary has all of the following keys\n"+ + str(important_keys)+"\nAn additional key 'method' exists: its single supported value, also default, is 0.") + return + if FilterScifiHits!="default" and any(k not in all_keys for k in FilterScifiHits): + logging.warning("Ignoring provided keys other than "+str(all_keys)) + + if 'simpleDisplay' not in h: + ut.bookCanvas(h,key='simpleDisplay',title='simple event display',nx=1200,ny=1016,cx=1,cy=2) + h['simpleDisplay'].cd(1) # TI18 coordinate system zStart = 250. @@ -325,6 +377,15 @@ def loopEvents( for d in ['xmin','xmax','ymin','ymax','zmin','zmax']: h['c'+d]=h[d] ut.bookHist(h,'xz','; z [cm]; x [cm]',500,h['czmin'],h['czmax'],100,h['cxmin'],h['cxmax']) ut.bookHist(h,'yz','; z [cm]; y [cm]',500,h['czmin'],h['czmax'],100,h['cymin'],h['cymax']) + if darkMode: + for hist in [h['xz'], h['yz']]: + # set axis lines, labels, titles to white in dark mode + hist.GetXaxis().SetAxisColor(ROOT.kWhite) + hist.GetYaxis().SetAxisColor(ROOT.kWhite) + hist.GetXaxis().SetLabelColor(ROOT.kWhite) + hist.GetYaxis().SetLabelColor(ROOT.kWhite) + hist.GetXaxis().SetTitleColor(ROOT.kWhite) + hist.GetYaxis().SetTitleColor(ROOT.kWhite) proj = {1:'xz',2:'yz'} h['xz'].SetStats(0) @@ -403,7 +464,42 @@ def loopEvents( if nAlltracks > 0: print('total number of tracks: ', nAlltracks) digis = [] - if event.FindBranch("Digi_ScifiHits"): digis.append(event.Digi_ScifiHits) + if event.FindBranch("Digi_ScifiHits"): + scifi_digis = event.Digi_ScifiHits + method = 0 + if FilterScifiHits!=None and FilterScifiHits!="default": + filter_parameters = {k: FilterScifiHits[k] for k in important_keys if k in FilterScifiHits} + method = FilterScifiHits.get("method", 0) # set to the default 0, if item is not provided + if FilterScifiHits and (Setup=="TI18" or Setup=="H8" or Setup=="H4"): + setup = Setup + # Only H8 is explicitly supported in the SciFi tools. However, the same baby SciFi + # system was reused in H4. It is then safe to use the SciFi tools for H4 as well. + if Setup =="H4": + setup ="H8" + # Convert the filter_parameters to the needed std.map format + selection_parameters = ROOT.std.map('string', 'float')() + selection_parameters["bins_x"] = float(filter_parameters["bins_x"]) + selection_parameters["min_x"] = float(filter_parameters["min_x"]) + selection_parameters["max_x"] = float(filter_parameters["max_x"]) + selection_parameters["time_lower_range"] = float(filter_parameters["time_lower_range"]) + selection_parameters["time_upper_range"] = float(filter_parameters["time_upper_range"]) + scifi_digis = ROOT.snd.analysis_tools.filterScifiHits(event.Digi_ScifiHits,selection_parameters,method,setup,mc) + else: + if FilterScifiHits: + logging.warning(Setup+" is not supported for the time-filtering of SciFi hits, using all hits instead.") + digis.append(scifi_digis) + run_conf = ROOT.snd.Configuration.Option.ti18_2022_2023 + if Setup=='TI18' or Setup=="H6": + run_conf = ROOT.snd.Configuration.Option.ti18_2022_2023 + elif Setup == "H8": + ROOT.snd.Configuration.Option.test_beam_2023 + elif Setup == "H4": + ROOT.snd.Configuration.Option.test_beam_2024 + else: + print("Going for default setup TI18") + configuration = ROOT.snd.Configuration(run_conf, geo.modules['Scifi'], geo.modules['MuFilter']) + scifi_planes = ROOT.snd.analysis_tools.FillScifi(configuration, scifi_digis, geo.modules['Scifi']) + us_planes = ROOT.snd.analysis_tools.FillUS(configuration, event.Digi_MuFilterHits, geo.modules['MuFilter'], mc) if event.FindBranch("Digi_MuFilterHits"): digis.append(event.Digi_MuFilterHits) if event.FindBranch("Digi_MuFilterHit"): digis.append(event.Digi_MuFilterHit) empty = True @@ -440,8 +536,12 @@ def loopEvents( rc = h[ 'simpleDisplay'].cd(p) h[proj[p]].Draw('b') + if options.drawCollAxis: + for k in proj: + drawCollisionAxis(h['simpleDisplay'], k) + if withDetector: - drawDetectors() + drawDetectors(darkMode=darkMode) for D in digis: for digi in D: detID = digi.GetDetectorID() @@ -467,10 +567,16 @@ def loopEvents( nav.MasterToLocal(globA, locA) Z = X[2] if digi.isVertical(): + # only using hits with positive qdc for centroids, so only show such + if options.drawShowerDir and system==0 and digi.GetSignal(0)<0: + continue collection = 'hitCollectionX' Y = locA[0] sY = detSize[system][0] else: + # only using hits with positive qdc for centroids, so only show such + if options.drawShowerDir and system==0 and digi.GetSignal(0)<0: + continue collection = 'hitCollectionY' Y = locA[1] sY = detSize[system][1] @@ -489,9 +595,12 @@ def loopEvents( this_qdc += qdc if this_qdc > max_QDC : this_qdc = max_QDC - fillNode(curPath, ROOT.TColor.GetPalette()[int(this_qdc/max_QDC*(len(ROOT.TColor.GetPalette())-1))]) + if not darkMode: + fillNode(curPath, ROOT.TColor.GetPalette()[int(this_qdc/max_QDC*(len(ROOT.TColor.GetPalette())-1))]) + else: # for dark mode, only use the lighter part of the colormap; the dark hues have bad contrast + fillNode(curPath, ROOT.TColor.GetPalette()[int((1/4 + 3/4*this_qdc/max_QDC) * (len(ROOT.TColor.GetPalette())-1))]) else : - fillNode(curPath) + fillNode(curPath, darkMode=darkMode) if digi.isVertical(): F = 'firedChannelsX' else: F = 'firedChannelsY' @@ -504,17 +613,20 @@ def loopEvents( h[F][systems[system]][0]+=1 if len(h[F][systems[system]]) < 2+side: continue h[F][systems[system]][2+side]+=qdc - h['hitCollectionY']['Scifi'][1].SetMarkerColor(ROOT.kBlue+2) - h['hitCollectionX']['Scifi'][1].SetMarkerColor(ROOT.kBlue+2) + h['hitCollectionY']['Scifi'][1].SetMarkerColor(ROOT.kBlue+2 if not darkMode else ROOT.kBlue-4) + h['hitCollectionX']['Scifi'][1].SetMarkerColor(ROOT.kBlue+2 if not darkMode else ROOT.kBlue-4) if hitColour == "q" : for orientation in ['X', 'Y']: max_density = 40 density = np.clip(0, max_density, getSciFiHitDensity(h['hitCollection'+orientation]['Scifi'][1])) for i in range(h['hitCollection'+orientation]['Scifi'][1].GetN()) : - h['hitColour'+orientation]['Scifi'].append(ROOT.TColor.GetPalette()[int(float(density[i])/max_density*(len(ROOT.TColor.GetPalette())-1))]) + if not darkMode: + h['hitColour'+orientation]['Scifi'].append(ROOT.TColor.GetPalette()[int(float(density[i])/max_density*(len(ROOT.TColor.GetPalette())-1))]) + else: # for dark mode, only use the lighter part of the colormap; the dark hues have bad contrast + h['hitColour'+orientation]['Scifi'].append(ROOT.TColor.GetPalette()[int((1/4 + 3/4*float(density[i])/max_density) * (len(ROOT.TColor.GetPalette())-1))]) - drawLegend(max_density, max_QDC, 5) + drawLegend(max_density, max_QDC, 5, darkMode=darkMode) k = 1 moreEventInfo = [] @@ -522,7 +634,7 @@ def loopEvents( for collection in ['hitCollectionX','hitCollectionY']: h['simpleDisplay'].cd(k) - drawInfo(h['simpleDisplay'], k, runId, N, T) + drawInfo(h['simpleDisplay'], k, runId, N, T, darkMode=darkMode) k+=1 for c in h[collection]: F = collection.replace('hitCollection','firedChannels') @@ -558,49 +670,76 @@ def loopEvents( k = 1 for collection in ['hitCollectionX','hitCollectionY']: h['simpleDisplay'].cd(k) - drawInfo(h['simpleDisplay'], k, runId, N, T,moreEventInfo) + drawInfo(h['simpleDisplay'], k, runId, N, T,moreEventInfo, darkMode=darkMode) k+=1 h['simpleDisplay'].Update() if withTiming: timingOfEvent() - addTrack(OT) + addTrack(OT, darkMode=darkMode) + + # try finding shower direction and intercept + if options.drawShowerDir: + sh_scifi_planes, sh_us_planes = ROOT.snd.analysis_tools.GetShoweringPlanes(scifi_planes, us_planes) + ref_point, shower_direction = ROOT.snd.analysis_tools.GetShowerInterceptAndDirection(configuration, sh_scifi_planes, sh_us_planes) + is_nan = np.isnan([ref_point.X(), ref_point.Y(), ref_point.Z(), shower_direction.X(), shower_direction.Y(), shower_direction.Z()]).any() + if is_nan: + print("Found shower direction and/or intercept contain NaN value") + # try finding the shower origin + else: + shower_start = 10*ROOT.snd.analysis_tools.GetScifiShowerStart(scifi_planes) + if shower_start<0: + print("Could not find shower start in SciFi, going for US") + shower_start = 100*ROOT.snd.analysis_tools.GetUSShowerStart(us_planes) + if shower_start<0: + print("Could not find shower start in SciFi or in US") + else: + for k in proj: + if k==1: + drawShowerAxis(h['simpleDisplay'], k, shower_start, ref_point.X(), + shower_direction.X()/shower_direction.Z()) + else: + drawShowerAxis(h['simpleDisplay'], k, shower_start, ref_point.Y(), + shower_direction.Y()/shower_direction.Z()) + h['simpleDisplay'].Update() if option == "2tracks": - rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=0.5) - if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=0.75) - if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=1.0) - if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=1.75) - if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=2.5) - if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=3.0) + rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=0.5, darkMode=darkMode) + if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=0.75, darkMode=darkMode) + if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=1.0, darkMode=darkMode) + if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=1.75, darkMode=darkMode) + if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=2.5, darkMode=darkMode) + if not rc: rc = twoTrackEvent(sMin=10,dClMin=7,minDistance=0.5,sepDistance=3.0, darkMode=darkMode) if verbose>0: dumpChannels() userProcessing(event) - if save: h['simpleDisplay'].Print('{:0>2d}-event_{:04d}'.format(runId,N)+'.png') + if save: + h['simpleDisplay'].Print('{:0>2d}-event_{:04d}'.format(runId, N) + '.' + options.extension) if auto: - h['simpleDisplay'].Print(options.storePic+str(runId)+'-event_'+str(event.EventHeader.GetEventNumber())+'.png') + h['simpleDisplay'].Print(options.storePic + str(runId) + '-event_' + str(event.EventHeader.GetEventNumber()) + '.' + options.extension) if not auto: rc = input("hit return for next event or p for print or q for quit: ") if rc=='p': - h['simpleDisplay'].Print(options.storePic+str(runId)+'-event_'+str(event.EventHeader.GetEventNumber())+'.png') + h['simpleDisplay'].Print(options.storePic + str(runId) + '-event_' + str(event.EventHeader.GetEventNumber()) + '.' + options.extension) elif rc == 'q': break else: eventComment[f"{runId}-event_{event.EventHeader.GetEventNumber()}"] = rc - if save: os.system("convert -delay 60 -loop 0 event*.png animated.gif") + if save: + os.system("convert -delay 60 -loop 0 event*." + options.extension + " animated.gif") -def addTrack(OT,scifi=False): +def addTrack(OT,scifi=False, darkMode=False): xax = h['xz'].GetXaxis() nTrack = 0 for aTrack in OT.Reco_MuonTracks: trackColor = ROOT.kRed if aTrack.GetUniqueID()==1: - trackColor = ROOT.kBlue+2 + trackColor = ROOT.kBlue+2 if not darkMode else ROOT.kBlue-4 flightDir = trackTask.trackDir(aTrack) print('flight direction: %5.3F significance: %5.3F'%(flightDir[0],flightDir[1])) - if aTrack.GetUniqueID()==3: trackColor = ROOT.kBlack - if aTrack.GetUniqueID()==11: trackColor = ROOT.kAzure-2 # HT scifi track - if aTrack.GetUniqueID()==13: trackColor = ROOT.kGray+2 # HT ds track + if aTrack.GetUniqueID()==3: trackColor = ROOT.kBlack if not darkMode else ROOT.kWhite + if aTrack.GetUniqueID()==11: trackColor = ROOT.kAzure-2 if not darkMode else ROOT.kAzure-3 # HT scifi track + if aTrack.GetUniqueID()==13: trackColor = ROOT.kGray+2 if not darkMode else ROOT.kGray # HT ds track # HT cross-system track fit if aTrack.GetUniqueID()==15: trackColor = ROOT.kOrange+7 S = aTrack.getFitStatus() @@ -610,7 +749,8 @@ def addTrack(OT,scifi=False): for p in [0,1]: h['aLine'+str(nTrack*10+p)] = ROOT.TGraph() - zEx = xax.GetBinCenter(1) + # draw the track line starting from the most upstream Veto plane + zEx = h['veto0_z'] mom = aTrack.getFittedState().getMom() pos = aTrack.getFittedState().getPos() lam = (zEx-pos.z())/mom.z() @@ -639,7 +779,7 @@ def addTrack(OT,scifi=False): h[ 'simpleDisplay'].Update() nTrack+=1 -def twoTrackEvent(sMin=10,dClMin=7,minDistance=1.5,sepDistance=0.5): +def twoTrackEvent(sMin=10,dClMin=7,minDistance=1.5,sepDistance=0.5, darkMode=False): trackTask.clusScifi.Clear() trackTask.scifiCluster() clusters = trackTask.clusScifi @@ -723,40 +863,40 @@ def twoTrackEvent(sMin=10,dClMin=7,minDistance=1.5,sepDistance=0.5): if len(tracks)==2: OT = sink.GetOutTree() OT.Reco_MuonTracks = tracks - addTrack(OT,True) + addTrack(OT,True, darkMode=darkMode) return passed -def drawDetectors(): - nodes = {'volMuFilter_1/volFeBlockEnd_1':ROOT.kGreen-6} +def drawDetectors(darkMode=False): + nodes = {'volMuFilter_1/volFeBlockEnd_1':ROOT.kGreen-6 if not darkMode else ROOT.kGreen-2} for i in range(mi.NVetoPlanes): nodes['volVeto_1/volVetoPlane_{}_{}'.format(i, i)]=ROOT.kRed for j in range(mi.NVetoBars): if i<2: nodes['volVeto_1/volVetoPlane_{}_{}/volVetoBar_1{}{:0>3d}'.format(i, i, i, j)]=ROOT.kRed if i==2: nodes['volVeto_1/volVetoPlane_{}_{}/volVetoBar_ver_1{}{:0>3d}'.format(i, i, i, j)]=ROOT.kRed - if i<2: nodes['volVeto_1/subVetoBox_{}'.format(i)]=ROOT.kGray+1 - if i==2: nodes['volVeto_1/subVeto3Box_{}'.format(i)]=ROOT.kGray+1 + if i<2: nodes['volVeto_1/subVetoBox_{}'.format(i)]=ROOT.kGray+1 if not darkMode else ROOT.kGray+2 + if i==2: nodes['volVeto_1/subVeto3Box_{}'.format(i)]=ROOT.kGray+1 if not darkMode else ROOT.kGray+2 for i in range(si.nscifi): # number of scifi stations - nodes['volTarget_1/ScifiVolume{}_{}000000'.format(i+1, i+1)]=ROOT.kBlue+1 + nodes['volTarget_1/ScifiVolume{}_{}000000'.format(i+1, i+1)]=ROOT.kBlue+1 if not darkMode else ROOT.kCyan-6 # iron blocks btw SciFi planes in the testbeam 2023-2024 det layout - nodes['volTarget_1/volFeTarget{}_1'.format(i+1)]=ROOT.kGreen-6 + nodes['volTarget_1/volFeTarget{}_1'.format(i+1)]=ROOT.kGreen-6 if not darkMode else ROOT.kGreen-2 for i in range(em.wall): # number of target walls - nodes['volTarget_1/volWallborder_{}'.format(i)]=ROOT.kGray + nodes['volTarget_1/volWallborder_{}'.format(i)]=ROOT.kGray if not darkMode else ROOT.kGray+2 for i in range(mi.NDownstreamPlanes): - nodes['volMuFilter_1/volMuDownstreamDet_{}_{}'.format(i, i+mi.NVetoPlanes+mi.NUpstreamPlanes)]=ROOT.kBlue+1 + nodes['volMuFilter_1/volMuDownstreamDet_{}_{}'.format(i, i+mi.NVetoPlanes+mi.NUpstreamPlanes)]=ROOT.kBlue+1 if not darkMode else ROOT.kCyan-6 for j in range(mi.NDownstreamBars): - nodes['volMuFilter_1/volMuDownstreamDet_{}_{}/volMuDownstreamBar_ver_3{}{:0>3d}'.format(i, i+mi.NVetoPlanes+mi.NUpstreamPlanes, i, j+mi.NDownstreamBars)]=ROOT.kBlue+1 + nodes['volMuFilter_1/volMuDownstreamDet_{}_{}/volMuDownstreamBar_ver_3{}{:0>3d}'.format(i, i+mi.NVetoPlanes+mi.NUpstreamPlanes, i, j+mi.NDownstreamBars)]=ROOT.kBlue+1 if not darkMode else ROOT.kCyan-6 if i < 3: - nodes['volMuFilter_1/volMuDownstreamDet_{}_{}/volMuDownstreamBar_hor_3{}{:0>3d}'.format(i, i+mi.NVetoPlanes+mi.NUpstreamPlanes, i, j)]=ROOT.kBlue+1 + nodes['volMuFilter_1/volMuDownstreamDet_{}_{}/volMuDownstreamBar_hor_3{}{:0>3d}'.format(i, i+mi.NVetoPlanes+mi.NUpstreamPlanes, i, j)]=ROOT.kBlue+1 if not darkMode else ROOT.kCyan-6 for i in range(mi.NDownstreamPlanes): - nodes['volMuFilter_1/subDSBox_{}'.format(i+mi.NVetoPlanes+mi.NUpstreamPlanes)]=ROOT.kGray+1 + nodes['volMuFilter_1/subDSBox_{}'.format(i+mi.NVetoPlanes+mi.NUpstreamPlanes)]=ROOT.kGray+1 if not darkMode else ROOT.kGray+2 for i in range(mi.NUpstreamPlanes): - nodes['volMuFilter_1/subUSBox_{}'.format(i+mi.NVetoPlanes)]=ROOT.kGray+1 - nodes['volMuFilter_1/volMuUpstreamDet_{}_{}'.format(i, i+mi.NVetoPlanes)]=ROOT.kBlue+1 + nodes['volMuFilter_1/subUSBox_{}'.format(i+mi.NVetoPlanes)]=ROOT.kGray+1 if not darkMode else ROOT.kGray+2 + nodes['volMuFilter_1/volMuUpstreamDet_{}_{}'.format(i, i+mi.NVetoPlanes)]=ROOT.kBlue+1 if not darkMode else ROOT.kCyan-6 for j in range(mi.NUpstreamBars): - nodes['volMuFilter_1/volMuUpstreamDet_{}_{}/volMuUpstreamBar_2{}00{}'.format(i, i+mi.NVetoPlanes, i, j)]=ROOT.kBlue+1 - nodes['volMuFilter_1/volFeBlock_{}'.format(i)]=ROOT.kGreen-6 + nodes['volMuFilter_1/volMuUpstreamDet_{}_{}/volMuUpstreamBar_2{}00{}'.format(i, i+mi.NVetoPlanes, i, j)]=ROOT.kBlue+1 if not darkMode else ROOT.kCyan-6 + nodes['volMuFilter_1/volFeBlock_{}'.format(i)]=ROOT.kGreen-6 if not darkMode else ROOT.kGreen-2 for i in range(mi.NVetoPlanes+mi.NUpstreamPlanes,mi.NVetoPlanes+mi.NUpstreamPlanes+mi.NDownstreamPlanes): - nodes['volMuFilter_1/volFeBlock_{}'.format(i)]=ROOT.kGreen-6 + nodes['volMuFilter_1/volFeBlock_{}'.format(i)]=ROOT.kGreen-6 if not darkMode else ROOT.kGreen-2 passNodes = {'Block', 'Wall', 'FeTarget'} xNodes = {'UpstreamBar', 'VetoBar', 'hor'} proj = {'X':0,'Y':1} @@ -793,6 +933,17 @@ def drawDetectors(): for C in P: M[C] = array('d',[0,0,0]) nav.LocalToMaster(P[C],M[C]) + if "volVetoPlane_0" in node: + h['veto0_z'] = M['LeftBottom'][2] + if "ScifiVolume" in node: + for st in range(1,si.nscifi+1): + if f"{st}_{st}000000" in node: + h[f"scifi{st}_z"] = M['LeftBottom'][2] + if "volMuUpstreamDet" in node: + for st in range(mi.NUpstreamPlanes): + if f"{st}_{st+mi.NVetoPlanes}" in node: + h[f"us{st+1}_z"] = M['LeftBottom'][2] + h[node+p] = ROOT.TPolyLine() X = h[node+p] c = proj[p] @@ -1056,12 +1207,12 @@ def dumpChannels(D='Digi_MuFilterHits'): keys.sort() for k in keys: print(text[k]) -def fillNode(node, color=None): +def fillNode(node, color=None, darkMode=False): xNodes = {'UpstreamBar', 'VetoBar', 'hor'} proj = {'X':0,'Y':1} if color == None : - hcal_color = ROOT.kBlack - veto_color = ROOT.kRed+1 + hcal_color = ROOT.kBlack if not darkMode else ROOT.kWhite + veto_color = ROOT.kRed+1 if not darkMode else ROOT.kRed-4 else : hcal_color = color veto_color = color @@ -1084,17 +1235,17 @@ def fillNode(node, color=None): X.Draw('f&&same') X.Draw('same') -def drawInfo(pad, k, run, event, timestamp,moreEventInfo=[]): +def drawInfo(pad, k, run, event, timestamp,moreEventInfo=[], darkMode=False): drawLogo = True drawText = True if drawLogo: padLogo = ROOT.TPad("logo","logo",0.1,0.1,0.2,0.3) padLogo.SetFillStyle(4000) padLogo.SetFillColorAlpha(0, 0) + if darkMode: padLogo.SetFillColor(ROOT.kBlack) padLogo.Draw() - logo = ROOT.TImage.Open('$SNDSW_ROOT/shipLHC/Large__SND_Logo_black_cut.png') + logo = ROOT.TImage.Open('$SNDSW_ROOT/shipLHC/Large__SND_Logo_black_cut.png') if not darkMode else ROOT.TImage.Open('$SNDSW_ROOT/shipLHC/Large__SND_Logo_white_cut.png') logo.SetConstRatio(True) - logo.DrawText(0, 0, 'SND', 98) padLogo.cd() logo.Draw() pad.cd(k) @@ -1114,6 +1265,7 @@ def drawInfo(pad, k, run, event, timestamp,moreEventInfo=[]): textInfo.SetTextAlign(11) textInfo.SetTextFont(42) textInfo.SetTextSize(.15) + if darkMode: textInfo.SetTextColor(ROOT.kWhite) textInfo.DrawLatex(0, 0.6, 'SND@LHC Experiment, CERN') if hasattr(eventTree.EventHeader,'GetEventNumber'): N = eventTree.EventHeader.GetEventNumber() else: N = event @@ -1130,10 +1282,42 @@ def drawInfo(pad, k, run, event, timestamp,moreEventInfo=[]): textInfo.SetTextAlign(11) textInfo.SetTextFont(42) textInfo.SetTextSize(.1) - textInfo.SetTextColor(ROOT.kMagenta+2) + textInfo.SetTextColor((ROOT.kMagenta+2) if not darkMode else (ROOT.kMagenta-4)) dely = 0.12 ystart = 0.85 for i in range(7): textInfo.DrawLatex(0.4, 0.9-dely*i, moreEventInfo[i]) pad.cd(k) +def drawCollisionAxis(pad, k): + line_name = "collision_axis_line_" + str(k) + h[line_name] = ROOT.TLine(h["zmin"], 0, h["zmax"], 0) + h[line_name].SetLineColor(ROOT.kRed) + h[line_name].SetLineStyle(2) + + text_name = "collision_axis_text_" + str(k) + h[text_name] = ROOT.TText(h["zmin"] + 8, 0 + 2, "Collision axis") + h[text_name].SetTextAlign(12) + h[text_name].SetTextFont(43) + h[text_name].SetTextSize(13 * resolution_factor) + h[text_name].SetTextColor(ROOT.kRed) + + pad.cd(k) + h[line_name].Draw() + h[text_name].Draw() + +def drawShowerAxis(pad, k, shower_start, ref_point, shower_direction): + line_name= "shower_dir_" + str(k) + if shower_start<100: + zpos_name='scifi'+str(shower_start//10) + else: + zpos_name='us'+str(shower_start//100) + h[line_name] = ROOT.TLine(h[zpos_name+'_z'], + ref_point+shower_direction*h[zpos_name+'_z'], + h['zmax'], + ref_point+shower_direction*h['zmax']) + h[line_name].SetLineColor(ROOT.kPink-6) + h[line_name].SetLineStyle(3) + h[line_name].SetLineWidth(5) + pad.cd(k) + h[line_name].Draw() diff --git a/shipLHC/scripts/2dMuEventBuilderDisplay.py b/shipLHC/scripts/2dMuEventBuilderDisplay.py index eac251b884..47372acfc3 100644 --- a/shipLHC/scripts/2dMuEventBuilderDisplay.py +++ b/shipLHC/scripts/2dMuEventBuilderDisplay.py @@ -33,13 +33,13 @@ mc = False if options.inputFile=="": f=ROOT.TFile('sndsw_raw_'+str(options.runNumber).zfill(6)+'.root') - eventTree = f.rawConv + eventTree = f.Get("rawConv") else: f=ROOT.TFile.Open(options.inputFile) if f.FindKey('cbmsim'): - eventTree = f.cbmsim + eventTree = f.Get("cbmsim") mc = True - else: eventTree = f.rawConv + else: eventTree = f.Get("rawConv") nav = ROOT.gGeoManager.GetCurrentNavigator() for p in [0,1]: diff --git a/shipLHC/scripts/FillingScheme.py b/shipLHC/scripts/FillingScheme.py index b43b0e3281..f54154ad96 100644 --- a/shipLHC/scripts/FillingScheme.py +++ b/shipLHC/scripts/FillingScheme.py @@ -117,6 +117,7 @@ def getNameOfFillingscheme(self,fillnr): Y = "2022" if fillnr > 8500 : Y = "2023" if fillnr >= 9323 : Y="2024" + if fillnr >= 10406 : Y="2025" if not self.lpcFillingscheme: with urlopen('https://lpc.web.cern.ch/cgi-bin/fillTable.py?year='+Y) as webpage: self.lpcFillingscheme = webpage.read().decode() @@ -159,10 +160,10 @@ def getFillNrFromRunNr(self,runNumber): options.rawData+"/run_"+str(runNumber).zfill(6)+"/data_0000.root") try: if R.Get('event'): - rc = R.event.GetEvent(R.event.GetEntries()-1) - flags = R.event.flags + rc = R.Get("event").GetEvent(R.Get("event").GetEntries()-1) + flags = R.Get("event").flags else: - event = R.data + event = R.Get("data") event.GetEvent(0) flags = event.evt_flags fillNumber = numpy.bitwise_and(FILL_NUMBER_MASK,flags) @@ -177,6 +178,7 @@ def getLumiAtIP1(self,fillnr=None,fromnxcals=False, fromAtlas=False): Y = "2022" if fillnr>8500: Y = "2023" if fillnr>=9323: Y="2024" + if fillnr>=10406 : Y="2025" if not fromnxcals and not fromAtlas: try: with urlopen('https://lpc.web.cern.ch/cgi-bin/fillAnalysis.py?year='+Y+'&action=fillData&exp=ATLAS&fillnr='+str(fillnr)) as webpage: @@ -259,7 +261,7 @@ def getLumiAtIP1(self,fillnr=None,fromnxcals=False, fromAtlas=False): def drawLumi(self,runNumber): R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root") ROOT.gROOT.cd() - bCanvas = R.daq.Get('T') + bCanvas = R.Get("daq").Get('T') Xt = {'time':None,'timeWtDS':None,'timeWt':None} for x in Xt: Xt[x] = bCanvas.FindObject(x) @@ -440,12 +442,12 @@ def extractFillingScheme(self,fillNr): print('Name of filling scheme: ',fs_name_table) if fs_name_table==0: return -1 - # since 2024 new there is new storage of FS in json files, no csv + # since 2024 there is new storage of FS in json files, no csv fs_url="https://gitlab.cern.ch/lhc-injection-scheme/injection-schemes/-/raw/master/"+\ fs_name_table+".json" F = ROOT.TFile(self.path+'fillingScheme-'+fillNr+'.root','recreate') nt = ROOT.TNtuple('fill'+fillNr,'b1 IP1 IP2','B1:IP1:IP2:IsB2') - # Get the 2024 data + # Get the data for year>=2024 # 9323 is first fillNr for 2024 if int(fillNr) >= 9323: # Get the JSON file content from the web @@ -582,9 +584,9 @@ def extractPhaseShift(self,fillNr,runNumber): R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root") ROOT.gROOT.cd() try: - self.h['bnr'] = R.daq.Get('bunchNumber').FindObject('bnr').Clone('bnr') + self.h['bnr'] = R.Get("daq").Get('bunchNumber').FindObject('bnr').Clone('bnr') except: - self.h['bnr'] = R.daq.Get('shifter/bunchNumber').FindObject('bnr').Clone('bnr') + self.h['bnr'] = R.Get("daq").Get('shifter/bunchNumber').FindObject('bnr').Clone('bnr') R.Close() # create the bunch number plot if offline monitoring file is missing except: @@ -684,9 +686,9 @@ def plotBunchStructure(self,fillNr,runNumber): try: R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root") ROOT.gROOT.cd() - bCanvas = R.daq.Get('bunchNumber') + bCanvas = R.Get("daq").Get('bunchNumber') if not bCanvas: - bCanvas = R.daq.shifter.Get('bunchNumber') + bCanvas = R.Get("daq").Get("shifter").Get('bunchNumber') h['bnr']= bCanvas.FindObject('bnr').Clone('bnr') # create the bunch number plot if offline monitoring file is missing except: @@ -765,7 +767,7 @@ def Xbunch(self): h['XIP2z'].SetStats(0) R = ROOT.TFile.Open(www+"offline/run"+str(runNumber).zfill(6)+".root") ROOT.gROOT.cd() - bCanvas = R.daq.Get('sndclock') + bCanvas = R.Get("daq").Get('sndclock') h['Xbnr']= bCanvas.FindObject('Xbnr').Clone('Xbnr') ROOT.gROOT.cd() h['c1'].cd() @@ -980,7 +982,7 @@ def FwBw(self,runNumber): xing = {'all':False,'B1only':False,'B2noB1':False,'noBeam':False} for T in ['Txing','TD','T']: - h[T] = self.F.daq.Get(T).Clone(T) + h[T] = self.F.Get("daq").Get(T).Clone(T) for X in ['time','timeWt','timeWtDS']: h[X] = h['T'].FindObject(X).Clone(X) for t in ['time','timeWt','timeWtDS','bnr']: @@ -1158,9 +1160,9 @@ def FwBw(self,runNumber): for B in ['B2noB1','B1only','noBeam']: for x in ['scifi-trackDir','scifi-TtrackPos']: T = x+B - tmp = self.F.scifi.Get(T) + tmp = self.F.Get("scifi").Get(T) if not tmp: - tmp = self.F.scifi.Get(B+'/'+T) + tmp = self.F.Get("scifi").Get(B+'/'+T) h[T] = tmp.Clone(T) if x.find('Dir')>0: for y in ['scifi-trackSlopesXL','slopeXL','slopeYL']: @@ -1189,9 +1191,9 @@ def FwBw(self,runNumber): for x in ['muonDSTracks','mufi-TtrackPos']: T = x+B - tmp = self.F.mufilter.Get(T) + tmp = self.F.Get("mufilter").Get(T) if not tmp: - tmp = self.F.mufilter.Get(B+'/'+T) + tmp = self.F.Get("mufilter").Get(B+'/'+T) h[T] = tmp.Clone(T) if x.find('DSTrack')>0: for y in ['mufi-slopes','slopeX','slopeY']: @@ -1643,7 +1645,7 @@ def runsWithBeam(self): runNumber = int(x[ir+4:ir+4+k+1]) if runNumber0: + # extract the em target run from the rawData path + em_run = options.rawData[options.rawData.find("run_"):] + options.convpath = "/eos/experiment/sndlhc/convertedData/physics/2025/"+em_run + options.rmin = 10919-1 #10587 fillN + offline =www+"offline.html" FS = fillingScheme() FS.Init(options) ut.bookCanvas(FS.h,'c1','c1',1800,900,1,1) diff --git a/shipLHC/scripts/Monitor.py b/shipLHC/scripts/Monitor.py index 9bf7b09c79..bdab6763ba 100644 --- a/shipLHC/scripts/Monitor.py +++ b/shipLHC/scripts/Monitor.py @@ -160,7 +160,7 @@ def __init__(self,options,FairTasks): f=ROOT.TFile.Open(options.fname) eventChain = f.Get('rawConv') if not eventChain: - eventChain = f.cbmsim + eventChain = f.Get("cbmsim") if eventChain.GetBranch('MCTrack'): self.MonteCarlo = True partitions = [] else: @@ -268,7 +268,7 @@ def __init__(self,options,FairTasks): def GetEntries(self): if self.options.online: if self.converter.newFormat: return self.converter.fiN.Get('data').GetEntries() - else: return self.converter.fiN.event.GetEntries() + else: return self.converter.fiN.Get("event").GetEntries() else: return self.eventTree.GetEntries() @@ -848,7 +848,8 @@ def ExecuteEvent(self,event): for item in track_container_list: for aTrack in item: i_muon += 1 - self.fittedTracks[i_muon] = aTrack + aTrack_TCA = self.fittedTracks.ConstructedAt(i_muon) + ROOT.std.swap(aTrack, aTrack_TCA) def Execute(self): for n in range(self.options.nStart,self.options.nStart+self.options.nEvents): @@ -867,13 +868,15 @@ def Execute(self): self.clusScifi.Delete() self.clusScifi.Expand(len(self.trackTask.clusScifi)) for index, aCl in enumerate(self.trackTask.clusScifi): - self.clusScifi[index] = aCl + aCl_TCA = self.clusScifi.ConstructedAt(index) + ROOT.std.swap(aCl, aCl_TCA) if self.options.simpleTracking and not self.options.trackType.find('DS')<0: if not self.eventTree.GetBranch("Cluster_Mufi"): self.clusMufi.Delete() self.clusMufi.Expand(len(self.trackTask.clusMufi)) for index, aCl in enumerate(self.trackTask.clusMufi): - self.clusMufi[index] = aCl + aCl_TCA = self.clusMufi.ConstructedAt(index) + ROOT.std.swap(aCl, aCl_TCA) # if using FairEventHeader, i.e. before sndlhc header was introduced if hasattr(self.OT.EventHeader, "SetMCEntryNumber"): diff --git a/shipLHC/scripts/MufiCTR.py b/shipLHC/scripts/MufiCTR.py index 572d82a6fc..57df6e7fd5 100644 --- a/shipLHC/scripts/MufiCTR.py +++ b/shipLHC/scripts/MufiCTR.py @@ -22,11 +22,11 @@ def start(): if d!='': dd = d+'/' # safety net if B2noB1 and B1Only plots don't exist try: - h['mufi-'+x+d] = f.mufilter.Get(dd+'mufi-'+x+d).Clone('mufi-'+x+d) + h['mufi-'+x+d] = f.Get("mufilter").Get(dd+'mufi-'+x+d).Clone('mufi-'+x+d) except: try: #2024: changed structure of offline monitoring file with shifter and expert subdirs - h['mufi-'+x+d] = f.mufilter.Get('expert/'+dd+'mufi-'+x+d).Clone('mufi-'+x+d) + h['mufi-'+x+d] = f.Get("mufilter").Get('expert/'+dd+'mufi-'+x+d).Clone('mufi-'+x+d) except: continue h['mufi-'+x+d].Draw() diff --git a/shipLHC/scripts/Survey-MufiScifi.py b/shipLHC/scripts/Survey-MufiScifi.py index 80c79e5d23..75f12dfaeb 100644 --- a/shipLHC/scripts/Survey-MufiScifi.py +++ b/shipLHC/scripts/Survey-MufiScifi.py @@ -326,8 +326,8 @@ def makeAnimation(histname,j0=1,j1=2,animated=True, findMinMax=True, lim = 50): else: # for MC data and other files f=ROOT.TFile.Open(options.fname) - if f.Get('rawConv'): eventChain = f.rawConv - else: eventChain = f.cbmsim + if f.Get('rawConv'): eventChain = f.Get("rawConv") + else: eventChain = f.Get("cbmsim") if options.remakeScifiClusters: eventChain.SetBranchStatus("Cluster_Scifi*",0) rc = eventChain.GetEvent(0) run = ROOT.FairRunAna() @@ -2832,7 +2832,7 @@ def analyze_EfficiencyAndResiduals(readHists=False,mode='S',local=True,zoom=Fals hist.GetYaxis().SetRangeUser(ymin,ymax) hist.Draw('colz') # get time x correlation, X = m*dt + b - h['gdtLRvsX_'+key] = ROOT.TGraph() + h['gdtLRvsX_'+key] = ROOT.TGraphErrors() g = h['gdtLRvsX_'+key] xproj = hist.ProjectionX('tmpx') if xproj.GetSumOfWeights()==0: continue @@ -2844,9 +2844,14 @@ def analyze_EfficiencyAndResiduals(readHists=False,mode='S',local=True,zoom=Fals X = hist.GetXaxis().GetBinCenter(nx) rc = tmp.Fit('gaus','NQS') res = rc.Get() - if not res: dt = tmp.GetMean() - else: dt = res.Parameter(1) + if not res: + dt = tmp.GetMean() + err = tmp.GetStd() + else: + dt = res.Parameter(1) + err =res.Parameter(2) g.SetPoint(np,X,dt) + g.SetPointError(np, 0, err) np+=1 g.SetLineColor(ROOT.kRed) g.SetLineWidth(2) diff --git a/shipLHC/thermalNeutrons.py b/shipLHC/thermalNeutrons.py index 373875a640..5cb57c82bd 100644 --- a/shipLHC/thermalNeutrons.py +++ b/shipLHC/thermalNeutrons.py @@ -140,7 +140,7 @@ def count(hFile = outFile): ut.bookHist(h,'dxyz','',100,-0.1,0.1,100,-0.1,0.1,200,-1.,1.) ut.bookHist(h,'Epassed',';log10(Ekin [MeV])',1000,-10.,4.) Npassed = 0 - for sTree in f.cbmsim: + for sTree in f.Get("cbmsim"): N = sTree.MCTrack[0] Ekin = N.GetP()**2/(2*N.GetMass())*1000. logEkin = ROOT.TMath.Log10(Ekin) diff --git a/shipgen/GenieGenerator.cxx b/shipgen/GenieGenerator.cxx index 9420abe41f..e1cdd0ed63 100644 --- a/shipgen/GenieGenerator.cxx +++ b/shipgen/GenieGenerator.cxx @@ -15,6 +15,7 @@ #include "TGeoCompositeShape.h" #include "TParticle.h" #include "TClonesArray.h" +#include "FairMCEventHeader.h" using std::cout; using std::endl; @@ -52,6 +53,8 @@ Bool_t GenieGenerator::Init(const char* fileName, const int firstEvent) { fTree = (TTree *)fInputFile->Get("gst"); fNevents = fTree->GetEntries(); fn = firstEvent; + fTree->SetBranchAddress("FlukaRun",&runN); // run number + fTree->SetBranchAddress("evtNumber",&eventN); // event number fTree->SetBranchAddress("Ev",&Ev); // incoming neutrino energy fTree->SetBranchAddress("pxv",&pxv); fTree->SetBranchAddress("pyv",&pyv); @@ -494,7 +497,12 @@ Bool_t GenieGenerator::ReadEvent(FairPrimaryGenerator* cpg) 0, true); } - + + // Set the generation run and event numbers through the MC event header + FairMCEventHeader* header = cpg->GetEvent(); + header->SetEventID(eventN); + header->SetRunID(runN); + return kTRUE; } else { @@ -706,6 +714,12 @@ Bool_t GenieGenerator::ReadEvent(FairPrimaryGenerator* cpg) //cout << "Info GenieGenerator Return from GenieGenerator" << endl; } } + + // Set the generation run and event numbers through the MC event header + FairMCEventHeader* header = cpg->GetEvent(); + header->SetEventID(eventN); + header->SetRunID(runN); + return kTRUE; } // ------------------------------------------------------------------------- diff --git a/shipgen/GenieGenerator.h b/shipgen/GenieGenerator.h index 35bfb23e0e..b087944791 100644 --- a/shipgen/GenieGenerator.h +++ b/shipgen/GenieGenerator.h @@ -1,5 +1,5 @@ #ifndef PNDGeGENERATOR_H -#define PNDGeGENERATOR_H 1 +#define PNDGeGENERATOR_H 3 #include "TROOT.h" #include "FairGenerator.h" @@ -52,6 +52,7 @@ class GenieGenerator : public FairGenerator private: protected: + int runN,eventN; Double_t FLUKA_x_cos, FLUKA_y_cos, FLUKA_x, FLUKA_y; Double_t fDeltaE_GenieFLUKA_nu; Double_t Yvessel,Xvessel,Lvessel,ztarget,startZ,endZ; @@ -78,7 +79,7 @@ class GenieGenerator : public FairGenerator TH1D* pyslice[3000][500];//! TClonesArray *ancstr;//! - ClassDef(GenieGenerator,2); + ClassDef(GenieGenerator,3); }; #endif /* !PNDGeGENERATOR_H */ diff --git a/shipgen/MuDISGenerator.cxx b/shipgen/MuDISGenerator.cxx index 211bfae8f9..68b790a7aa 100644 --- a/shipgen/MuDISGenerator.cxx +++ b/shipgen/MuDISGenerator.cxx @@ -13,6 +13,7 @@ #include "TVectorD.h" #include "TParticle.h" #include "TGeoCompositeShape.h" +#include "FairMCEventHeader.h" using std::cout; using std::endl; @@ -49,6 +50,8 @@ Bool_t MuDISGenerator::Init(const char* fileName, const int firstEvent) { fTree = (TTree *)fInputFile->Get("DIS"); fNevents = fTree->GetEntries(); fn = firstEvent; + fTree->SetBranchAddress("run",&runN); // run number + fTree->SetBranchAddress("event",&eventN); // event number fTree->SetBranchAddress("InMuon",&iMuon); // incoming muon fTree->SetBranchAddress("Particles",&dPart); return kTRUE; @@ -286,6 +289,12 @@ Bool_t MuDISGenerator::ReadEvent(FairPrimaryGenerator* cpg) TParticle* Part = dynamic_cast(dPart->AddrAt(i)); cpg->AddTrack(Part->GetPdgCode(),Part->Px(),Part->Py(),Part->Pz(),xmu,ymu,zmu,0,true,Part->Energy(),mu->T(),w); } + + // Set the generation run and event numbers through the MC event header + FairMCEventHeader* header = cpg->GetEvent(); + header->SetEventID(eventN); + header->SetRunID(runN); + return kTRUE; } // ------------------------------------------------------------------------- diff --git a/shipgen/MuDISGenerator.h b/shipgen/MuDISGenerator.h index 6d6dd4606e..5d2916f586 100644 --- a/shipgen/MuDISGenerator.h +++ b/shipgen/MuDISGenerator.h @@ -1,5 +1,5 @@ #ifndef PNDMuGENERATOR_H -#define PNDMuGENERATOR_H 1 +#define PNDMuGENERATOR_H 2 #include "TROOT.h" #include "FairGenerator.h" @@ -40,13 +40,14 @@ class MuDISGenerator : public FairGenerator protected: Double_t startZ,endZ; TClonesArray* iMuon ; - TClonesArray* dPart ; + TClonesArray* dPart ; + int runN, eventN; FairLogger* fLogger; //! don't make it persistent, magic ROOT command TFile* fInputFile; TTree* fTree; int fNevents; int fn; bool fFirst; - ClassDef(MuDISGenerator,1); + ClassDef(MuDISGenerator,2); }; #endif /* !PNDMuGENERATOR_H */ diff --git a/shipgen/NtupleGenerator_FLUKA.cxx b/shipgen/NtupleGenerator_FLUKA.cxx index e01cd2bc8a..404f6ba531 100644 --- a/shipgen/NtupleGenerator_FLUKA.cxx +++ b/shipgen/NtupleGenerator_FLUKA.cxx @@ -5,6 +5,7 @@ #include "FairPrimaryGenerator.h" #include "NtupleGenerator_FLUKA.h" #include "FairLogger.h" +#include "FairMCEventHeader.h" // read events from ntuples produced by FLUKA @@ -39,6 +40,8 @@ Bool_t NtupleGenerator_FLUKA::Init(const char* fileName, const int firstEvent) { fTree->SetBranchAddress("primaries",&primaries); } else { + fTree->SetBranchAddress("run",&runN); // run number + fTree->SetBranchAddress("event",&eventN); // event number fTree->SetBranchAddress("id",&id); // particle id fTree->SetBranchAddress("generation",&generation); // origin generation number fTree->SetBranchAddress("t",&t); // time of flight @@ -83,7 +86,7 @@ Bool_t NtupleGenerator_FLUKA::ReadEvent(FairPrimaryGenerator* cpg) for(int i=0; i(primaries->AddrAt(i)); - // what to do with generation info? convert time from ns to sec + //convert time from ns to sec cpg->AddTrack(int(primary->id),primary->px,primary->py,primary->pz, primary->x,primary->y,primary->z-SND_Z,-1,true,primary->E, primary->t/1E9,primary->w,(TMCProcess)primary->generation); @@ -94,11 +97,16 @@ Bool_t NtupleGenerator_FLUKA::ReadEvent(FairPrimaryGenerator* cpg) } } else { - // what to do with generation info? convert time from ns to sec + // convert time from ns to sec cpg->AddTrack(int(id[0]),px[0],py[0],pz[0],x[0],y[0],z[0]-SND_Z,-1,true,E[0],t[0]/1E9,w[0],(TMCProcess)generation[0]); LOG(DEBUG) << "NtupleGenerator_FLUKA: add muon " << id[0]<<","<GetEvent(); + header->SetEventID(int(eventN)); + header->SetRunID(int(runN)); + return kTRUE; } diff --git a/shipgen/NtupleGenerator_FLUKA.h b/shipgen/NtupleGenerator_FLUKA.h index f7756d287b..fd1852bce1 100644 --- a/shipgen/NtupleGenerator_FLUKA.h +++ b/shipgen/NtupleGenerator_FLUKA.h @@ -1,5 +1,5 @@ #ifndef FLUKAntGENERATOR_H -#define FLUKAntGENERATOR_H 1 +#define FLUKAntGENERATOR_H 4 #include "TROOT.h" #include "FairGenerator.h" @@ -32,13 +32,13 @@ class NtupleGenerator_FLUKA : public FairGenerator private: protected: - Double_t id[1],generation[1],E[1],t[1],px[1],py[1],pz[1],x[1],y[1],z[1],w[1],SND_Z; + Double_t runN,eventN,id[1],generation[1],E[1],t[1],px[1],py[1],pz[1],x[1],y[1],z[1],w[1],SND_Z; TFile* fInputFile; TTree* fTree; int fNevents; int fn; TClonesArray* primaries; - ClassDef(NtupleGenerator_FLUKA,3); + ClassDef(NtupleGenerator_FLUKA,4); }; #endif /* !FLUKAntGENERATOR_H */ diff --git a/shipgen/sndlhcGSimpleNtpConverter.cxx b/shipgen/sndlhcGSimpleNtpConverter.cxx index 5faff044c0..10adbf2868 100644 --- a/shipgen/sndlhcGSimpleNtpConverter.cxx +++ b/shipgen/sndlhcGSimpleNtpConverter.cxx @@ -15,199 +15,262 @@ Converts neutrino rays generated by FLUKA into GSimpleNtpFlux format for GENIE */ -int main(int argc, char** argv){ - - if (argc != 4) { - std::cout << "Three arguments required: path to FLUKA file AND output file name AND number of pp collisions used to generate FLUKA file." << std::endl; - return -1; - } - - std::string inFileName = std::string(argv[1]); - std::string outFileName = std::string(argv[2]); - double pp_collision_number = std::stod(argv[3]); - - // Convert FLUKA to PDG particle IDs. (ONLY NEUTRINOS IMPLEMENTED!!!) - std::map FLUKAtoPDG { {27, 14}, - {28, -14}, - {5, 12}, - {6, -12}, - {43, 16}, - {44, -16} }; - - // To generate a random seed which is stored in the metadata. - TRandom3 * ran = new TRandom3(); - - bool verbose = false; - - // Scoring plane info - // Plane z in FLUKA coordinates: 48386 cm - // FLUKA z = 0 in sndsw coordinates : -48000 cm - double z = (48386-48000); // In cm for consistency with the FLUKA file. - - double plane_corner[] = {-70., 5., z}; - double plane_dir1[] = {140., 0, 0}; - double plane_dir2[] = {0, 65., 0}; - - // Set up input file - std::ifstream in_file(inFileName); - - // Input variables - /* - # Scoring particles entering Region No 2935 - # Col 1: FLUKA run number - # Col 2: primary event number - # -- Particle information -- - # Col 3: FLUKA particle type ID - # Col 4: Generation number - # Col 5: Kinetic energy (GeV) - # Col 6: Statistical weight - # -- Crossing at scoring plane -- - # Col 7: x coord (cm) - # Col 8: y coord (cm) - # Col 9: x dir cosine - # Col 10: y dir cosine - # Col 11: Particle age since primary event (sec) - # Col 12: ... +int main(int argc, char **argv) +{ + + if (argc != 4) { + std::cout << "Three arguments required: path to FLUKA file AND output file name AND number of pp collisions used " + "to generate FLUKA file." + << std::endl; + return -1; + } + + std::string inFileName = std::string(argv[1]); + std::string outFileName = std::string(argv[2]); + double pp_collision_number = std::stod(argv[3]); + + // Convert FLUKA to PDG particle IDs. (ONLY NEUTRINOS IMPLEMENTED!!!) + std::map FLUKAtoPDG{{27, 14}, {28, -14}, {5, 12}, {6, -12}, {43, 16}, {44, -16}}; + + // To generate a random seed which is stored in the metadata. + TRandom3 *ran = new TRandom3(); + + bool verbose = false; + + std::cout << "Start converting" << std::endl; + + // Set up input file + std::ifstream in_file(inFileName); + + // Input variables + /* + # Scoring particles entering Region No 2935 + # Col 1: FLUKA run number + # Col 2: primary event number + # -- Particle information -- + # Col 3: FLUKA particle type ID + # Col 4: Generation number + # Col 5: Kinetic energy (GeV) + # Col 6: Statistical weight + # -- Crossing at scoring plane -- + # Col 7: x coord (cm) + # Col 8: y coord (cm) + # Col 9: x dir cosine + # Col 10: y dir cosine + # Col 11: z coord (cm) + # Col 12: Particle age since primary event (sec) + # Col 13: Last decay x cooord (cm) + # Col 14: Last decay y cooord (cm) + # Col 15: Last decay z cooord (cm) + # Col 16: Last decay ID + # Col 17: Kinetic energy of decay parent (GeV) + # Col 18: Last interaction x cooord (cm) + # Col 19: Last interaction y cooord (cm) + # Col 20: Last interaction z cooord (cm) + # Col 21: Last interaction ID + # Col 22: Kinetic energy of parent (GeV) */ - int FlukaRun, evtNumber, FlukaID, generationN; - float Ekin, weight, x, y, x_cos, y_cos, age; - - TFile * fOut = new TFile(outFileName.c_str(), "RECREATE"); - TTree * tOut = new TTree("flux", "a simple flux n-tuple"); - TTree * metaOut = new TTree("meta","metadata for flux n-tuple"); - - genie::flux::GSimpleNtpEntry * gsimple_entry = new genie::flux::GSimpleNtpEntry; - tOut->Branch("entry", &gsimple_entry); - genie::flux::GSimpleNtpMeta* meta_entry = new genie::flux::GSimpleNtpMeta; - metaOut->Branch("meta", &meta_entry); - - genie::flux::GSimpleNtpAux* aux_entry = new genie::flux::GSimpleNtpAux; - tOut->Branch("aux", &aux_entry); - - - UInt_t metakey = TString::Hash(outFileName.c_str(),strlen(outFileName.c_str())); - metakey &= 0x7FFFFFFF; - - // Metadata accumulators: - std::set pdglist; - double min_weight = 1e10; - double max_weight = -1e10; - double max_energy = 0; - - unsigned long int counter = 0; - - if (in_file.is_open()){ - string in_line; - while (!in_file.eof()){ - getline(in_file, in_line); - - // Skip lines containing "#" - if (in_line.find("#") == in_line.npos){ - in_file >> FlukaRun - >> evtNumber - >> FlukaID - >> generationN - >> Ekin - >> weight - >> x - >> y - >> x_cos - >> y_cos - >> age; - - gsimple_entry->Reset(); - aux_entry->Reset(); - - - if (verbose) std::cout << "Got entry " << counter++ << std::endl; - - if (verbose){ - std::cout << FlukaRun << "\n" - << evtNumber << "\n" - << FlukaID << "\n" - << generationN << "\n" - << Ekin << "\n" - << weight << "\n" - << x << "\n" - << y << "\n" - << z << "\n" - << x_cos << "\n" - << y_cos << "\n" - << age << "\n" - << "--------------------------------------------------" << std::endl; - } - - gsimple_entry->metakey = metakey; - gsimple_entry->pdg = FLUKAtoPDG[FlukaID]; - gsimple_entry->wgt = weight; - gsimple_entry->vtxx = x*ShipUnit::cm/ShipUnit::m; // in m - gsimple_entry->vtxy = y*ShipUnit::cm/ShipUnit::m; - gsimple_entry->vtxz = z*ShipUnit::cm/ShipUnit::m; - gsimple_entry->dist = 0.; // Distance from hadron decay point to neutrino "vertex", to use for oscillations, for example. Don't use. - - // I'm assuming x_cos, y_cos are normalized. - double z_cos = sqrt(1 - pow(x_cos, 2) - pow(y_cos, 2)); - // And massless neutrinos - gsimple_entry->px = Ekin*x_cos; // in GeV/c - gsimple_entry->py = Ekin*y_cos; - gsimple_entry->pz = Ekin*z_cos; - gsimple_entry->E = Ekin; - - // Set auxiliary data - aux_entry->auxint.push_back(FlukaRun); - aux_entry->auxint.push_back(evtNumber); - aux_entry->auxint.push_back(FlukaID); - aux_entry->auxint.push_back(generationN); - - aux_entry->auxdbl.push_back(age); - aux_entry->auxdbl.push_back(weight); - - // Accumulate metadata - pdglist.insert(gsimple_entry->pdg); - min_weight = TMath::Min(min_weight, gsimple_entry->wgt); - max_weight = TMath::Max(max_weight, gsimple_entry->wgt); - max_energy = TMath::Max(max_energy, gsimple_entry->E); - - // All done! - tOut->Fill(); + int FlukaRun, evtNumber, FlukaID, generationN, last_decay_ID, last_interaction_ID; + double Ekin, weight, x, y, z, x_cos, y_cos, age, Ekin_parent; + double last_decay_x, last_decay_y, last_decay_z, Ekin_decay; + double last_interaction_x, last_interaction_y, last_interaction_z; + + // Scoring plane info + // Find a Minimum and Maximum of the x, y, z axis (FLUKA coordinates) cm + // Max z, min z + double min_z = 999; + double max_z = -999; + // Max x, min x + double min_x = 999; + double max_x = -999; + // Max y, min y + double min_y = 999; + double max_y = -999; + + TFile *fOut = new TFile(outFileName.c_str(), "RECREATE"); + TTree *tOut = new TTree("flux", "a simple flux n-tuple"); + TTree *metaOut = new TTree("meta", "metadata for flux n-tuple"); + + genie::flux::GSimpleNtpEntry *gsimple_entry = new genie::flux::GSimpleNtpEntry; + tOut->Branch("entry", &gsimple_entry); + genie::flux::GSimpleNtpMeta *meta_entry = new genie::flux::GSimpleNtpMeta; + metaOut->Branch("meta", &meta_entry); + + genie::flux::GSimpleNtpAux *aux_entry = new genie::flux::GSimpleNtpAux; + tOut->Branch("aux", &aux_entry); + + UInt_t metakey = TString::Hash(outFileName.c_str(), strlen(outFileName.c_str())); + metakey &= 0x7FFFFFFF; + + // Metadata accumulators: + std::set pdglist; + double min_weight = 1e10; + double max_weight = -1e10; + double max_energy = 0; + + unsigned long int counter = 0; + + if (in_file.is_open()) { + string in_line; + + while (!in_file.eof()) { + getline(in_file, in_line); + + // Skip lines containing "#" + if (in_line.find("#") == in_line.npos) { + in_file >> FlukaRun >> evtNumber >> FlukaID >> generationN >> Ekin >> weight >> x >> y >> x_cos >> y_cos >> + z >> age >> last_decay_x >> last_decay_y >> last_decay_z >> last_decay_ID >> Ekin_decay >> + last_interaction_x >> last_interaction_y >> last_interaction_z >> last_interaction_ID >> Ekin_parent; + + gsimple_entry->Reset(); + aux_entry->Reset(); + + // FLUKA z = 0 in sndsw coordinates : -48000 cm + z -= 48000; // In cm for consistency with the FLUKA file. + + min_z = std::min(min_z, z); + max_z = std::max(max_z, z); + min_x = std::min(min_x, x); + max_x = std::max(max_x, x); + min_y = std::min(min_y, y); + max_y = std::max(max_y, y); + + if (verbose) + std::cout << "Got entry " << counter++ << std::endl; + + if (verbose) { + std::cout << FlukaRun << "\n" + << evtNumber << "\n" + << FlukaID << "\n" + << generationN << "\n" + << Ekin << "\n" + << weight << "\n" + << x << "\n" + << y << "\n" + << x_cos << "\n" + << y_cos << "\n" + << z << "\n" + << age << "\n" + << last_decay_x << "\n" + << last_decay_y << "\n" + << last_decay_z << "\n" + << last_decay_ID << "\n" + << Ekin_decay << "\n" + << last_interaction_x << "\n" + << last_interaction_y << "\n" + << last_interaction_z << "\n" + << last_interaction_ID << "\n" + << Ekin_parent << "\n" + << "--------------------------------------------------" << std::endl; + } + + gsimple_entry->metakey = metakey; + gsimple_entry->pdg = FLUKAtoPDG[FlukaID]; + gsimple_entry->wgt = weight; + gsimple_entry->vtxx = x * ShipUnit::cm / ShipUnit::m; // in m + gsimple_entry->vtxy = y * ShipUnit::cm / ShipUnit::m; + gsimple_entry->vtxz = z * ShipUnit::cm / ShipUnit::m; + gsimple_entry->dist = 0.; // Distance from hadron decay point to neutrino "vertex", to use for oscillations, + // for example. Don't use. + + // I'm assuming x_cos, y_cos are normalized. + double z_cos = sqrt(1 - pow(x_cos, 2) - pow(y_cos, 2)); + // And massless neutrinos + gsimple_entry->px = Ekin * x_cos; // in GeV/c + gsimple_entry->py = Ekin * y_cos; + gsimple_entry->pz = Ekin * z_cos; + gsimple_entry->E = Ekin; + + // Set auxiliary data + aux_entry->auxint.push_back(FlukaRun); + aux_entry->auxint.push_back(evtNumber); + aux_entry->auxint.push_back(FlukaID); + aux_entry->auxint.push_back(generationN); + aux_entry->auxint.push_back(last_decay_ID); + aux_entry->auxint.push_back(last_interaction_ID); + aux_entry->auxdbl.push_back(age); + aux_entry->auxdbl.push_back(weight); + aux_entry->auxdbl.push_back(Ekin); + aux_entry->auxdbl.push_back(x); + aux_entry->auxdbl.push_back(y); + aux_entry->auxdbl.push_back(z); + aux_entry->auxdbl.push_back(x_cos); + aux_entry->auxdbl.push_back(y_cos); + aux_entry->auxdbl.push_back(last_decay_x); + aux_entry->auxdbl.push_back(last_decay_y); + aux_entry->auxdbl.push_back(last_decay_z); + aux_entry->auxdbl.push_back(Ekin_decay); + aux_entry->auxdbl.push_back(last_interaction_x); + aux_entry->auxdbl.push_back(last_interaction_y); + aux_entry->auxdbl.push_back(last_interaction_z); + aux_entry->auxdbl.push_back(Ekin_parent); + + // Accumulate metadata + pdglist.insert(gsimple_entry->pdg); + min_weight = TMath::Min(min_weight, gsimple_entry->wgt); + max_weight = TMath::Max(max_weight, gsimple_entry->wgt); + max_energy = TMath::Max(max_energy, gsimple_entry->E); + + // All done! + tOut->Fill(); + } } - } - } - - // Sort out metadata - // Copy pdg list - meta_entry->pdglist.clear(); - for (std::set::const_iterator pdg_iterator = pdglist.begin(); - pdg_iterator != pdglist.end(); - ++pdg_iterator) meta_entry->pdglist.push_back(*pdg_iterator); - meta_entry->maxEnergy = max_energy; - meta_entry->maxWgt = max_weight; - meta_entry->minWgt = min_weight; - - meta_entry->protons = pp_collision_number; // Number of pp collisions. - for (int i = 0; i < 3; i++) meta_entry->windowBase[i] = plane_corner[i]*ShipUnit::cm/ShipUnit::m; - for (int i = 0; i < 3; i++) meta_entry->windowDir1[i] = plane_dir1[i]*ShipUnit::cm/ShipUnit::m; - for (int i = 0; i < 3; i++) meta_entry->windowDir2[i] = plane_dir2[i]*ShipUnit::cm/ShipUnit::m; - - meta_entry->infiles.push_back(inFileName); - meta_entry->seed = ran->GetSeed(); - meta_entry->metakey = metakey; - - meta_entry->auxintname.push_back("FlukaRun"); - meta_entry->auxintname.push_back("evtNumber"); - meta_entry->auxintname.push_back("FlukaID"); - meta_entry->auxintname.push_back("generationN"); - - meta_entry->auxdblname.push_back("age"); - meta_entry->auxdblname.push_back("FLUKA_weight"); - - metaOut->Fill(); - - fOut->cd(); - metaOut->Write(); - tOut->Write(); - fOut->Close(); - + } + + // plane corner and plane direction of the scoring plane (The scoring plane is tilted) + double plane_corner[] = {min_x, min_y, min_z}; + double plane_dir1[] = {max_x - min_x, 0, max_z - min_z}; + double plane_dir2[] = {0, max_y - min_y, max_z - min_z}; + + // Sort out metadata + // Copy pdg list + meta_entry->pdglist.clear(); + for (std::set::const_iterator pdg_iterator = pdglist.begin(); pdg_iterator != pdglist.end(); ++pdg_iterator) + meta_entry->pdglist.push_back(*pdg_iterator); + meta_entry->maxEnergy = max_energy; + meta_entry->maxWgt = max_weight; + meta_entry->minWgt = min_weight; + + meta_entry->protons = pp_collision_number; // Number of pp collisions. + for (int i = 0; i < 3; i++) + meta_entry->windowBase[i] = plane_corner[i] * ShipUnit::cm / ShipUnit::m; + for (int i = 0; i < 3; i++) + meta_entry->windowDir1[i] = plane_dir1[i] * ShipUnit::cm / ShipUnit::m; + for (int i = 0; i < 3; i++) + meta_entry->windowDir2[i] = plane_dir2[i] * ShipUnit::cm / ShipUnit::m; + + meta_entry->infiles.push_back(inFileName); + meta_entry->seed = ran->GetSeed(); + meta_entry->metakey = metakey; + + meta_entry->auxintname.push_back("FlukaRun"); + meta_entry->auxintname.push_back("evtNumber"); + meta_entry->auxintname.push_back("FlukaID"); + meta_entry->auxintname.push_back("generationN"); + meta_entry->auxintname.push_back("Last decay ID"); + meta_entry->auxintname.push_back("Last interaction ID"); + meta_entry->auxdblname.push_back("age"); + meta_entry->auxdblname.push_back("FLUKA_weight"); + meta_entry->auxdblname.push_back("Kinetic energy"); + meta_entry->auxdblname.push_back("X coordinate"); + meta_entry->auxdblname.push_back("Y coordinate"); + meta_entry->auxdblname.push_back("X cosine"); + meta_entry->auxdblname.push_back("Y cosine"); + meta_entry->auxdblname.push_back("Last decay X coordinate"); + meta_entry->auxdblname.push_back("Last decay Y coordinate"); + meta_entry->auxdblname.push_back("Last decay Z coordinate"); + meta_entry->auxdblname.push_back("Kinetic Energy of decay parent"); + meta_entry->auxdblname.push_back("Last interaction in Z coordinate"); + meta_entry->auxdblname.push_back("Last in interaction Y coordinate"); + meta_entry->auxdblname.push_back("Last in interaction Z coordinate"); + meta_entry->auxdblname.push_back("Kinetic energy of parent"); + + metaOut->Fill(); + + fOut->cd(); + metaOut->Write(); + tOut->Write(); + fOut->Close(); + std::cout << "Done converting" << std::endl; } - diff --git a/sndFairTasks/CMakeLists.txt b/sndFairTasks/CMakeLists.txt index dc730ad889..ef4d3f8bd8 100644 --- a/sndFairTasks/CMakeLists.txt +++ b/sndFairTasks/CMakeLists.txt @@ -29,6 +29,7 @@ set(SRCS DigiTaskSND.cxx ConvRawData.cxx boardMappingParser.cxx +MCEventBuilder.cxx ) Set(HEADERS) diff --git a/sndFairTasks/ConvRawData.cxx b/sndFairTasks/ConvRawData.cxx index 900b335d7d..c2b2fffeae 100644 --- a/sndFairTasks/ConvRawData.cxx +++ b/sndFairTasks/ConvRawData.cxx @@ -49,8 +49,8 @@ ConvRawData::ConvRawData() : FairTask("ConvRawData") , fEventTree(nullptr) , boards{} - , fEventHeader(nullptr) , fSNDLHCEventHeader(nullptr) + , fEventHeader(nullptr) , fDigiSciFi(nullptr) , fDigiMuFilter(nullptr) {} @@ -145,9 +145,6 @@ InitStatus ConvRawData::Init() timerBMap.Stop(); LOG (info) << "Time to set the board mapping " << timerBMap.RealTime(); - // Get the FairLogger - FairLogger *logger = FairLogger::GetLogger(); - eventNumber = fnStart; return kSUCCESS; @@ -336,7 +333,9 @@ void ConvRawData::Process0() system = MufiSystem[board_id][tofpet_id]; key = (tofpet_id%2)*1000 + tofpet_channel; tmp = boardMapsMu["MuFilter"][board.first][slots[tofpet_id]]; - if (debug || !tmp.find("not") == string::npos ) + // mini DTs are labelled DS 5V, board is 32 + if (tmp == "DS_5Vert") continue; + if (debug || !(tmp.find("not") == string::npos)) { LOG (info) << system << " " << key << " " << board.first << " " << tofpet_id << " " << tofpet_id%2 << " " << tofpet_channel; @@ -630,7 +629,9 @@ void ConvRawData::Process1() system = MufiSystem[board_id][tofpet_id]; key = (tofpet_id%2)*1000 + tofpet_channel; tmp = boardMapsMu["MuFilter"][board_name][slots[tofpet_id]]; - if (debug || !tmp.find("not") == string::npos ) + // mini DTs are labelled DS 5V, board is 32 + if (tmp == "DS_5Vert") continue; + if (debug || !(tmp.find("not") == string::npos)) { LOG (info) << system << " " << key << " " << board_name << " " << tofpet_id << " " << tofpet_id%2 << " " << tofpet_channel; @@ -778,7 +779,7 @@ void ConvRawData::StartTimeofRun(string Path) char *buff = new char[size]; status = file.Read(offset, size, buff, bytesRead); vector vec; - for (int i = 0; i < size; i++){vec.push_back(buff[i]);} + for (size_t i = 0; i < size; i++){vec.push_back(buff[i]);} j = json::parse(vec); status = file.Close(); delete info; @@ -826,7 +827,7 @@ void ConvRawData::DetMapping(string Path) char *buff = new char[size]; status = file.Read(offset, size, buff, bytesRead); vector vec; - for (int i = 0; i < size; i++){vec.push_back(buff[i]);} + for (size_t i = 0; i < size; i++){vec.push_back(buff[i]);} j = json::parse(vec); status = file.Close(); delete info; diff --git a/sndFairTasks/DigiTaskSND.cxx b/sndFairTasks/DigiTaskSND.cxx index 26d9bf5b82..df442c844f 100644 --- a/sndFairTasks/DigiTaskSND.cxx +++ b/sndFairTasks/DigiTaskSND.cxx @@ -62,6 +62,10 @@ InitStatus DigiTaskSND::Init() if ( fMCEventHeader == nullptr ) { fMCEventHeader = static_cast(gROOT->FindObjectAny("MCEventHeader.")); } + // To ensure the header is updated per newly read event. + ioman->GetInTree()->SetBranchAddress("MCEventHeader.", &fMCEventHeader); + LOG(INFO) << "MCEventHeader. branch is found"; + // Get input MC points fScifiPointArray = static_cast(ioman->GetObject("ScifiPoint")); fvetoPointArray = static_cast(ioman->GetObject("vetoPoint")); diff --git a/sndFairTasks/MCEventBuilder.cxx b/sndFairTasks/MCEventBuilder.cxx new file mode 100644 index 0000000000..b60403c2e0 --- /dev/null +++ b/sndFairTasks/MCEventBuilder.cxx @@ -0,0 +1,633 @@ +#include "MCEventBuilder.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "FairMCEventHeader.h" +#include "FairFileHeader.h" +#include "MuFilter.h" +#include "Scifi.h" +#include "FairLink.h" +#include "FairRunSim.h" +#include "FairRunAna.h" +#include "FairRootManager.h" +#include "MuFilterPoint.h" +#include "ScifiPoint.h" +#include "ShipMCTrack.h" +#include "ShipUnit.h" + +class vetoPoint; +class EmulsionDetPoint; + +namespace { + int n_clockcycles = 4; + + float t_electronics = 0.0; //if known + float timeWindow = n_clockcycles * (1 /(ShipUnit::snd_freq)) + t_electronics; + + MuFilter* MuFilterDet = nullptr; + float DsPropSpeed = 0.0; + float VandUpPropSpeed = 0.0; + + Scifi* ScifiDet = nullptr; + float ScifisignalSpeed = 0.0; + std::map>> siPMFibres; + + //Scifi Fast & Advanced Noise Filter Params + int min_scifi_boards = 3; + int min_us_boards = 5; + int min_ds_boards = 2; + + int min_veto_hits = 5; + int min_scifi_hits = 10; + int min_ds_hits =2; + + int min_veto_planes = 1; + int min_scifi_planes = 4; + int min_us_planes = 2; + int min_ds_planes = 2; +} + +MCEventBuilder::MCEventBuilder(const std::string& outputFileName, + bool saveFirst25nsOnly) + : FairTask("MCEventBuilder"), + fOutputFileName(outputFileName), + fSaveFirst25nsOnly(saveFirst25nsOnly), + fOutFile(nullptr), + fOutTree(nullptr), + fInMufiArray(nullptr), + fInSciFiArray(nullptr), + fInMCTrackArray(nullptr), + fInHeader(nullptr), + fOutMufiArray(nullptr), + fOutSciFiArray(nullptr), + fOutMCTrackArray(nullptr), + fOutHeader(nullptr) +{} + +MCEventBuilder::~MCEventBuilder() {} + +InitStatus MCEventBuilder::Init() { + LOG(INFO) << "Initializing MCEventBuilder"; + + FairRootManager* ioman = FairRootManager::Instance(); + if (!ioman) { + LOG(FATAL) << "MCEventBuilder::Init: RootManager not instantiated!"; + } + + // —---------- INPUT BRANCHES —----------- + // Only standard ROOT way of reading branches works for the header, otherwise nullptr + ioman->GetInTree()->SetBranchAddress("MCEventHeader.",&fInHeader); + fInMufiArray = static_cast(ioman->GetObject("MuFilterPoint")); + fInSciFiArray = static_cast(ioman->GetObject("ScifiPoint")); + fInMCTrackArray = static_cast(ioman->GetObject("MCTrack")); + + if (!fInMufiArray && !fInSciFiArray) { + LOG(ERROR) << "No Scifi and no MuFilter MC points array!"; + return kERROR; + } + + // —---------- OUTPUT FILE & TREE —---------- + fOutFile = new TFile(fOutputFileName.c_str(), "RECREATE"); + fOutTree = new TTree("cbmsim", "RebuiltEvents"); + fOutTree->SetTitle("/cbmroot_0"); + + fOutHeader = new FairMCEventHeader(); + fOutMufiArray = new TClonesArray("MuFilterPoint"); + fOutSciFiArray = new TClonesArray("ScifiPoint"); + fOutMCTrackArray = new TClonesArray("ShipMCTrack"); + fOutMufiArray->SetName("MuFilterPoint"); + fOutSciFiArray->SetName("ScifiPoint"); + fOutMCTrackArray->SetName("ShipMCTrack"); + + fOutTree->Branch("MuFilterPoint", &fOutMufiArray, 32000, 1); + fOutTree->Branch("ScifiPoint", &fOutSciFiArray, 32000, 1); + fOutTree->Branch("MCTrack", &fOutMCTrackArray, 32000, 1); + fOutTree->Branch("MCEventHeader.", &fOutHeader); + + // --------------Global Variables-------------------- + //Muon Filter + MuFilterDet = dynamic_cast(gROOT->GetListOfGlobals()->FindObject("MuFilter")); + if (!MuFilterDet) { + LOG(ERROR) << "MuFilter detector not found in gROOT"; + return kERROR; + } + DsPropSpeed = MuFilterDet->GetConfParF("MuFilter/DsPropSpeed"); + VandUpPropSpeed = MuFilterDet->GetConfParF("MuFilter/VandUpPropSpeed"); + + //Scifi + ScifiDet = dynamic_cast(gROOT->GetListOfGlobals()->FindObject("Scifi")); + if (!ScifiDet) { + LOG(ERROR) << "Scifi detector not found in gROOT"; + return kERROR; + } + ScifisignalSpeed = ScifiDet->GetConfParF("Scifi/signalSpeed"); + ScifiDet->SiPMmapping(); + siPMFibres = ScifiDet->GetFibresMap(); + + LOG(INFO) << "MCEventBuilder initialized successfully."; + return kSUCCESS; +} + +void MCEventBuilder::Exec(Option_t*) { + ProcessEvent(); +} + +std::vector MCEventBuilder::OrderedIds(const std::vector& times, + double firstTime) const { + size_t n = times.size(); + std::vector bins(n); + for (size_t i = 0; i < n; ++i) { + bins[i] = static_cast((times[i] - firstTime) / timeWindow); + } + + std::vector ids(n, 0); + for (size_t i = 1; i < n; ++i) { + if ((bins[i] - bins[i - 1]) < 1) { + ids[i] = ids[i - 1]; + } else { + ids[i] = ids[i - 1] + 1; + } + } + return ids; +} + +void MCEventBuilder::ProcessEvent() { + + fOutHeader->SetRunID(fInHeader->GetRunID()); + fOutHeader->SetEventID(fInHeader->GetEventID()); + fOutHeader->SetVertex(fInHeader->GetX(), fInHeader->GetY(), fInHeader->GetZ()); + fOutHeader->SetTime(fInHeader->GetT()); + fOutHeader->SetB(fInHeader->GetB()); + fOutHeader->SetNPrim(fInHeader->GetNPrim()); + fOutHeader->MarkSet(fInHeader->IsSet()); + fOutHeader->SetRotX(fInHeader->GetRotX()); + fOutHeader->SetRotY(fInHeader->GetRotY()); + fOutHeader->SetRotZ(fInHeader->GetRotZ()); + + //---------------------------Muon filter------------------------------------- + std::vector muFilterPoints; + std::vector muArrivalTimes; + std::vector muTrackIDs; + + for (auto* p : ROOT::RRangeCast(*fInMufiArray)) { + muFilterPoints.push_back(p); + muTrackIDs.push_back(p->GetTrackID()); + + int detID = p->GetDetectorID(); + + float propspeed; + if (floor(detID / 10000) == 3) + propspeed = DsPropSpeed; + else + propspeed = VandUpPropSpeed; + + TVector3 vLeft,vRight; + TVector3 impact(p->GetX(), p->GetY(), p->GetZ()); + MuFilterDet->GetPosition(detID, vLeft, vRight); + TVector3 vTop = vLeft; //Used for vertical bars + + //Vertical bars with only 1 readout at the top + if ( (floor(detID/10000)==3&&detID%1000>59) || + (floor(detID/10000)==1&&int(detID/1000)%10==2) ) { + double arrivalTime = p->GetTime() + (vTop - impact).Mag() / propspeed; + muArrivalTimes.push_back(arrivalTime); + } + //Horizontal + else{ + double tLeft = p->GetTime() + (vLeft - impact).Mag() / propspeed; + double tRight = p->GetTime() + (vRight - impact).Mag() / propspeed; + double arrivalTime = std::min(tLeft, tRight); + muArrivalTimes.push_back(arrivalTime); + } + } + + std::vector idxM(muArrivalTimes.size()); + std::iota(idxM.begin(), idxM.end(), 0); + std::sort(idxM.begin(), idxM.end(), [&](size_t a, size_t b) { + return muArrivalTimes[a] < muArrivalTimes[b]; + }); + + std::vector sortedMuPoints; + std::vector sortedMuArrivalTimes; + std::vector sortedMuTrackIDs; + + sortedMuPoints.reserve(muFilterPoints.size()); + sortedMuArrivalTimes.reserve(muArrivalTimes.size()); + sortedMuTrackIDs.reserve(muTrackIDs.size()); + + for (auto i : idxM) { + sortedMuPoints.push_back(muFilterPoints[i]); + sortedMuArrivalTimes.push_back(muArrivalTimes[i]); + sortedMuTrackIDs.push_back(muTrackIDs[i]); + } + + //---------------------------Scifi------------------------------------- + std::vector scifiPoints; + std::vector scifiArrivalTimes; + std::vector scifiTrackIDs; + + float signalSpeed = ScifisignalSpeed; + + for (auto* p : ROOT::RRangeCast(*fInSciFiArray)) { + scifiPoints.push_back(p); + scifiTrackIDs.push_back(p->GetTrackID()); + + TVector3 impact(p->GetX(), p->GetY(), p->GetZ()); + int point_detID = p->GetDetectorID(); + int localFiberID = (point_detID)%100000; + int a_sipmChan = static_cast(siPMFibres[localFiberID].begin()->first); + int detID_geo = int(point_detID/100000)*100000+a_sipmChan; + + TVector3 a, b; + ScifiDet->GetSiPMPosition(detID_geo, a, b); + bool verticalHit = int(detID_geo / 100000) % 10 == 1; + double distance; + if (verticalHit) { + distance = (b - impact).Mag(); + } else { + distance = (impact - a).Mag(); + } + double arrivalTime = p->GetTime() + distance / signalSpeed; + scifiArrivalTimes.push_back(arrivalTime); + } + + std::vector idxS(scifiArrivalTimes.size()); + std::iota(idxS.begin(), idxS.end(), 0); + std::sort(idxS.begin(), idxS.end(), [&](size_t a, size_t b) { + return scifiArrivalTimes[a] < scifiArrivalTimes[b]; + }); + + std::vector sortedScifiPoints; + std::vector sortedScifiArrivalTimes; + std::vector sortedScifiTrackIDs; + + sortedScifiPoints.reserve(scifiPoints.size()); + sortedScifiArrivalTimes.reserve(scifiArrivalTimes.size()); + sortedScifiTrackIDs.reserve(scifiTrackIDs.size()); + + for (auto i : idxS) { + sortedScifiPoints.push_back(scifiPoints[i]); + sortedScifiArrivalTimes.push_back(scifiArrivalTimes[i]); + sortedScifiTrackIDs.push_back(scifiTrackIDs[i]); + } + + //----------------------------------Tracks------------------------------------- + std::vector mcTrackClones; + for (auto* t : ROOT::RRangeCast(*fInMCTrackArray)) { + mcTrackClones.push_back(t); + } + + //---------Finding the earliest time between Scifi and MuFilter----------------- + double tMufi = sortedMuArrivalTimes.empty() ? -1 : sortedMuArrivalTimes.front(); + double tScifi = sortedScifiArrivalTimes.empty() ? -1 : sortedScifiArrivalTimes.front(); + bool hasMCPoints = (tMufi >= 0 || tScifi >= 0); + double firstT = hasMCPoints ? (tMufi < 0 ? tScifi : (tScifi < 0 ? tMufi : std::min(tMufi, tScifi))): 0; + + if (!hasMCPoints) { + fOutMufiArray->Delete(); + fOutSciFiArray->Delete(); + fOutMCTrackArray->Delete(); + fOutTree->Fill(); + return; + } + + //------------------Preparations before chunking the data--------------------- + auto idsMufi = OrderedIds(sortedMuArrivalTimes, firstT); + auto idsScifi = OrderedIds(sortedScifiArrivalTimes, firstT); + + bool FirstEvent = true; + + std::vector used; + int i_mufi = 0, i_scifi = 0, sliceMufi = 0, sliceScifi = 0; + + while (i_mufi < (int)sortedMuArrivalTimes.size() || i_scifi < (int)sortedScifiArrivalTimes.size()) { + // To avoid re-copy of points or tracks from chunk N to following chunks, one needs to free the memory + // e.g. chunk N has points with track ids(indices) T-10, T-5, T, chunk N+1 holds points with track id T+1. + // If memory is not de-allocated, all tracks of smaller index than T+1 will be also written to chunk N+1. + // This will be incorrect since there are no points of tracks with id < T in the N+1 chunck. + fOutMufiArray->Delete(); + fOutSciFiArray->Delete(); + fOutMCTrackArray->Delete(); + + std::vector muSlicePoints; + std::vector scifiSlicePoints; + fOutMCTrackArray->ExpandCreate(mcTrackClones.size()); + + //Adding the mother track to the first event + if (sliceMufi==0 && sliceScifi==0){ + ShipMCTrack* newTrack = new ((*fOutMCTrackArray)[0]) ShipMCTrack(*mcTrackClones[0]); + // In NC neutrino events, the outgoing neutrino is also saved in the first chunk, + // otherwise it will be lost. This is to facilitate NC event tagging in analysis. + if (mcTrackClones.size()>1 && std::set({12,14,16}).count(fabs(mcTrackClones[1]->GetPdgCode()))){ + ShipMCTrack* newTrack = new ((*fOutMCTrackArray)[1]) ShipMCTrack(*mcTrackClones[1]); + } + } + + //MUON FILTER POINTS CHUNKING + while (i_mufi < (int)sortedMuArrivalTimes.size() && idsMufi[i_mufi] == sliceMufi) { + MuFilterPoint* origMu = sortedMuPoints[i_mufi]; + muSlicePoints.push_back(origMu); + Int_t trackID = origMu->GetTrackID(); + Int_t detID = origMu->GetDetectorID(); + TVector3 pos; origMu->Position(pos); + TVector3 mom; origMu->Momentum(mom); + Double_t time = origMu->GetTime(); + Double_t length = origMu->GetLength(); + Double_t eLoss = origMu->GetEnergyLoss(); + Int_t pdgCode = origMu->PdgCode(); + + new ((*fOutMufiArray)[fOutMufiArray->GetEntriesFast()]) + MuFilterPoint(trackID, detID, pos, mom, time, length, eLoss, pdgCode); + + int tid = sortedMuTrackIDs[i_mufi++]; + //MC tracks having ID=-2 are below the pre-set energy threshold (default is 100MeV) and cannot be linked to their corresponding MC point. Such tracks are not saved to the chunked events. + if (tid != -2) { + ShipMCTrack* newTrack = new ((*fOutMCTrackArray)[tid]) ShipMCTrack(*mcTrackClones[tid]); + } + } + ++sliceMufi; + + //SCIFI POINTS CHUNKING + while (i_scifi < (int)sortedScifiArrivalTimes.size() && idsScifi[i_scifi] == sliceScifi) { + ScifiPoint* origSci = sortedScifiPoints[i_scifi]; + scifiSlicePoints.push_back(origSci); + Int_t trackID = origSci->GetTrackID(); + Int_t detID = origSci->GetDetectorID(); + TVector3 pos; origSci->Position(pos); + TVector3 mom; origSci->Momentum(mom); + Double_t time = origSci->GetTime(); + Double_t length = origSci->GetLength(); + Double_t eLoss = origSci->GetEnergyLoss(); + Int_t pdgCode = origSci->PdgCode(); + + new ((*fOutSciFiArray)[fOutSciFiArray->GetEntriesFast()]) + ScifiPoint(trackID, detID, pos, mom, time, length, eLoss, pdgCode); + + int tid = sortedScifiTrackIDs[i_scifi++]; + //MC tracks having ID=-2 are below the pre-set energy threshold (default is 100MeV) and cannot be linked to their corresponding MC point. Such tracks are not saved to the chunked events. + if (tid != -2) { + ShipMCTrack* newTrack = new ((*fOutMCTrackArray)[tid]) ShipMCTrack(*mcTrackClones[tid]); + } + } + ++sliceScifi; + + //Noise Filters + if (!(FastNoiseFilterScifi_Hits(fOutSciFiArray, siPMFibres) || FastNoiseFilterMu_Hits(fOutMufiArray))) { + fOutMufiArray->Delete(); + fOutSciFiArray->Delete(); + fOutMCTrackArray->Delete(); + FirstEvent = false; + continue; + } + else if(!(FastNoiseFilterScifi_Boards(fOutSciFiArray, siPMFibres) || FastNoiseFilterMu_Boards(fOutMufiArray))){ + fOutMufiArray->Delete(); + fOutSciFiArray->Delete(); + fOutMCTrackArray->Delete(); + FirstEvent = false; + continue; + } + + if (!(AdvancedNoiseFilterScifi(fOutSciFiArray, siPMFibres) || AdvancedNoiseFilterMu(fOutMufiArray))){ + fOutMufiArray->Delete(); + fOutSciFiArray->Delete(); + fOutMCTrackArray->Delete(); + FirstEvent = false; + continue; + } + + //--------Save only first 25ns chunks mode------- + if (fSaveFirst25nsOnly) { + if (FirstEvent) { + fOutTree->Fill(); + FirstEvent = false; + break; + } + } + //-----Save all chunks mode------- + else { + fOutTree->Fill(); + } + } +} + +//-----------------------Mu Fast Noise Filtering------------------------------- +bool MCEventBuilder::FastNoiseFilterMu_Hits(TClonesArray* muArray) { + std::set veto, ds; + for (int i = 0; i < muArray->GetEntriesFast(); ++i) { + auto* p = static_cast(muArray->At(i)); + int detID = p->GetDetectorID(); + int system = detID / 10000; + if (system == 1) { + veto.insert(detID); + if (veto.size() >= min_veto_hits) { + return true; + } + } + else if (system == 3) { + ds.insert(detID); + if (ds.size() >= min_ds_hits) { + return true; + } + } + } + return false; +} + +bool MCEventBuilder::FastNoiseFilterMu_Boards(TClonesArray* muArray) { + //The digits used here to assign IDs to systems are set to match the DAQ board IDs from the 2022 board_mapping.json file + std::set us, ds; + for (int i = 0; i < muArray->GetEntriesFast(); ++i) { + auto* p = static_cast(muArray->At(i)); + int detID = p->GetDetectorID(); + int system = detID / 10000; + int TypeofPlane = (detID / 1000) % 10; + if (system == 2) { + if (TypeofPlane == 1 || TypeofPlane == 2) { + us.insert(7); + } else if (TypeofPlane == 3 || TypeofPlane == 4) { + us.insert(60); + } else if (TypeofPlane == 5) { + us.insert(52); + } + if (us.size() >= min_us_boards) { + return true; + } + } else if (system == 3) { + if (TypeofPlane == 1) { + ds.insert(41); + } else if (TypeofPlane == 2) { + ds.insert(35); + } else if (TypeofPlane == 3 || TypeofPlane == 4) { + ds.insert(55); + } + if (ds.size() >= min_ds_boards) { + return true; + } + } + } + return false; +} + +//-----------------------Mu Advanced Noise Filtering------------------------------- +bool MCEventBuilder::AdvancedNoiseFilterMu(TClonesArray* muArray) { + std::set veto_planes, us_planes, ds_planes; + for (int i = 0; i < muArray->GetEntriesFast(); ++i) { + auto* p = static_cast(muArray->At(i)); + int detID = p->GetDetectorID(); + int system = detID / 10000; + if (system == 1) { + veto_planes.insert(detID / 1000); + if (veto_planes.size() >= min_veto_planes) { + return true; + } + } else if (system == 2) { + us_planes.insert(detID / 1000); + if (us_planes.size() >= min_us_planes) { + return true; + } + } else if (system == 3) { + ds_planes.insert(detID / 1000); + if (ds_planes.size() >= min_ds_planes) { + return true; + } + } + } + return false; +} + + +//-----------------------Scifi Fast Noise Filtering------------------------------- +bool MCEventBuilder::FastNoiseFilterScifi_Hits( + TClonesArray* scifiArray, + const std::map>>& siPMFibresMap) + { + std::set detIDs; + for (int i = 0; i < scifiArray->GetEntriesFast(); ++i) { + auto* p = static_cast(scifiArray->At(i)); + int point_detID = p->GetDetectorID(); + int locFibreID = point_detID % 100000; + + for (const auto& sipmChan : siPMFibresMap.at(locFibreID)) { + int channel = sipmChan.first; + int detID_geo = (point_detID / 100000) * 100000 + channel; + detIDs.insert(detID_geo); + + if (detIDs.size() >= min_scifi_hits){ + return true; + } + } + } + return false; + } + +bool MCEventBuilder::FastNoiseFilterScifi_Boards( + TClonesArray* scifiArray, + const std::map>>& siPMFibresMap) + { + //A DAQ board reads one Scifi mat and we use the first three digits of the detID_geo, representing the unique combination of StationPlaneMat, to count boards with fired channels. + std::set dividedDetIds; + + for (int i = 0; i < scifiArray->GetEntriesFast(); ++i) { + auto* p = static_cast(scifiArray->At(i)); + int point_detID = p->GetDetectorID(); + int locFibreID = point_detID % 100000; + + for (const auto& sipmChan : siPMFibresMap.at(locFibreID)) { + int channel = sipmChan.first; + int detID_geo = (point_detID / 100000) * 100000 + channel; + dividedDetIds.insert(detID_geo / 10000); + + if (dividedDetIds.size() >= min_scifi_boards){ + return true; + } + } + } + return false; + } + +//-----------------------Scifi Advanced Noise Filtering------------------------------- +bool MCEventBuilder::AdvancedNoiseFilterScifi( + TClonesArray* scifiArray, + const std::map>>& siPMFibresMap) + { + std::set UniquePlanes; + for (int i = 0; i < scifiArray->GetEntriesFast(); ++i) { + auto* p = static_cast(scifiArray->At(i)); + int point_detID = p->GetDetectorID(); + int locFibreID = point_detID % 100000; + + for (const auto& sipmChan : siPMFibresMap.at(locFibreID)) { + int channel = sipmChan.first; + int detID_geo = (point_detID / 100000) * 100000 + channel; + UniquePlanes.insert(detID_geo / 100000); + + if (UniquePlanes.size() >= min_scifi_planes) { + return true; + } + } + } + return false; + } + +void MCEventBuilder::FinishTask() { + // Make the output file structure compatible with FairTasks to facilitate + // runing the digitization task. + fOutFile->cd(); + TFolder* folder = new TFolder("cbmroot", "Main Folder"); + TFolder* fol_Stack = folder->AddFolder("Stack", "Stack Folder"); + auto mcTrackArr = new TClonesArray("ShipMCTrack"); + mcTrackArr->SetName("MCTrack"); // avoid default plural names + fol_Stack->Add(mcTrackArr); + TFolder* fol_veto = folder->AddFolder("veto", "veto Folder"); + auto vetoArr = new TClonesArray("vetoPoint"); + vetoArr->SetName("vetoPoint"); + fol_veto->Add(vetoArr); + TFolder* fol_EmulsionDet = folder->AddFolder("EmulsionDet", "EmulsionDetector Folder"); + auto emuArr = new TClonesArray("EmulsionDetPoint"); + emuArr->SetName("EmulsionDetPoint"); + fol_EmulsionDet->Add(emuArr); + TFolder* fol_Scifi = folder->AddFolder("Scifi", "Scifi Folder"); + auto scifiArr = new TClonesArray("ScifiPoint"); + scifiArr->SetName("ScifiPoint"); + fol_Scifi->Add(scifiArr); + TFolder* fol_MuFilter = folder->AddFolder("MuFilter", "MuFilter Folder"); + auto mufArr = new TClonesArray("MuFilterPoint"); + mufArr->SetName("MuFilterPoint"); + fol_MuFilter->Add(mufArr); + TFolder* fol_Event = folder->AddFolder("Event", "Event Folder"); + folder->Write(); + + fOutFile->cd(); + TList* branch_list = new TList(); + branch_list->SetName("BranchList"); + branch_list->Add(new TObjString("MCTrack")); + branch_list->Add(new TObjString("vetoPoint")); + branch_list->Add(new TObjString("EmulsionDetPoint")); + branch_list->Add(new TObjString("ScifiPoint")); + branch_list->Add(new TObjString("MuFilterPoint")); + branch_list->Add(new TObjString("MCEventHeader.")); + fOutFile->WriteObject(branch_list, "BranchList"); + auto* timebased_branch_list = new TList(); + timebased_branch_list->SetName("TimeBasedBranchList"); + fOutFile->WriteObject(timebased_branch_list, "TimeBasedBranchList"); + fOutFile->cd(); + auto* fileHeader = new FairFileHeader(); + fileHeader->SetTitle("FileHeader"); + fileHeader->Write("FileHeader"); + fOutFile->cd(); + if (fOutTree) { + fOutTree->Write(); + } + LOG(INFO) << "Writing and closing output file: " << fOutputFileName; +} diff --git a/sndFairTasks/MCEventBuilder.h b/sndFairTasks/MCEventBuilder.h new file mode 100644 index 0000000000..0d4069db5f --- /dev/null +++ b/sndFairTasks/MCEventBuilder.h @@ -0,0 +1,73 @@ +#ifndef MCEVENTBUILDER_H +#define MCEVENTBUILDER_H +#include +#include "FairTask.h" +#include +#include +#include +#include "FairMCPoint.h" +#include +#include + +class TFile; +class TTree; +class MuFilterPoint; +class ScifiPoint; +class ShipMCTrack; +class FairMCEventHeader; + +class MCEventBuilder : public FairTask { +public: + MCEventBuilder(const std::string& outputFileName, bool saveOnlyFirst25); + ~MCEventBuilder(); + + virtual InitStatus Init(); + virtual void Exec(Option_t* opt); + virtual void FinishTask(); + +private: + //Function I need later for ordering the mc points + std::vector OrderedIds(const std::vector& times, double firstTime) const; + + //Function that does all the processing + void ProcessEvent(); + + // Fast filter functions + bool FastNoiseFilterMu_Hits(TClonesArray* muArray); + bool FastNoiseFilterMu_Boards(TClonesArray* muArray); + + bool FastNoiseFilterScifi_Hits( + TClonesArray* scifiArray, + const std::map>>& siPMFibres); + bool FastNoiseFilterScifi_Boards( + TClonesArray* scifiArray, + const std::map>>& siPMFibres); + + //Advanced Noise Filter + bool AdvancedNoiseFilterScifi( + TClonesArray* scifiArray, + const std::map>>& siPMFibres); + + bool AdvancedNoiseFilterMu(TClonesArray* muArray); + + //Input + bool fSaveFirst25nsOnly; + FairMCEventHeader* fInHeader; + TClonesArray* fInMufiArray; + TClonesArray* fInSciFiArray; + TClonesArray* fInMCTrackArray; + + //Output + std::string fOutputFileName; + TFile* fOutFile; + TTree* fOutTree; + FairMCEventHeader* fOutHeader; + TClonesArray* fOutMufiArray; + TClonesArray* fOutSciFiArray; + TClonesArray* fOutMCTrackArray; + + + ClassDef(MCEventBuilder, 1) +}; + +#endif // MCEVENTBUILDER_H diff --git a/sndFairTasks/sndFairTasksLinkDef.h b/sndFairTasks/sndFairTasksLinkDef.h index 59168d5240..7ea592291b 100755 --- a/sndFairTasks/sndFairTasksLinkDef.h +++ b/sndFairTasks/sndFairTasksLinkDef.h @@ -6,6 +6,7 @@ #pragma link C++ class DigiTaskSND; #pragma link C++ class ConvRawData; +#pragma link C++ class MCEventBuilder; #endif