snd::analysis_tools Namespace Reference

Classes
class	ScifiPlane
class	USPlane

Functions
std::string	GetGeoPath (int run_number, std::string csv_file_path="")
std::pair< Scifi *, MuFilter * >	GetGeometry (const std::string &geometry_path)
std::pair< Scifi *, MuFilter * >	GetGeometry (int run_number, const std::string &csv_file_path="")
std::vector< ScifiPlane >	FillScifi (const Configuration &configuration, TClonesArray *sf_hits, Scifi *scifi_geometry)
std::vector< USPlane >	FillUS (const Configuration &configuration, TClonesArray *mufi_hits, MuFilter *mufilter_geometry, bool use_small_sipms=false)
void	getSciFiHitsPerStation (const TClonesArray *digiHits, std::vector< int > &horizontal_hits, std::vector< int > &vertical_hits)
int	getTotalSciFiHits (std::vector< int > &horizontal_hits, std::vector< int > &vertical_hits)
int	getTotalSciFiHits (const TClonesArray *digiHits)
std::vector< float >	getFractionalHitsPerScifiPlane (std::vector< int > &horizontal_hits, std::vector< int > &vertical_hits)
std::vector< float >	getFractionalHitsPerScifiPlane (const TClonesArray *digiHits)
int	findScifiStation (std::vector< int > &horizontal_hits, std::vector< int > &vertical_hits, float threshold)
int	findScifiStation (const TClonesArray *digiHits, float threshold)
float	peakScifiTiming (const TClonesArray &digiHits, int bins, float min_x, float max_x, bool isMC=false)
std::unique_ptr< TClonesArray >	getScifiHits (const TClonesArray &digiHits, int station, bool orientation)
std::unique_ptr< TClonesArray >	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< TClonesArray >	selectScifiHits (const TClonesArray &digiHits, int station, bool orientation, const std::map< std::string, float > &selection_parameters, bool make_selection=true, bool isMC=false)
std::unique_ptr< TClonesArray >	filterScifiHits (const TClonesArray &digiHits, const std::map< std::string, float > &selection_parameters, int method=0, std::string setup="TI18", bool isMC=false)
std::unique_ptr< TClonesArray >	filterScifiHits (const TClonesArray &digiHits, int method=0, std::string setup="TI18", bool isMC=false)
int	calculateSiPMNumber (int reference_SiPM)
int	densityScifi (int reference_SiPM, const TClonesArray &digiHits, int radius, int min_hit_density, bool min_check)
bool	densityCheck (const TClonesArray &digiHits, int radius=64, int min_hit_density=36, int station=1, bool orientation=false)
int	showerInteractionWall (const TClonesArray &digiHits, const std::map< std::string, float > &selection_parameters, int method=0, std::string setup="TI18")
int	showerInteractionWall (const TClonesArray &digiHits, int method=0, std::string setup="TI18")
std::pair< double, double >	findCentreOfGravityPerStation (const TClonesArray *digiHits, int station, Scifi *ScifiDet)
std::pair< std::vector< double >, std::vector< double > >	hitPositionVectorsPerStation (const TClonesArray *digiHits, int station, Scifi *ScifiDet)
std::pair< double, double >	hitDensityPerStation (const TClonesArray *digiHits, int station, Scifi *ScifiDet)
int	GetScifiShowerStart (const std::vector< ScifiPlane > &scifi_planes)
int	GetScifiShowerEnd (const std::vector< ScifiPlane > &scifi_planes)
int	GetUSShowerStart (const std::vector< USPlane > &us_planes)
int	GetUSShowerEnd (const std::vector< USPlane > &us_planes)
std::pair< ROOT::Math::XYZPoint, ROOT::Math::XYZVector >	GetShowerInterceptAndDirection (const Configuration &configuration, const std::vector< ScifiPlane > &scifi_planes, const std::vector< USPlane > &us_planes)
std::pair< std::vector< ScifiPlane >, std::vector< USPlane > >	GetShoweringPlanes (const std::vector< ScifiPlane > &scifi_planes, const std::vector< USPlane > &us_planes)
std::pair< double, double >	GetSciFiSpatialAnisotropy (const std::vector< ScifiPlane > &scifi_planes, bool use_all_centroids=false)
double	GetUSSpatialAnisotropy (const std::vector< USPlane > &us_planes, bool use_all_centroids=false)
std::unique_ptr< TChain >	GetTChain (int run_number, int n_files=-1, const std::string &csv_file_path="")
std::unique_ptr< TChain >	GetTChain (const std::string &file_name)
std::string	GetDataBasePath (int run_number, std::string csv_file_path="")

Function Documentation

calculateSiPMNumber()

int snd::analysis_tools::calculateSiPMNumber

(

int

reference_SiPM

)

Definition at line 381 of file sndSciFiTools.cxx.

{
 
   int ref_matAux = reference_SiPM % 100000;
   int ref_mat = ref_matAux / 10000;
   int ref_arrayAux = reference_SiPM % 10000;
   int ref_array = ref_arrayAux / 1000;
   int ref_channel = reference_SiPM % 1000;
   int referenceChannel = ref_mat * 4 * 128 + ref_array * 128 + ref_channel;
 
   return referenceChannel;
}

densityCheck()

bool snd::analysis_tools::densityCheck	(	const TClonesArray &	digiHits,
		int	radius = 64,
		int	min_hit_density = 36,
		int	station = 1,
		bool	orientation = false
	)

Definition at line 433 of file sndSciFiTools.cxx.

{
 
   if (digiHits.GetEntries() <= 0) {
      return false;
   }
 
   if (radius <= 0) {
      LOG(FATAL) << "Radius<=0. Please provide a radius bigger than 0.";
      return false;
   }
 
   if (min_hit_density < 0) {
      LOG(FATAL) << "Min_hit_density < 0. Please provide a min_hit_density >= 0.";
      return false;
   }
 
   if (min_hit_density > 2 * radius) {
      LOG(warning) << "Warning! Radius of density check does not allow for the required minimum density!";
      return false;
   }
 
   // Creating a vector that stores the fired SiPM channels (should already be ordered but we sort
   // afterwards to make sure)
   std::vector<int> fired_channels{};
 
   // Only fill the set with hits from the same station and with the same orientation as given in
   // argument (false==horizontal and true==vertical)
   for (auto *p : digiHits) {
      auto *hit = dynamic_cast<sndScifiHit *>(p);
      if (!validateHit(hit, station, orientation)) {
         continue;
      }
      fired_channels.push_back(calculateSiPMNumber(hit->GetChannelID()));
   }
 
   int n_fired_channels = fired_channels.size();
 
   if (n_fired_channels < min_hit_density) {
      return false;
   }
 
   // Looping over the ordered hits, checking whether within an interval of "min_hit_density" the
   // difference between the channel IDs is smaller than "radius". If so, then we have the required
   // density. Else, we check the next combination until we meet the criterion, or end the loop,
   // returning false
   std::sort(fired_channels.begin(), fired_channels.end());
   for (int i = 0; i < n_fired_channels - min_hit_density; i++) {
 
      if (fired_channels[i + min_hit_density - 1] - fired_channels[i] <= radius * 2) {
         return true;
      }
   }
   return false;
}

densityScifi()

int snd::analysis_tools::densityScifi	(	int	reference_SiPM,
		const TClonesArray &	digiHits,
		int	radius,
		int	min_hit_density,
		bool	min_check
	)

Definition at line 394 of file sndSciFiTools.cxx.

{
 
   int hit_density = 0;
 
   bool orientation = false;
   if (int(reference_SiPM / 100000) % 10 == 1) {
      orientation = true;
   }
   int ref_station = reference_SiPM / 1000000;
   int referenceChannel = calculateSiPMNumber(reference_SiPM);
 
   for (auto *p : digiHits) {
      auto *hit = dynamic_cast<sndScifiHit *>(p);
      if (!validateHit(hit, ref_station, orientation)) {
         continue;
      }
      int hitChannel = calculateSiPMNumber(hit->GetChannelID());
      if (radius == -1) {
         hit_density++;
      } else {
         if (hitChannel > referenceChannel + radius) {
            break;
         }
         if (abs(referenceChannel - hitChannel) <= radius) {
            hit_density++;
         }
      }
 
      if (min_check && (hit_density >= min_hit_density)) {
         break;
      }
   }
 
   return hit_density;
}

FillScifi()

std::vector< snd::analysis_tools::ScifiPlane > snd::analysis_tools::FillScifi	(	const Configuration &	configuration,
		TClonesArray *	sf_hits,
		Scifi *	scifi_geometry
	)

Definition at line 14 of file sndPlaneTools.cxx.

{
 
  std::vector<snd::analysis_tools::ScifiPlane> scifi_planes;
  int n_sf_hits{sf_hits->GetEntries()};
 
  const int max_station = configuration.scifi_n_stations;
  std::vector<std::vector<sndScifiHit*>> stations_hits(max_station);
 
  for (int i{0}; i < n_sf_hits; ++i) {
      auto hit = static_cast<sndScifiHit*>(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;
}

FillUS()

std::vector< snd::analysis_tools::USPlane > snd::analysis_tools::FillUS	(	const Configuration &	configuration,
		TClonesArray *	mufi_hits,
		MuFilter *	mufilter_geometry,
		bool	use_small_sipms = false
	)

Definition at line 39 of file sndPlaneTools.cxx.

{
 
  std::vector<snd::analysis_tools::USPlane> us_planes;
  int n_mufi_hits{mufi_hits->GetEntries()};
 
  const int n_station = configuration.us_n_stations;
  std::vector<std::vector<MuFilterHit*>> plane_hits(n_station);
 
  for (int i{0}; i < n_mufi_hits; ++i) {
    auto hit = static_cast<MuFilterHit*>(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;
}

filterScifiHits() [1/2]

std::unique_ptr< TClonesArray > snd::analysis_tools::filterScifiHits	(	const TClonesArray &	digiHits,
		const std::map< std::string, float > &	selection_parameters,
		int	method = 0,
		std::string	setup = "TI18",
		bool	isMC = false
	)

Definition at line 308 of file sndSciFiTools.cxx.

{
   TClonesArray supportArray("sndScifiHit", 0);
   auto filteredHits = std::make_unique<TClonesArray>("sndScifiHit", digiHits.GetEntries());
   int filteredHitsIndex = 0;
   int ScifiStations = 5;
   if (setup == "H8") {
      ScifiStations = 4;
   } else {
      LOG(info) << "\"TI18\" setup will be used by default, please provide \"H8\" for the Testbeam setup.";
   }
 
   TIter hitIterator(&supportArray);
 
   if (method == 0) {
 
      if ((selection_parameters.find("bins_x") == selection_parameters.end()) ||
          (selection_parameters.find("min_x") == selection_parameters.end()) ||
          (selection_parameters.find("max_x") == selection_parameters.end()) ||
          (selection_parameters.find("time_lower_range") == selection_parameters.end())) {
         LOG(FATAL) << "In order to use method 0 please provide the correct selection_parameters. Consider the default "
                       "= {{\"bins_x\", 52.}, {\"min_x\", 0.}, {\"max_x\", 26.}, {\"time_lower_range\", "
                       "1E9/(2*ShipUnit::snd_freq/ShipUnit::hertz)}, {\"time_upper_range\", "
                       "1.2E9/(ShipUnit::snd_freq/ShipUnit::hertz)}}";
      }
 
      // This is overwriting the previous arrays with the newest one
      for (auto station : ROOT::MakeSeq(1, ScifiStations + 1)) {
         for (auto orientation : {false, true}) {
 
            auto supportArray_p = selectScifiHits(digiHits, station, orientation, selection_parameters, true, isMC);
            for (auto *p : *supportArray_p) {
               auto *hit = dynamic_cast<sndScifiHit *>(p);
               if (hit->isValid()) {
                  new ((*filteredHits)[filteredHitsIndex++]) sndScifiHit(*hit);
               }
            }
         }
      }
   } else {
      LOG(error) << "Please provide a valid time filter method from:";
      LOG(error) << "(0): Events within \\mp time_lower_range time_upper_range of the peak of the time distribution "
                    "for Scifi Hits within each station and orientation";
      LOG(FATAL) << "selection_parameters = {bins_x, min_x, max_x, time_lower_range, time_upper_range}.";
   }
   return filteredHits;
}

filterScifiHits() [2/2]

std::unique_ptr< TClonesArray > snd::analysis_tools::filterScifiHits	(	const TClonesArray &	digiHits,
		int	method = 0,
		std::string	setup = "TI18",
		bool	isMC = false
	)

Definition at line 359 of file sndSciFiTools.cxx.

{
 
   std::map<std::string, float> selection_parameters;
 
   if (method == 0) {
 
      selection_parameters["bins_x"] = 52.0;
      selection_parameters["min_x"] = 0.0;
      selection_parameters["max_x"] = 26.0;
      selection_parameters["time_lower_range"] = 1E9 / (2 * ShipUnit::snd_freq / ShipUnit::hertz);
      selection_parameters["time_upper_range"] = 1.2E9 / (ShipUnit::snd_freq / ShipUnit::hertz);
 
   } else {
      LOG(FATAL) << "Please use method=0. No other methods implemented so far.";
   }
 
   return filterScifiHits(digiHits, selection_parameters, method, setup, isMC);
}

findCentreOfGravityPerStation()

std::pair< double, double > snd::analysis_tools::findCentreOfGravityPerStation	(	const TClonesArray *	digiHits,
		int	station,
		Scifi *	ScifiDet
	)

Definition at line 573 of file sndSciFiTools.cxx.

{
   if (!digiHits) {
      LOG(ERROR) << "Error: digiHits is null in findCentreOfGravityPerStation";
   }
   
   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};
}

findScifiStation() [1/2]

int snd::analysis_tools::findScifiStation	(	const TClonesArray *	digiHits,
		float	threshold
	)

Definition at line 96 of file sndSciFiTools.cxx.

{
   std::vector<int> horizontal_hits = std::vector<int>(5);
   std::vector<int> vertical_hits = std::vector<int>(5);
 
   getSciFiHitsPerStation(digiHits, horizontal_hits, vertical_hits);
 
   return findScifiStation(horizontal_hits, vertical_hits, threshold);
}

findScifiStation() [2/2]

int snd::analysis_tools::findScifiStation	(	std::vector< int > &	horizontal_hits,
		std::vector< int > &	vertical_hits,
		float	threshold
	)

Definition at line 80 of file sndSciFiTools.cxx.

{
 
   std::vector<float> frac = getFractionalHitsPerScifiPlane(horizontal_hits, vertical_hits);
 
   std::vector<float> frac_sum = std::vector<float>(frac.size());
 
   std::partial_sum(frac.begin(), frac.end(), frac_sum.begin());
 
   std::vector<float>::iterator station =
      std::find_if(frac_sum.begin(), frac_sum.end(), [&threshold](const auto &f) { return f > threshold; });
 
   return station - frac_sum.begin() + 1;
}

GetDataBasePath()

std::string snd::analysis_tools::GetDataBasePath	(	int	run_number,
		std::string	csv_file_path = ""
	)

Definition at line 13 of file sndTchainGetter.cxx.

                                                                                    {
    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);
    }
 
    std::string line;
    std::getline(file, line); // skip header
 
    while (std::getline(file, line)) {
        std::istringstream ss(line);
        std::string token;
 
        std::getline(ss, token, ',');
        int min_run = std::stoi(token);
 
        std::getline(ss, token, ',');
        int max_run = std::stoi(token);
 
        std::getline(ss, token);
        std::string path = token;
 
        if (run_number >= min_run && run_number <= max_run) {
            return path;
        }
    }
    throw std::runtime_error("Run number not found in CSV mapping.");
}

getFractionalHitsPerScifiPlane() [1/2]

std::vector< float > snd::analysis_tools::getFractionalHitsPerScifiPlane

(

const TClonesArray *

digiHits

)

Definition at line 70 of file sndSciFiTools.cxx.

{
   std::vector<int> horizontal_hits = std::vector<int>(5);
   std::vector<int> vertical_hits = std::vector<int>(5);
 
   getSciFiHitsPerStation(digiHits, horizontal_hits, vertical_hits);
 
   return getFractionalHitsPerScifiPlane(horizontal_hits, vertical_hits);
}

getFractionalHitsPerScifiPlane() [2/2]

std::vector< float > snd::analysis_tools::getFractionalHitsPerScifiPlane	(	std::vector< int > &	horizontal_hits,
		std::vector< int > &	vertical_hits
	)

Definition at line 56 of file sndSciFiTools.cxx.

{
 
   int total_hits = getTotalSciFiHits(horizontal_hits, vertical_hits);
 
   std::vector<float> fractional_hits_per_station = std::vector<float>(horizontal_hits.size());
 
   std::transform(horizontal_hits.begin(), horizontal_hits.end(), vertical_hits.begin(),
                  fractional_hits_per_station.begin(),
                  [&total_hits](const auto &hor, const auto &ver) { return ((float)hor + ver) / total_hits; });
 
   return fractional_hits_per_station;
}

GetGeometry() [1/2]

std::pair< Scifi *, MuFilter * > snd::analysis_tools::GetGeometry

(

const std::string &

geometry_path

)

Definition at line 49 of file sndGeometryGetter.cxx.

{
    TPython::Exec("import SndlhcGeo");
    TPython::Exec(("SndlhcGeo.GeoInterface('" + geometry_path + "')").c_str());
 
    // Init detectors
    Scifi *scifi = new Scifi("Scifi", kTRUE);
    MuFilter *mufilter = new MuFilter("MuFilter", kTRUE);
 
    // Retrieve the detectors from ROOT's global list
    scifi = dynamic_cast<Scifi *>(gROOT->GetListOfGlobals()->FindObject("Scifi"));
    mufilter = dynamic_cast<MuFilter *>(gROOT->GetListOfGlobals()->FindObject("MuFilter"));
 
    return std::make_pair(scifi, mufilter);
}

GetGeometry() [2/2]

std::pair< Scifi *, MuFilter * > snd::analysis_tools::GetGeometry	(	int	run_number,
		const std::string &	csv_file_path = ""
	)

Definition at line 66 of file sndGeometryGetter.cxx.

{
    std::string geometry_path = GetGeoPath(run_number, csv_file_path);
 
    return GetGeometry(geometry_path);
}

GetGeoPath()

std::string snd::analysis_tools::GetGeoPath	(	int	run_number,
		std::string	csv_file_path = ""
	)

Definition at line 15 of file sndGeometryGetter.cxx.

{
    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);
    }
 
    std::string line;
    std::getline(file, line); // skip header
 
    while (std::getline(file, line)) {
        std::istringstream ss(line);
        std::string token;
 
        std::getline(ss, token, ',');
        int min_run = std::stoi(token);
 
        std::getline(ss, token, ',');
        int max_run = std::stoi(token);
 
        std::getline(ss, token);
        std::string path = token;
 
        if (run_number >= min_run && run_number <= max_run) {
            return path;
        }
    }
    throw std::runtime_error("Run number not found in CSV mapping.");
}

getScifiHits()

std::unique_ptr< TClonesArray > snd::analysis_tools::getScifiHits	(	const TClonesArray &	digiHits,
		int	station,
		bool	orientation
	)

Definition at line 171 of file sndSciFiTools.cxx.

{
 
   auto selectedHits = std::make_unique<TClonesArray>("sndScifiHit");
 
   int i = 0;
   for (auto *p : digiHits) {
      auto *hit = dynamic_cast<sndScifiHit *>(p);
      if (!validateHit(hit, station, orientation)) {
         continue;
      }
      new ((*selectedHits)[i++]) sndScifiHit(*hit);
   }
 
   return selectedHits;
}

getSciFiHitsPerStation()

void snd::analysis_tools::getSciFiHitsPerStation	(	const TClonesArray *	digiHits,
		std::vector< int > &	horizontal_hits,
		std::vector< int > &	vertical_hits
	)

Definition at line 15 of file sndSciFiTools.cxx.

{
 
   // Clear hits per plane vectors
   std::fill(horizontal_hits.begin(), horizontal_hits.end(), 0);
   std::fill(vertical_hits.begin(), vertical_hits.end(), 0);
 
   // Add valid hits to hits per plane vectors
   sndScifiHit *hit;
   TIter hitIterator(digiHits);
 
   while ((hit = dynamic_cast<sndScifiHit *>(hitIterator.Next()))) {
      if (hit->isValid()) {
         int station = hit->GetStation();
         if (hit->isVertical()) {
            vertical_hits[station - 1]++;
         } else {
            horizontal_hits[station - 1]++;
         }
      }
   }
}

GetScifiShowerEnd()

int snd::analysis_tools::GetScifiShowerEnd

(

const std::vector< ScifiPlane > &

scifi_planes

)

GetScifiShowerStart()

int snd::analysis_tools::GetScifiShowerStart

(

const std::vector< ScifiPlane > &

scifi_planes

)

GetSciFiSpatialAnisotropy()

std::pair< double, double > snd::analysis_tools::GetSciFiSpatialAnisotropy	(	const std::vector< ScifiPlane > &	scifi_planes,
		bool	use_all_centroids = false
	)

Definition at line 209 of file sndShowerTools.cxx.

{
  std::vector<double> 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};
}

GetShoweringPlanes()

std::pair< std::vector< ScifiPlane >, std::vector< USPlane > > snd::analysis_tools::GetShoweringPlanes	(	const std::vector< ScifiPlane > &	scifi_planes,
		const std::vector< USPlane > &	us_planes
	)

GetShowerInterceptAndDirection()

std::pair< ROOT::Math::XYZPoint, ROOT::Math::XYZVector > snd::analysis_tools::GetShowerInterceptAndDirection	(	const Configuration &	configuration,
		const std::vector< ScifiPlane > &	scifi_planes,
		const std::vector< USPlane > &	us_planes
	)

GetTChain() [1/2]

std::unique_ptr< TChain > snd::analysis_tools::GetTChain

(

const std::string &

file_name

)

Definition at line 59 of file sndTchainGetter.cxx.

                                                                            {
    auto tchain = std::make_unique<TChain>("rawConv");
    tchain->Add(file_name.c_str());
    return tchain;
}

GetTChain() [2/2]

std::unique_ptr< TChain > snd::analysis_tools::GetTChain	(	int	run_number,
		int	n_files = -1,
		const std::string &	csv_file_path = ""
	)

Definition at line 45 of file sndTchainGetter.cxx.

                                                                                                             {
    std::string base_folder = GetDataBasePath(run_number, csv_file_path);
    auto tchain = std::make_unique<TChain>("rawConv");
    if (n_files == -1) {
        tchain->Add(Form("%srun_%06d/sndsw_raw-*", base_folder.c_str(), run_number));
    }
    else {
        for (int i = 0; i<n_files; ++i){
            tchain->Add(Form("%srun_%06d/sndsw_raw-%04d.root", base_folder.c_str(), run_number, i));
        }
    }
    return tchain;
};

getTotalSciFiHits() [1/2]

int snd::analysis_tools::getTotalSciFiHits

(

const TClonesArray *

digiHits

)

Definition at line 45 of file sndSciFiTools.cxx.

{
   std::vector<int> horizontal_hits = std::vector<int>(5);
   std::vector<int> vertical_hits = std::vector<int>(5);
 
   getSciFiHitsPerStation(digiHits, horizontal_hits, vertical_hits);
 
   return getTotalSciFiHits(horizontal_hits, vertical_hits);
}

getTotalSciFiHits() [2/2]

int snd::analysis_tools::getTotalSciFiHits	(	std::vector< int > &	horizontal_hits,
		std::vector< int > &	vertical_hits
	)

Definition at line 39 of file sndSciFiTools.cxx.

{
   return std::accumulate(horizontal_hits.begin(), horizontal_hits.end(),
                          std::accumulate(vertical_hits.begin(), vertical_hits.end(), 0));
}

GetUSShowerEnd()

int snd::analysis_tools::GetUSShowerEnd

(

const std::vector< USPlane > &

us_planes

)

GetUSShowerStart()

int snd::analysis_tools::GetUSShowerStart

(

const std::vector< USPlane > &

us_planes

)

GetUSSpatialAnisotropy()

double snd::analysis_tools::GetUSSpatialAnisotropy	(	const std::vector< USPlane > &	us_planes,
		bool	use_all_centroids = false
	)

Definition at line 252 of file sndShowerTools.cxx.

{
  std::vector<double> 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;
}

hitDensityPerStation()

std::pair< double, double > snd::analysis_tools::hitDensityPerStation	(	const TClonesArray *	digiHits,
		int	station,
		Scifi *	ScifiDet
	)

Definition at line 706 of file sndSciFiTools.cxx.

                                                                                                                        {
    /* 
     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<double>, std::vector<double>> hit_position_vec = hitPositionVectorsPerStation(digiHits, station, ScifiDet);
    std::vector<double> hit_position_x = hit_position_vec.first; // Assuming we are interested in X positions
    std::vector<double> 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};
}

hitPositionVectorsPerStation()

std::pair< std::vector< double >, std::vector< double > > snd::analysis_tools::hitPositionVectorsPerStation	(	const TClonesArray *	digiHits,
		int	station,
		Scifi *	ScifiDet
	)

Definition at line 592 of file sndSciFiTools.cxx.

                                                                                                                                                      {
     /*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<double> x_positions;
    std::vector<double> y_positions;
 
    TVector3 A, B;
    for (auto* hit : ROOT::RRangeCast<sndScifiHit*, false, decltype(*digiHits)>(*digiHits)) {
      if (!hit || !hit->isValid()) {
         continue;
      }
      if (hit->GetStation() != station) {
         continue;
      }
      ScifiDet->GetSiPMPosition(hit->GetDetectorID(), A, B);
      if (hit->isVertical()) {
         x_positions.push_back((A.X() + B.X()) * 0.5);
      }
      else {
         y_positions.push_back((A.Y() + B.Y()) * 0.5);
      }
    }
    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};
}

peakScifiTiming()

float snd::analysis_tools::peakScifiTiming	(	const TClonesArray &	digiHits,
		int	bins,
		float	min_x,
		float	max_x,
		bool	isMC = false
	)

Definition at line 127 of file sndSciFiTools.cxx.

{
 
   if (digiHits.GetEntries() <= 0) {
      LOG(warning) << "digiHits has no valid SciFi Hits and as such no maximum for the timing distribution.";
      return -1.;
   }
 
   TH1F ScifiTiming("Timing", "Scifi Timing", bins, min_x, max_x);
 
   Scifi *ScifiDet = dynamic_cast<Scifi*> (gROOT->GetListOfGlobals()->FindObject("Scifi") );
   auto* first_hit = static_cast<sndScifiHit*>(digiHits[0]);
   int refStation = first_hit->GetStation();
   bool refOrientation = first_hit->isVertical();
   float hitTime = -1.0;
   float timeConversion = 1.;
   if (!isMC) {
      timeConversion = 1E9 / (ShipUnit::snd_freq / ShipUnit::hertz);
   }
 
   for (auto *p : digiHits) {
      auto *hit = dynamic_cast<sndScifiHit *>(p);
      if (!validateHit(hit, refStation, refOrientation)) {
         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;
      }
      ScifiTiming.Fill(hitTime);
      hitTime = -1.0;
   }
 
   float peakTiming = (ScifiTiming.GetMaximumBin() - 0.5) * (max_x - min_x) / bins + min_x;
 
   return peakTiming;
}

selectScifiHits() [1/2]

std::unique_ptr< TClonesArray > snd::analysis_tools::selectScifiHits	(	const TClonesArray &	digiHits,
		int	station,
		bool	orientation,
		const std::map< std::string, float > &	selection_parameters,
		bool	make_selection = true,
		bool	isMC = false
	)

Definition at line 278 of file sndSciFiTools.cxx.

{
 
   if ((selection_parameters.find("bins_x") == selection_parameters.end()) ||
       (selection_parameters.find("min_x") == selection_parameters.end()) ||
       (selection_parameters.find("max_x") == selection_parameters.end()) ||
       (selection_parameters.find("time_lower_range") == selection_parameters.end())) {
      LOG(FATAL) << "In order to use method 0 please provide the correct selection_parameters. Consider the default = "
                    "{{\"bins_x\", 52.}, {\"min_x\", 0.}, {\"max_x\", 26.}, {\"time_lower_range\", "
                    "1E9/(2*ShipUnit::snd_freq/ShipUnit::hertz))}, {\"time_upper_range\", "
                    "2E9/(ShipUnit::snd_freq/ShipUnit::hertz)}}";
   }
 
   float time_upper_range = -1.;
 
   if (selection_parameters.find("time_upper_range") == selection_parameters.end()) {
      time_upper_range = selection_parameters.at("time_lower_range");
   } else {
      time_upper_range = selection_parameters.at("time_upper_range");
   }
 
   return selectScifiHits(digiHits, station, orientation, int(selection_parameters.at("bins_x")),
                          selection_parameters.at("min_x"), selection_parameters.at("max_x"),
                          selection_parameters.at("time_lower_range"), time_upper_range, make_selection,
              isMC);
}

selectScifiHits() [2/2]

std::unique_ptr< TClonesArray > snd::analysis_tools::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
	)

Definition at line 195 of file sndSciFiTools.cxx.

{
 
   if (bins_x < 1) {
      LOG(FATAL) << "bins_x in selection_parameters cannot be <1. Consider using the default value of 52 instead.";
   }
   if (min_x > max_x) {
      LOG(warning) << "In selection_parameters min_x > max_x. Values will be swapped.";
      float aux_float = min_x;
      min_x = max_x;
      max_x = aux_float;
   }
   if (min_x < 0.0) {
      LOG(warning) << "In selection_parameters min_x < 0.0. Consider using the default value of 0.0.";
   }
   if (max_x < 0.0) {
      LOG(FATAL) << "In selection_parameters max_x < 0.0. Consider using the default value of 26.0.";
   }
   if (time_lower_range <= 0.0) {
      LOG(FATAL) << "In selection_parameters time_lower_range <= 0.0. Value should always be positive, in ns.";
   }
   if (time_upper_range <= 0.0) {
      LOG(FATAL) << "In selection_parameters time_upper_range <= 0.0. Value should always be positive, in ns.";
   }
 
   auto filteredHits = std::make_unique<TClonesArray>("sndScifiHit", digiHits.GetEntries());
 
   float peakTiming = -1.0;
 
   float timeConversion = 1.;
   if (!isMC) {
      timeConversion = 1E9 / (ShipUnit::snd_freq / ShipUnit::hertz);
   }
 
   if (make_selection) {
 
      auto selectedHits = getScifiHits(digiHits, station, orientation);
 
      peakTiming = peakScifiTiming(*selectedHits, bins_x, min_x, max_x, isMC);
 
      int i = 0;
      for (auto *p : *selectedHits) {
         auto *hit = dynamic_cast<sndScifiHit *>(p);
         if (!validateHit(hit, station, orientation)) {
            continue;
         }
         if ((peakTiming - time_lower_range > hit->GetTime() * timeConversion) ||
             (hit->GetTime() * timeConversion > peakTiming + time_upper_range)) {
            continue;
         }
         new ((*filteredHits)[i++]) sndScifiHit(*hit);
      }
 
   } else {
      // Does not create selectedHits and just uses digiHits (not unique_ptr)
 
      peakTiming = peakScifiTiming(digiHits, bins_x, min_x, max_x, isMC);
 
      int i = 0;
      for (auto *p : digiHits) {
         auto *hit = dynamic_cast<sndScifiHit *>(p);
         if (!validateHit(hit, station, orientation)) {
            continue;
         }
         if ((peakTiming - time_lower_range > hit->GetTime() * timeConversion) ||
             (hit->GetTime() * timeConversion > peakTiming + time_upper_range)) {
            continue;
         }
         new ((*filteredHits)[i++]) sndScifiHit(*hit);
      }
   }
 
   return filteredHits;
}

showerInteractionWall() [1/2]

int snd::analysis_tools::showerInteractionWall	(	const TClonesArray &	digiHits,
		const std::map< std::string, float > &	selection_parameters,
		int	method = 0,
		std::string	setup = "TI18"
	)

Definition at line 496 of file sndSciFiTools.cxx.

{
 
   int totalScifiStations = 5;
   if (setup == "H8") {
      totalScifiStations = 4;
   } else {
      LOG(info) << "\"TI18\" setup will be used by default, please provide \"H8\" for the Testbeam setup.";
   }
 
   // There is always 1 more Scifi Station than a target block. As such, showerStart == totalScifiStations
   // means that the shower did not start developing in the target, before the last Scifi Station
   int showerStart = totalScifiStations;
 
   if (method == 0) {
 
      if ((selection_parameters.find("radius") == selection_parameters.end()) ||
          (selection_parameters.find("min_hit_density") == selection_parameters.end())) {
         LOG(FATAL)
            << "Argument of select_parameters is incorrect. Please provide a map with the arguments \"radius\" and "
               "\"min_hit_density\". Consider using the default values of {{\"radius\",64}, {\"min_hit_density\",36}}.";
      }
 
      for (int scifiStation = 1; scifiStation <= totalScifiStations; scifiStation++) {
 
         // For each ScifiStation we check whether we see clusters with the desired parameters
         // By default, we require passing the check on both the horizontal and vertical mats
         // In case selection_parameters["orientation"] was provided, we set the NOT CHOSEN orientation
         // as true by default, so as to only need to pass the chosen orientation check
         bool horizontalCheck = false;
         bool verticalCheck = false;
         if (selection_parameters.find("orientation") != selection_parameters.end()) {
            if (int(selection_parameters.at("orientation")) == 0) {
               verticalCheck = true;
            }
            if (int(selection_parameters.at("orientation")) == 1) {
               horizontalCheck = true;
            }
         }
         horizontalCheck = (densityCheck(digiHits, int(selection_parameters.at("radius")),
                                         int(selection_parameters.at("min_hit_density")), scifiStation, false) ||
                            horizontalCheck);
         verticalCheck = (densityCheck(digiHits, int(selection_parameters.at("radius")),
                                       int(selection_parameters.at("min_hit_density")), scifiStation, true) ||
                          verticalCheck);
         if ((horizontalCheck) && (verticalCheck)) {
            showerStart = scifiStation - 1;
            return showerStart;
         }
      }
   }
 
   return showerStart;
}

showerInteractionWall() [2/2]

int snd::analysis_tools::showerInteractionWall	(	const TClonesArray &	digiHits,
		int	method = 0,
		std::string	setup = "TI18"
	)

Definition at line 553 of file sndSciFiTools.cxx.

{
 
   if (method != 0) {
      LOG(FATAL) << "Please use method=0. No other methods implemented so far.";
   }
 
   std::map<std::string, float> selection_parameters = {{"radius", 64.}, {"min_hit_density", 36.}};
 
   return showerInteractionWall(digiHits, selection_parameters, method, setup);
}