snd::analysis_tools Namespace Reference

Functions
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)
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), float time_upper_range=2E9/(ShipUnit::snd_freq/ShipUnit::hertz), bool make_selection=true)
std::unique_ptr< TClonesArray >	selectScifiHits (const TClonesArray &digiHits, int station, bool orientation, const std::map< std::string, float > &selection_parameters, bool make_selection=true)
std::unique_ptr< TClonesArray >	filterScifiHits (const TClonesArray &digiHits, const std::map< std::string, float > &selection_parameters, int method=0, std::string setup="TI18")
std::unique_ptr< TClonesArray >	filterScifiHits (const TClonesArray &digiHits, int method=0, std::string setup="TI18")
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")

Function Documentation

calculateSiPMNumber()

int snd::analysis_tools::calculateSiPMNumber

(

int

reference_SiPM

)

Definition at line 361 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 413 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 374 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;
}

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"
	)

Definition at line 287 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.";
   }
 
   sndScifiHit *hit;
   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\", "
                       "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 = selectScifiHits(digiHits, station, orientation, selection_parameters, true);
            for (auto *p : *supportArray) {
               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"
	)

Definition at line 339 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"] = 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);
}

findScifiStation() [1/2]

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

Definition at line 93 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 77 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;
}

getFractionalHitsPerScifiPlane() [1/2]

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

(

const TClonesArray *

digiHits

)

Definition at line 67 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 53 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;
}

getScifiHits()

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

Definition at line 158 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 12 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]++;
         }
      }
   }
}

getTotalSciFiHits() [1/2]

int snd::analysis_tools::getTotalSciFiHits

(

const TClonesArray *

digiHits

)

Definition at line 42 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 36 of file sndSciFiTools.cxx.

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

peakScifiTiming()

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

Definition at line 124 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);
 
   int refStation = ((sndScifiHit *)digiHits.At(0))->GetStation();
   bool refOrientation = ((sndScifiHit *)digiHits.At(0))->isVertical();
   float hitTime = -1.0;
 
   for (auto *p : digiHits) {
      auto *hit = dynamic_cast<sndScifiHit *>(p);
      if (!validateHit(hit, refStation, refOrientation)) {
         continue;
      }
      hitTime = hit->GetTime() * 1E9 / (ShipUnit::snd_freq / ShipUnit::hertz);
      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
	)

Definition at line 259 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);
}

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),
		float	time_upper_range = 2E9/(ShipUnit::snd_freq/ShipUnit::hertz),
		bool	make_selection = true
	)

Definition at line 182 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;
 
   if (make_selection) {
 
      auto selectedHits = getScifiHits(digiHits, station, orientation);
 
      peakTiming = peakScifiTiming(*selectedHits, bins_x, min_x, max_x);
 
      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() * 1E9 / (ShipUnit::snd_freq / ShipUnit::hertz)) ||
             (hit->GetTime() * 1E9 / (ShipUnit::snd_freq / ShipUnit::hertz) > 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);
 
      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() * 1E9 / (ShipUnit::snd_freq / ShipUnit::hertz)) ||
             (hit->GetTime() * 1E9 / (ShipUnit::snd_freq / ShipUnit::hertz) > 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 476 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 533 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);
}