Merge pull request #15 from Kinexity/develJZ

another pull request

DataFormats/CMakeLists.txt

@@ -25,6 +25,7 @@ ROOT_GENERATE_DICTIONARY(G__${MODULE_NAME}
   ${PROJECT_SOURCE_DIR}/${MODULE_NAME}/include/GeometryTPC.h 
   ${PROJECT_SOURCE_DIR}/${MODULE_NAME}/include/Geometry_Strip.h  
   ${PROJECT_SOURCE_DIR}/${MODULE_NAME}/include/SigClusterTPC.h  
+  ${PROJECT_SOURCE_DIR}/${MODULE_NAME}/include/Event_Information.h  
   LINKDEF ${PROJECT_SOURCE_DIR}/${MODULE_NAME}/include/LinkDef.h)
 
 # The data is just added to the executable, because in some IDEs (QtCreator) 

DataFormats/include/CommonDefinitions.h

@@ -25,7 +25,7 @@ class RV_Storage {
 	};
 };
 
-const double pi = 4 * atan(1);
+const double pi = 4 * std::atan(1);
 const double deg_to_rad = pi / 180.0;
 
 enum class direction : unsigned {

DataFormats/include/EventTPC.h

@@ -57,13 +57,13 @@ class EventTPC {
 
 	// helper methods for extracting data points
 	// they return 0.0 for non-existing data points
-	double GetValByStrip(direction strip_dir, int strip_number, int time_cell);  // valid range [0-2][1-1024][0-511]
+	double GetValByStrip(direction strip_dir, int strip_number, int time_cell) const;  // valid range [0-2][1-1024][0-511]
 
 	double GetMaxCharge();                   // maximal charge from all strips
 
 	std::shared_ptr<SigClusterTPC> GetOneCluster(double thr, int delta_strips, int delta_timecells); // applies clustering threshold to all space-time data points 
 
-	std::shared_ptr<TH2D> GetStripVsTime(direction strip_dir);                               // whole event, all strip dirs
+	TH2D&& GetStripVsTime(direction strip_dir);                               // whole event, all strip dirs
 };
 
 #endif

DataFormats/include/GeometryTPC.h

@@ -59,7 +59,7 @@ class GeometryTPC {
   std::map<std::tuple<int, int, int, int>, std::shared_ptr<Geometry_Strip>> arrayByAget;/*channel_idx[0-63] AGET_idx[0-3] ASAD_idx[0-3] COBO_idx[0-1]*/
   std::map<std::tuple<direction, int> ,std::shared_ptr<Geometry_Strip>> stripArray; // key = {U,V,W}, STRIP_NUMBER [1-1024]
   std::vector<int> ASAD_N;       // pair=(COBO_idx, number of ASAD boards)
-  std::vector<int> FPN_chanId;     // FPN channels in AGET chips
+  const std::vector<int> FPN_chanId = { 11, 22, 45, 56 };     // FPN channels in AGET chips
   double pad_size;                 // in [mm]
   double pad_pitch;                // in [mm]
   double strip_pitch;              // in [mm]

DataFormats/include/LinkDef.h

@@ -8,5 +8,7 @@
 #pragma link C++ class GeometryTPC+;
 #pragma link C++ class Geometry_Strip+;
 #pragma link C++ class SigClusterTPC+;
+#pragma link C++ class Event_Information+;
+#pragma link C++ class property<long>+;
 
 #endif

DataFormats/include/SigClusterTPC.h

@@ -35,12 +35,6 @@ class SigClusterTPC {
 	std::set< std::tuple<int, direction, int>> hitListByTimeDir; // same list of selected space-time cells as a map
 														  // with key=(REBIN_TIME_CELL [0-511], STRIP_DIR [0-2]) 
 														  // that returns vector of STRIP_NUMBERS 
-	std::set<std::tuple<direction, int, int>> envelope_hitList; // temporary list of selected envelope space-time cells for further analysis where return
-								 // value=key(STRIP_DIR [0-2], STRIP_NUM [1-1024], REBIN_TIME_CELL [0-511])
-	std::set< std::tuple<int, direction, int>> envelope_hitListByTimeDir; // same temporary list of selected envelope space-time cells as a map
-														  // with key=(REBIN_TIME_CELL [0-511], STRIP_DIR [0-2]) 
-														  // that returns vector of STRIP_NUMBERS 
-
 	// statistics variables
 	std::map<direction, size_t> nhits;   // number of space-time cells in a given U,V,W direction
 

DataFormats/src/EventTPC.cpp

@@ -10,11 +10,7 @@ void EventTPC::Clear() {
 bool EventTPC::AddValByStrip(std::shared_ptr<Geometry_Strip> strip, int time_cell, double val) {  // valid range [0-2][1-1024][0-511]
 	auto op = (*strip)();
 	if (time_cell < 0 || time_cell >= Geometry().GetAgetNtimecells()) return false;
-	auto& ref = chargeMap[{op.dir, op.num, time_cell}];
-	ref += val; //update hit or add new one
-
-	// upadate charge maxima
-	glb_max_charge = std::max(glb_max_charge, ref);
+	chargeMap[{op.dir, op.num, time_cell}] += val; //update hit or add new one
 
 	return true;
 }
@@ -23,7 +19,7 @@ bool EventTPC::AddValByAgetChannel(int cobo_idx, int asad_idx, int aget_idx, int
 	return AddValByStrip(Geometry().GetStripByAget(cobo_idx, asad_idx, aget_idx, channel_idx), time_cell, val);
 }
 
-double EventTPC::GetValByStrip(direction strip_dir, int strip_number, int time_cell) {  // valid range [0-2][1-1024][0-511]
+double EventTPC::GetValByStrip(direction strip_dir, int strip_number, int time_cell) const {  // valid range [0-2][1-1024][0-511]
 	// check if hit is unique
 	auto it = chargeMap.find({ strip_dir,strip_number,time_cell });
 	if (it != chargeMap.end()) {
@@ -87,25 +83,25 @@ std::shared_ptr<SigClusterTPC> EventTPC::GetOneCluster(double thr, int delta_str
 	return cluster;
 }
 
-std::shared_ptr<TH2D> EventTPC::GetStripVsTime(direction strip_dir) {  // valid range [0-2]
+TH2D&& EventTPC::GetStripVsTime(direction strip_dir) {  // valid range [0-2]
 
-	std::shared_ptr<TH2D> result = std::make_shared<TH2D>(Form("hraw_%s_vs_time_evt%lld", Geometry().GetDirName(strip_dir).c_str(), event_info.EventId()),
+	TH2D result{ Form("hraw_%s_vs_time_evt%lld", Geometry().GetDirName(strip_dir).c_str(), event_info.EventId()),
 		Form("Event-%lld: Raw signals from %s strips;Time bin [arb.u.];%s strip no.;Charge/bin [arb.u.]",
 			event_info.EventId(), Geometry().GetDirName(strip_dir).c_str(), Geometry().GetDirName(strip_dir).c_str()),
 		Geometry().GetAgetNtimecells(),
 		0.0 - 0.5,
 		1. * Geometry().GetAgetNtimecells() - 0.5, // ends at 511.5 (cells numbered from 0 to 511)
 		Geometry().GetDirNstrips(strip_dir),
 		1.0 - 0.5,
-		1. * Geometry().GetDirNstrips(strip_dir) + 0.5);
+		1. * Geometry().GetDirNstrips(strip_dir) + 0.5 };
 	// fill new histogram
 	auto min_strip = chargeMap.lower_bound({strip_dir, std::numeric_limits<int>::min(), std::numeric_limits<int>::min()});
 	auto max_strip = chargeMap.upper_bound({strip_dir, std::numeric_limits<int>::max(), std::numeric_limits<int>::max()});
 	for (auto charge = min_strip; charge != max_strip; charge++) {
 		auto strip_num = std::get<1>(charge->first);
 		auto time_cell = std::get<2>(charge->first);
 		auto val = charge->second;
-		result->Fill(time_cell, strip_num, val);
+		result.Fill(time_cell, strip_num, val);
 	}
-	return result;
+	return std::move(result);
 }

DataFormats/src/GeometryTPC.cpp

@@ -114,20 +114,14 @@ GeometryTPC::GeometryTPC(std::string  fname, bool debug)
 	trigger_delay(0.),
 	drift_zmin(0.),
 	drift_zmax(0.),
-	_debug(debug)
+	_debug(debug),
+	reference_point(0., 0.) // default REFERENCE POINT position on XY plane
 {
 	if (_debug) {
 		std::cout << "GeometryTPC::Constructor - Started..." << std::endl;
 	}
 
-	// default REFERENCE POINT position on XY plane
-	reference_point.Set(0., 0.);
-
 	// FPN channels in each AGET chip
-	FPN_chanId.push_back(11);
-	FPN_chanId.push_back(22);
-	FPN_chanId.push_back(45);
-	FPN_chanId.push_back(56);
 	AGET_Nchan_fpn = FPN_chanId.size();
 	AGET_Nchan_raw = AGET_Nchan + AGET_Nchan_fpn;
 
@@ -180,20 +174,16 @@ bool GeometryTPC::Load(std::string fname) {
 				for (auto dir : dir_vec_UVW) {
 					angle[dir] = fmod(fmod(angle[dir], 360.) + 360., 360.); // convert to range [0, 360 deg[ range 
 				}
-				if (fmod(angle[direction::U], 180.) != fmod(angle[direction::V], 180.) &&  // reject parallel / anti-parallel duplicates 
-					fmod(angle[direction::V], 180.) != fmod(angle[direction::W], 180.) &&  // reject parallel / anti-parallel duplicates 
-					fmod(angle[direction::W], 180.) != fmod(angle[direction::U], 180.)) { // reject parallel / anti-parallel duplicates 
-					std::cout << Form("Angle of U/V/W strips wrt X axis = %lf / %lf / %lf deg", angle[direction::U], angle[direction::V], angle[direction::W]) << std::endl;
-				}
-				else {
+				std::cout << Form("Angle of U/V/W strips wrt X axis = %lf / %lf / %lf deg", angle[direction::U], angle[direction::V], angle[direction::W]) << std::endl;
+				if (fmod(angle[direction::U], 180.) == fmod(angle[direction::V], 180.) ||  // reject parallel / anti-parallel duplicates 
+					fmod(angle[direction::V], 180.) == fmod(angle[direction::W], 180.) ||  // reject parallel / anti-parallel duplicates 
+					fmod(angle[direction::W], 180.) == fmod(angle[direction::U], 180.)) { // reject parallel / anti-parallel duplicates 
 					std::cerr << "ERROR: Wrong U/V/W angles !!!" << std::endl;
-					std::cout << Form("Angle of U/V/W strips wrt X axis = %lf / %lf / %lf deg", angle[direction::U], angle[direction::V], angle[direction::W]) << std::endl;
 					return false;
 				}
 			}
 			else {
-				std::cerr << "ERROR: Wrong U/V/W angles !!!" << std::endl;
-				std::cout << Form("Angle of U/V/W strips wrt X axis = %lf / %lf / %lf deg", angle[direction::U], angle[direction::V], angle[direction::W]) << std::endl;
+				std::cerr << "ERROR: U/V/W angles not provided !!!" << std::endl;
 				return false;
 			}
 		}
@@ -258,7 +248,7 @@ bool GeometryTPC::Load(std::string fname) {
 
 						// update maximal ASAD index (by: COBO board) 
 						if (cobo >= ASAD_N.size()) { //resize ASAD_N if necessary
-							ASAD_N.resize(cobo + 1);
+							ASAD_N.resize(static_cast<size_t>(cobo) + 1);
 						}
 						ASAD_N[cobo] = std::max(ASAD_N[cobo], asad + 1); // ASAD indexing starts from 0
 
@@ -318,8 +308,6 @@ bool GeometryTPC::Load(std::string fname) {
 	std::cout << "Number of ASAD boards = " << GetAsadNboards() << std::endl;
 	std::cout << "Number of COBO boards = " << GetCoboNboards() << std::endl;
 
-	// now initialize TH2Poly (while initOK=true) and set initOK flag according to the result 
-
 	std::cout << "\n==== INITIALIZING TPC GEOMETRY - END ====\n\n";
 
 	if (_debug) {
@@ -358,7 +346,7 @@ int GeometryTPC::Aget_fpn2raw(int FPN_idx) { // valid range [0-3]
 }
 
 int GeometryTPC::Global_normal2normal(int COBO_idx, int ASAD_idx, int aget_idx, int channel_idx) {   // valid range [0-1][0-3][0-3][0-63]
-	if (COBO_idx < 0 || COBO_idx >= COBO_N || ASAD_idx < 0 || ASAD_idx > ASAD_N[COBO_idx] || aget_idx < 0 || aget_idx >= AGET_Nchips ||
+	if (ASAD_idx < 0 || ASAD_idx > ASAD_N[COBO_idx] || aget_idx < 0 || aget_idx >= AGET_Nchips ||
 		channel_idx < 0 || channel_idx >= AGET_Nchan) return ERROR;
 	int ASAD_offset = std::accumulate(ASAD_N.begin(), ASAD_N.begin() + COBO_idx, 0);
 	return ((ASAD_offset + ASAD_idx) * AGET_Nchips + aget_idx) * AGET_Nchan + channel_idx;
@@ -396,8 +384,8 @@ bool GeometryTPC::MatchCrossPoint(std::array<int, 3> strip_nums, double radius,
 	std::transform(strip_nums.begin(), strip_nums.end(), dir_vec_UVW.begin(), strips.begin(), [&](int strip_num, direction dir) { return GetStripByDir(dir, strip_num); });
 	const double rad2 = pow(radius, 2);
 	if (!GetCrossPoint(strips[1], strips[2], points[0]) ||
-		!GetCrossPoint(strips[2], strips[3], points[1]) ||
-		!GetCrossPoint(strips[3], strips[1], points[2]) ||
+		!GetCrossPoint(strips[2], strips[0], points[1]) ||
+		!GetCrossPoint(strips[0], strips[1], points[2]) ||
 		(points[0] - points[1]).Mod2() > rad2 ||
 		(points[1] - points[2]).Mod2() > rad2 ||
 		(points[2] - points[0]).Mod2() > rad2) return false;

DataFormats/src/SigClusterTPC.cpp

@@ -1,5 +1,16 @@
 #include "SigClusterTPC.h"
 #include "EventTPC.h"
+
+class { //object keeping temporary envelope containers outside of SigClusterTPC
+	friend class SigClusterTPC;
+private:
+	std::set<std::tuple<direction, int, int>> hitList; // temporary list of selected envelope space-time cells for further analysis where return
+								 // value=key(STRIP_DIR [0-2], STRIP_NUM [1-1024], REBIN_TIME_CELL [0-511])
+	std::set< std::tuple<int, direction, int>> hitListByTimeDir; // same temporary list of selected envelope space-time cells as a map
+														  // with key=(REBIN_TIME_CELL [0-511], STRIP_DIR [0-2]) 
+														  // that returns vector of STRIP_NUMBERS 
+} envelope;
+
 /* ============= SPACE-TIME CLUSTER CLASS ===========*/
 
 SigClusterTPC::SigClusterTPC(EventTPC& e)
@@ -24,8 +35,8 @@ bool SigClusterTPC::AddEnvelopeHit(std::tuple<direction, int, int> hit) {  // va
 	if (time_cell < 0 || time_cell >= Geometry().GetAgetNtimecells() ||
 		strip_number < 1 || strip_number > Geometry().GetDirNstrips(strip_dir)) return false;
 
-	envelope_hitList.insert(hit);
-	envelope_hitListByTimeDir.insert({ time_cell,strip_dir,strip_number });
+	envelope.hitList.insert(hit);
+	envelope.hitListByTimeDir.insert({ time_cell,strip_dir,strip_number });
 	return true;
 }
 
@@ -37,10 +48,10 @@ void SigClusterTPC::UpdateStats() {
 }
 
 void SigClusterTPC::Combine() {
-	hitList.insert(envelope_hitList.begin(), envelope_hitList.end()); // insert envelope hit list
-	hitListByTimeDir.insert(envelope_hitListByTimeDir.begin(), envelope_hitListByTimeDir.end());
-	envelope_hitList.clear(); //clear envelope hit list
-	envelope_hitListByTimeDir.clear();
+	hitList.insert(envelope.hitList.begin(), envelope.hitList.end()); // insert envelope hit list
+	hitListByTimeDir.insert(envelope.hitListByTimeDir.begin(), envelope.hitListByTimeDir.end());
+	envelope.hitList.clear(); //clear envelope hit list
+	envelope.hitListByTimeDir.clear();
 }
 
 size_t SigClusterTPC::GetNhits(direction strip_dir) const {   // # of hits in a given direction
@@ -88,18 +99,17 @@ std::shared_ptr<TH2D> SigClusterTPC::GetStripVsTimeInMM(direction strip_dir) {
 // get three directions on: XY, XZ, YZ planes
 Reconstr_hist&& SigClusterTPC::Get(double radius, int rebin_space, int rebin_time, int method) {
 
-	static std::function<std::pair<bool, TVector2>(int, int, int, double)> fn = [](int s1_num, int s2_num, int s3_num, double rad)->std::pair<bool, TVector2> {
+	static std::function<std::pair<bool, TVector2>(std::array<int, 3>, double)> fn = [](std::array<int, 3> nums, double rad)->std::pair<bool, TVector2> {
 		TVector2 pos;
-		return { Geometry().MatchCrossPoint({s1_num, s2_num, s3_num}, rad, pos), pos };
+		return { Geometry().MatchCrossPoint(nums, rad, pos), pos };
 	};
-	static RV_Storage<std::pair<bool, TVector2>, int, int, int, double> MatchCrossPoint_functor(fn);
+	static RV_Storage<std::pair<bool, TVector2>, std::array<int, 3>, double> MatchCrossPoint_functor(fn); //return value manager, saves already found crossing strips
 
 	Reconstr_hist h_all;
 	bool err_flag = false;
 	if (std::any_of(dir_vec_UVW.begin(), dir_vec_UVW.end(), [&](direction dir_) { return GetNhits(dir_) < 1; })) return std::move(h_all);
 
 	//std::cout << Form(">>>> EventId = %lld", event_id) << std::endl;
-	//std::cout << Form(">>>> Time cell range = [%d, %d]", time_cell_min, time_cell_max) << std::endl;
 
 	for (auto hits_in_time_cell_begin = hitListByTimeDir.begin(); hits_in_time_cell_begin != hitListByTimeDir.end(); hits_in_time_cell_begin = hitListByTimeDir.upper_bound({ std::get<0>(*hits_in_time_cell_begin), std::numeric_limits<direction>::max(), std::numeric_limits<int>::max() })) { //iterate over time cells
 		auto icell = std::get<0>(*hits_in_time_cell_begin);
@@ -112,28 +122,21 @@ Reconstr_hist&& SigClusterTPC::Get(double radius, int rebin_space, int rebin_tim
 			}
 		}
 		// check if there is at least one hit in each direction
-		if (std::any_of(dir_vec_UVW.begin(), dir_vec_UVW.end(), [&](direction dir_) { return hits_strip_nums_in_dir[dir_].size() == 0; })) continue;
+		if (std::any_of(dir_vec_UVW.begin(), dir_vec_UVW.end(), [&](direction dir) { return hits_strip_nums_in_dir[dir].size() == 0; })) continue;
 
 		//   std::cout << Form(">>>> Number of hits: time cell=%d: U=%d / V=%d / W=%d",
 		//		      icell, (int)hits[int(direction::U)].size(), (int)hits[int(direction::V)].size(), (int)hits[int(direction::W)].size()) << std::endl;
 
-		std::map<std::pair<direction, int>, int> n_match; // map of number of matched points for each strip, key=STRIP_NUM [1-1024]
 		std::map<std::array<int, 3>, TVector2> hitPos; // key=(STRIP_NUM_U, STRIP_NUM_V, STRIP_NUM_W), value=(X [mm],Y [mm])
+		std::array<int, 3> pos;
 		// loop over hits and confirm matching in space
 		for (auto& hit_strip_num_U : hits_strip_nums_in_dir[direction::U]) {
 			for (auto& hit_strip_num_V : hits_strip_nums_in_dir[direction::V]) {
 				for (auto& hit_strip_num_W : hits_strip_nums_in_dir[direction::W]) {
-					TVector2 pos;
-					auto result = MatchCrossPoint_functor(hit_strip_num_U, hit_strip_num_V, hit_strip_num_W, radius);
+					pos = { hit_strip_num_U, hit_strip_num_V, hit_strip_num_W };
+					auto result = MatchCrossPoint_functor(pos, radius);
 					if (std::get<0>(result)) {
-						n_match[{direction::U, hit_strip_num_U}]++;
-						n_match[{direction::V, hit_strip_num_V}]++;
-						n_match[{direction::W, hit_strip_num_W}]++;
-						hitPos[{hit_strip_num_U, hit_strip_num_V, hit_strip_num_W}] = std::get<1>(result);
-						//	    std::cout << Form(">>>> Checking triplet: result = TRUE" ) << std::endl;
-					}
-					else {
-						//	    std::cout << Form(">>>> Checking triplet: result = FALSE" ) << std::endl;
+						hitPos[pos] = std::get<1>(result);
 					}
 				}
 			}
@@ -218,26 +221,12 @@ Reconstr_hist&& SigClusterTPC::Get(double radius, int rebin_space, int rebin_tim
 		// loop over matched hits and fill histograms
 		for (auto&& it : hitPos) {
 
-			double val = 0.0;
-
-			switch (method) {
-
-			case 0: // mehtod #1 - divide charge equally
-				val = std::inner_product(dir_vec_UVW.begin(), dir_vec_UVW.end(), it.first.begin(), 0.0, std::plus<>(), [&](direction dir, int key) {
-					return evt_ref.GetValByStrip(dir, key, icell) / n_match.at({ dir, key });
-				});
-				break;
-
-			case 1: // method #2 - divide charge according to charge-fraction in two other directions
-				val = 0.5 *
-					std::inner_product(dir_vec_UVW.begin(), dir_vec_UVW.end(), it.first.begin(), 0.0, std::plus<>(), [&](direction dir, int key) {
-					return evt_ref.GetValByStrip(dir, key, icell) * std::inner_product(fraction.begin(), fraction.end(), dir_vec_UVW.begin(), 0.0, std::plus<>(), [&](auto& map_, direction dir2) {
-						return (dir2 != dir ? map_.second.at(it.first) : 0.0);
-					});
+			double val = 0.5 * // divide charge according to charge-fraction in two other directions
+				std::inner_product(dir_vec_UVW.begin(), dir_vec_UVW.end(), it.first.begin(), 0.0, std::plus<>(), [&](direction dir, int key) {
+				return evt_ref.GetValByStrip(dir, key, icell) * std::inner_product(fraction.begin(), fraction.end(), dir_vec_UVW.begin(), 0.0, std::plus<>(), [&](auto& map_, direction dir2) {
+					return (dir2 != dir ? map_.second.at(it.first) : 0.0);
 				});
-				break;
-
-			}; // end of switch (method)...
+			});
 			double z = Geometry().Timecell2pos(icell);
 			h_all.second->Fill((it.second).X(), (it.second).Y(), z, val);
 			h_all.first->Fill((it.second).X(), (it.second).Y(), val);