Modified DeadChannels such that there is now a Sabre and Anasen version. Implemented in AnasenEfficiency. Small changes to StripDetector.

etc/AnasenDeadChannels.txt

@@ -0,0 +1,290 @@
+6 FQQQ 1 RING NONE 6
+16 FQQQ 0 RING NONE 0
+17 FQQQ 0 RING NONE 1
+18 FQQQ 0 RING NONE 2
+19 FQQQ 0 RING NONE 3
+20 FQQQ 0 RING NONE 4
+23 FQQQ 0 RING NONE 7
+24 FQQQ 0 RING NONE 8
+27 FQQQ 0 RING NONE 11
+28 FQQQ 0 RING NONE 12
+29 FQQQ 0 RING NONE 13
+30 FQQQ 0 RING NONE 14
+32 FQQQ 1 WEDGE NONE 15
+37 FQQQ 1 WEDGE NONE 10
+38 FQQQ 1 WEDGE NONE 9
+40 FQQQ 1 WEDGE NONE 7
+48 FQQQ 0 WEDGE NONE 15
+49 FQQQ 0 WEDGE NONE 14
+50 FQQQ 0 WEDGE NONE 13
+51 FQQQ 0 WEDGE NONE 12
+52 FQQQ 0 WEDGE NONE 11
+59 FQQQ 0 WEDGE NONE 4
+64 FQQQ 2 RING NONE 0
+65 FQQQ 2 RING NONE 1
+66 FQQQ 2 RING NONE 2
+67 FQQQ 2 RING NONE 3
+68 FQQQ 2 RING NONE 4
+69 FQQQ 2 RING NONE 5
+70 FQQQ 2 RING NONE 6
+71 FQQQ 2 RING NONE 7
+72 FQQQ 2 RING NONE 8
+73 FQQQ 2 RING NONE 9
+74 FQQQ 2 RING NONE 10
+75 FQQQ 2 RING NONE 11
+76 FQQQ 2 RING NONE 12
+77 FQQQ 2 RING NONE 13
+78 FQQQ 2 RING NONE 14
+79 FQQQ 2 RING NONE 15
+80 FQQQ 3 RING NONE 0
+81 FQQQ 3 RING NONE 1
+82 FQQQ 3 RING NONE 2
+83 FQQQ 3 RING NONE 3
+84 FQQQ 3 RING NONE 4
+85 FQQQ 3 RING NONE 5
+86 FQQQ 3 RING NONE 6
+87 FQQQ 3 RING NONE 7
+88 FQQQ 3 RING NONE 8
+89 FQQQ 3 RING NONE 9
+90 FQQQ 3 RING NONE 10
+91 FQQQ 3 RING NONE 11
+92 FQQQ 3 RING NONE 12
+93 FQQQ 3 RING NONE 13
+94 FQQQ 3 RING NONE 14
+95 FQQQ 3 RING NONE 15
+96 FQQQ 2 WEDGE NONE 15
+97 FQQQ 2 WEDGE NONE 14
+98 FQQQ 2 WEDGE NONE 13
+99 FQQQ 2 WEDGE NONE 12
+100 FQQQ 2 WEDGE NONE 11
+101 FQQQ 2 WEDGE NONE 10
+102 FQQQ 2 WEDGE NONE 9
+103 FQQQ 2 WEDGE NONE 8
+104 FQQQ 2 WEDGE NONE 7
+105 FQQQ 2 WEDGE NONE 6
+106 FQQQ 2 WEDGE NONE 5
+107 FQQQ 2 WEDGE NONE 4
+108 FQQQ 2 WEDGE NONE 3
+109 FQQQ 2 WEDGE NONE 2
+110 FQQQ 2 WEDGE NONE 1
+111 FQQQ 2 WEDGE NONE 0
+114 FQQQ 3 WEDGE NONE 13
+127 FQQQ 3 WEDGE NONE 0
+128 BARREL1B B FRONT DOWN 7
+129 BARREL1B B FRONT UP 6
+130 BARREL1B B FRONT DOWN 5
+131 BARREL1B B FRONT UP 4
+132 BARREL1B B FRONT UP 3
+133 BARREL1B B FRONT DOWN 2
+134 BARREL1B B FRONT UP 1
+135 BARREL1B B FRONT DOWN 0
+136 BARREL1B A FRONT DOWN 7
+137 BARREL1B A FRONT UP 6
+138 BARREL1B A FRONT DOWN 5
+139 BARREL1B A FRONT UP 4
+140 BARREL1B A FRONT UP 3
+141 BARREL1B A FRONT DOWN 2
+142 BARREL1B A FRONT UP 1
+143 BARREL1B A FRONT DOWN 0
+144 BARREL1B C FRONT DOWN 7
+145 BARREL1B C FRONT UP 6
+168 BARREL1B E FRONT DOWN 7
+169 BARREL1B E FRONT UP 6
+170 BARREL1B E FRONT DOWN 5
+171 BARREL1B E FRONT UP 4
+172 BARREL1B E FRONT UP 3
+173 BARREL1B E FRONT DOWN 2
+174 BARREL1B E FRONT UP 1
+175 BARREL1B E FRONT DOWN 0
+176 BARREL1A B FRONT DOWN 7
+177 BARREL1A B FRONT UP 6
+178 BARREL1A B FRONT DOWN 5
+179 BARREL1A B FRONT UP 4
+180 BARREL1A B FRONT UP 3
+181 BARREL1A B FRONT DOWN 2
+182 BARREL1A B FRONT UP 1
+183 BARREL1A B FRONT DOWN 0
+184 BARREL1A A FRONT DOWN 7
+185 BARREL1A A FRONT UP 6
+186 BARREL1A A FRONT DOWN 5
+187 BARREL1A A FRONT UP 4
+188 BARREL1A A FRONT UP 3
+189 BARREL1A A FRONT DOWN 2
+190 BARREL1A A FRONT UP 1
+191 BARREL1A A FRONT DOWN 0
+192 BQQQ 3 RING NONE 0
+193 BQQQ 3 RING NONE 1
+194 BQQQ 3 RING NONE 2
+195 BQQQ 3 RING NONE 3
+196 BQQQ 3 RING NONE 4
+197 BQQQ 3 RING NONE 5
+198 BQQQ 3 RING NONE 6
+199 BQQQ 3 RING NONE 7
+200 BQQQ 3 RING NONE 8
+201 BQQQ 3 RING NONE 9
+202 BQQQ 3 RING NONE 10
+203 BQQQ 3 RING NONE 11
+204 BQQQ 3 RING NONE 12
+205 BQQQ 3 RING NONE 13
+206 BQQQ 3 RING NONE 14
+207 BQQQ 3 RING NONE 15
+208 BQQQ 2 RING NONE 0
+209 BQQQ 2 RING NONE 1
+210 BQQQ 2 RING NONE 2
+211 BQQQ 2 RING NONE 3
+212 BQQQ 2 RING NONE 4
+213 BQQQ 2 RING NONE 5
+214 BQQQ 2 RING NONE 6
+215 BQQQ 2 RING NONE 7
+216 BQQQ 2 RING NONE 8
+217 BQQQ 2 RING NONE 9
+218 BQQQ 2 RING NONE 10
+219 BQQQ 2 RING NONE 11
+220 BQQQ 2 RING NONE 12
+221 BQQQ 2 RING NONE 13
+222 BQQQ 2 RING NONE 14
+223 BQQQ 2 RING NONE 15
+224 BARREL2A B FRONT DOWN 7
+225 BARREL2A B FRONT UP 6
+226 BARREL2A B FRONT DOWN 5
+227 BARREL2A B FRONT UP 4
+228 BARREL2A B FRONT UP 3
+229 BARREL2A B FRONT DOWN 2
+230 BARREL2A B FRONT UP 1
+231 BARREL2A B FRONT DOWN 0
+232 BARREL2A A FRONT DOWN 7
+233 BARREL2A A FRONT UP 6
+234 BARREL2A A FRONT DOWN 5
+235 BARREL2A A FRONT UP 4
+236 BARREL2A A FRONT UP 3
+237 BARREL2A A FRONT DOWN 2
+238 BARREL2A A FRONT UP 1
+239 BARREL2A A FRONT DOWN 0
+240 BARREL2A F FRONT DOWN 7
+241 BARREL2A F FRONT UP 6
+242 BARREL2A F FRONT DOWN 5
+243 BARREL2A F FRONT UP 4
+244 BARREL2A F FRONT UP 3
+245 BARREL2A F FRONT DOWN 2
+246 BARREL2A F FRONT UP 1
+247 BARREL2A F FRONT DOWN 0
+248 BARREL2A E FRONT DOWN 7
+249 BARREL2A E FRONT UP 6
+250 BARREL2A E FRONT DOWN 5
+251 BARREL2A E FRONT UP 4
+252 BARREL2A E FRONT UP 3
+253 BARREL2A E FRONT DOWN 2
+254 BARREL2A E FRONT UP 1
+255 BARREL2A E FRONT DOWN 0
+265 BQQQ 1 RING NONE 9
+266 BQQQ 1 RING NONE 10
+304 BARREL2B B FRONT DOWN 7
+305 BARREL2B B FRONT UP 6
+306 BARREL2B B FRONT DOWN 5
+307 BARREL2B B FRONT UP 4
+308 BARREL2B B FRONT UP 3
+309 BARREL2B B FRONT DOWN 2
+310 BARREL2B B FRONT UP 1
+311 BARREL2B B FRONT DOWN 0
+312 BARREL2B A FRONT DOWN 7
+313 BARREL2B A FRONT UP 6
+314 BARREL2B A FRONT DOWN 5
+315 BARREL2B A FRONT UP 4
+316 BARREL2B A FRONT UP 3
+317 BARREL2B A FRONT DOWN 2
+318 BARREL2B A FRONT UP 1
+319 BARREL2B A FRONT DOWN 0
+334 BARREL1B F BACK NONE 1
+335 BARREL1B F BACK NONE 0
+343 BARREL1A D BACK NONE 0
+347 BARREL1A E BACK NONE 0
+358 BARREL2B F FRONT UP 1
+359 BARREL2B F FRONT DOWN 0
+368 BARREL2A C FRONT DOWN 7
+369 BARREL2A C FRONT UP 6
+370 BARREL2A C FRONT DOWN 5
+371 BARREL2A C FRONT UP 4
+372 BARREL2A C FRONT UP 3
+373 BARREL2A C FRONT DOWN 2
+374 BARREL2A C FRONT UP 1
+375 BARREL2A C FRONT DOWN 0
+376 BARREL2A D FRONT DOWN 7
+377 BARREL2A D FRONT UP 6
+378 BARREL2A D FRONT DOWN 5
+379 BARREL2A D FRONT UP 4
+380 BARREL2A D FRONT UP 3
+381 BARREL2A D FRONT DOWN 2
+382 BARREL2A D FRONT UP 1
+383 BARREL2A D FRONT DOWN 0
+386 BARREL2B C BACK NONE 1
+394 BARREL2B E BACK NONE 1
+398 BARREL2B F BACK NONE 1
+402 BARREL2A C BACK NONE 1
+406 BARREL2A D BACK NONE 1
+408 BARREL2A E BACK NONE 3
+409 BARREL2A E BACK NONE 2
+410 BARREL2A E BACK NONE 1
+411 BARREL2A E BACK NONE 0
+412 BARREL2A F BACK NONE 3
+413 BARREL2A F BACK NONE 2
+414 BARREL2A F BACK NONE 1
+415 BARREL2A F BACK NONE 0
+416 BARREL1A F FRONT DOWN 7
+417 BARREL1A F FRONT UP 6
+418 BARREL1A F FRONT DOWN 5
+419 BARREL1A F FRONT UP 4
+420 BARREL1A F FRONT UP 3
+421 BARREL1A F FRONT DOWN 2
+422 BARREL1A F FRONT UP 1
+423 BARREL1A F FRONT DOWN 0
+424 BARREL1A E FRONT DOWN 7
+425 BARREL1A E FRONT UP 6
+426 BARREL1A E FRONT DOWN 5
+427 BARREL1A E FRONT UP 4
+428 BARREL1A E FRONT UP 3
+429 BARREL1A E FRONT DOWN 2
+430 BARREL1A E FRONT UP 1
+431 BARREL1A E FRONT DOWN 0
+432 BARREL1A C FRONT DOWN 7
+433 BARREL1A C FRONT UP 6
+434 BARREL1A C FRONT DOWN 5
+435 BARREL1A C FRONT UP 4
+436 BARREL1A C FRONT UP 3
+437 BARREL1A C FRONT DOWN 2
+438 BARREL1A C FRONT UP 1
+439 BARREL1A C FRONT DOWN 0
+440 BARREL1A D FRONT DOWN 7
+441 BARREL1A D FRONT UP 6
+442 BARREL1A D FRONT DOWN 5
+443 BARREL1A D FRONT UP 4
+444 BARREL1A D FRONT UP 3
+445 BARREL1A D FRONT DOWN 2
+446 BARREL1A D FRONT UP 1
+447 BARREL1A D FRONT DOWN 0
+448 BARREL2A A BACK NONE 3
+449 BARREL2A A BACK NONE 2
+450 BARREL2A A BACK NONE 1
+451 BARREL2A A BACK NONE 0
+454 BARREL2A B BACK NONE 1
+457 BARREL1A A BACK NONE 2
+461 BARREL1A B BACK NONE 2
+462 BARREL1A B BACK NONE 1
+463 BARREL1A B BACK NONE 0
+464 BARREL2B A BACK NONE 3
+465 BARREL2B A BACK NONE 2
+466 BARREL2B A BACK NONE 1
+467 BARREL2B A BACK NONE 0
+468 BARREL2B B BACK NONE 3
+469 BARREL2B B BACK NONE 2
+470 BARREL2B B BACK NONE 1
+471 BARREL2B B BACK NONE 0
+472 BARREL1B A BACK NONE 3
+473 BARREL1B A BACK NONE 2
+474 BARREL1B A BACK NONE 1
+475 BARREL1B A BACK NONE 0
+476 BARREL1B B BACK NONE 3
+477 BARREL1B B BACK NONE 2
+478 BARREL1B B BACK NONE 1
+479 BARREL1B B BACK NONE 0
+528 BQQQ 3 WEDGE NONE 15
+537 BQQQ 3 WEDGE NONE 6

etc/convertedAnasenCoords.txt

@@ -1,30 +1,30 @@
 For SX3's: det id | mid rho | mid z | mid phi | name
 For QQQ's: det id | mid phi | z | name
-0 5.49779 0 QQQ_0
+0 2.30659 0 QQQ_0
-1 0.785398 0.8028 QQQ_1
+1 7.01897 0.8028 QQQ_1
-2 2.35619 0 QQQ_2
+2 5.44818 0 QQQ_2
-3 3.92699 0 QQQ_3
+3 3.87738 0 QQQ_3
-6 8.90601 6.25 0.785795 R1a_A
+6 8.90601 6.25 5.49739 R1a_A
-7 8.89871 6.25 0.262014 R1a_B
+7 8.89871 6.25 6.02117 R1a_B
-8 8.90354 6.25 6.02132 R1a_C
+8 8.90354 6.25 0.261868 R1a_C
-9 8.90247 6.25 5.49779 R1a_D
+9 8.90247 6.25 0.785398 R1a_D
-10 8.90354 6.25 4.97426 R1a_E
+10 8.90354 6.25 1.30893 R1a_E
-11 8.90354 6.25 4.45052 R1a_F
+11 8.90354 6.25 1.83266 R1a_F
-12 8.90247 6.25 3.92699 R1b_A
+12 8.90247 6.25 2.35619 R1b_A
-13 8.90354 6.25 3.40346 R1b_B
+13 8.90354 6.25 2.87972 R1b_B
-14 8.90354 6.25 2.87972 R1b_C
+14 8.90354 6.25 3.40346 R1b_C
-15 8.90247 6.25 2.35619 R1b_D
+15 8.90247 6.25 3.92699 R1b_D
-16 8.90354 6.25 1.83266 R1b_E
+16 8.90354 6.25 4.45052 R1b_E
-17 8.90354 6.25 1.30893 R1b_F
+17 8.90354 6.25 4.97426 R1b_F
-18 8.90601 18.65 0.785795 R2a_A
+18 8.90601 18.65 5.49739 R2a_A
-19 8.89871 18.65 0.262014 R2a_B
+19 8.89871 18.65 6.02117 R2a_B
-20 8.90354 18.65 6.02132 R2a_C
+20 8.90354 18.65 0.261868 R2a_C
-21 8.90247 18.65 5.49779 R2a_D
+21 8.90247 18.65 0.785398 R2a_D
-22 8.90354 18.65 4.97426 R2a_E
+22 8.90354 18.65 1.30893 R2a_E
-23 8.90354 18.65 4.45052 R2a_F
+23 8.90354 18.65 1.83266 R2a_F
-24 8.90247 18.65 3.92699 R2b_A
+24 8.90247 18.65 2.35619 R2b_A
-25 8.90354 18.65 3.40346 R2b_B
+25 8.90354 18.65 2.87972 R2b_B
-26 8.90354 18.65 2.87972 R2b_C
+26 8.90354 18.65 3.40346 R2b_C
-27 8.90247 18.65 2.35619 R2b_D
+27 8.90247 18.65 3.92699 R2b_D
-28 8.90354 18.65 1.83266 R2b_E
+28 8.90354 18.65 4.45052 R2b_E
-29 8.90354 18.65 1.30893 R2b_F
+29 8.90354 18.65 4.97426 R2b_F

include/AnasenDeadChannelMap.h

@@ -0,0 +1,53 @@
+#ifndef ANASENDEADCHANNELMAP_H
+#define ANASENDEADCHANNELMAP_H
+
+#include <unordered_map>
+#include <string>
+
+enum class AnasenDetectorType
+{
+	Barrel1,
+	Barrel2,
+	FQQQ,
+	BQQQ
+};
+
+enum AnasenDetectorSide
+{
+	Front=0,
+	Back=1
+};
+
+class AnasenDeadChannelMap
+{
+public:
+	AnasenDeadChannelMap();
+	AnasenDeadChannelMap(const std::string& filename);
+	~AnasenDeadChannelMap();
+	void LoadMapfile(const std::string& filename);
+	inline const bool IsValid() const { return valid_flag; }
+	const bool IsDead(AnasenDetectorType type, int detIndex, int channel, AnasenDetectorSide side) const;
+private:
+	void InitMap();
+	int ConvertStringTypeIndexToOffset(const std::string& type, const std::string& index, const std::string& side);
+	bool valid_flag;
+	std::unordered_map<int, bool> dcMap;
+
+	/*
+		global channel calculated as 
+		detTypeOffset + detIndex*total + detSideOffset + channel
+	*/
+
+	const int nfronts_sx3 = 4;
+	const int nfronts_qqq = 16;
+	const int nchannels_sx3 = 8;
+	const int nchannels_qqq = 32;
+
+	const int barrel1_offset = 0;
+	const int barrel2_offset = 12*nchannels_sx3;
+	const int fqqq_offset = barrel2_offset + 12*nchannels_sx3;
+	const int bqqq_offset = fqqq_offset + 4*nchannels_qqq;
+
+};
+
+#endif

include/AnasenEfficiency.h

@@ -9,6 +9,7 @@
 #include "Target.h"
 #include "Nucleus.h"
 #include "MaskFile.h"
+#include "AnasenDeadChannelMap.h"
 
 struct DetectorResult {
 	bool detectFlag = false;
@@ -24,6 +25,7 @@ public:
 	void CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) override;
 	void DrawDetectorSystem(const std::string& filename) override;
 	double RunConsistencyCheck() override;
+	inline void SetDeadChannelMap(const std::string& filename) { dmap.LoadMapfile(filename); }
 
 private:
 	DetectorResult IsRing1(Mask::Nucleus& nucleus);
@@ -38,6 +40,8 @@ private:
 
 	Mask::Target det_silicon;
 
+	AnasenDeadChannelMap dmap;
+
 	/**** ANASEN geometry constants *****/
 	const int n_sx3_per_ring = 12;
 	const int n_qqq = 4;
@@ -54,7 +58,7 @@ private:
 	const double qqq_z[4] = {qqq_nom_z, qqq_nom_z, qqq_nom_z, qqq_nom_z};
 	const double qqq_phi[4] = {5.49779, 0.785398, 2.35619, 3.92699};
 	const double ring_rho[12] = {0.0890601, 0.0889871, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354, 0.0890247, 0.0890354, 0.0890354};
-	const double ring_phi[12] = {0.785795, 0.262014, 6.02132, 5.49779, 4.97426, 4.45052, 3.92699, 3.40346, 2.87972, 2.35619, 1.83266, 1.30893};
+	const double ring_phi[12] = {4.97426, 5.49739, 6.02132, 0.261868, 0.785398, 1.30893, 1.83266, 2.35619, 2.87972, 3.40346, 3.92699, 4.45052};
 	/*************************/
 
 	static constexpr double threshold = 0.6; //MeV

include/QQQDetector.h

@@ -3,7 +3,7 @@
 	Class implementing geometry for QQQDetector where the detector is perpendicular to the beam axis.
 	Detector is first generated centered on the x-axis (phi=0)
 
-	Coordinate convention : +z is downstream, -z is upstream. +y is vertically up in the lab.
+	Coordinate convention : +z is downstream, -z is upstream. +y is vertically down in the lab.
 */
 #ifndef QQQDETECTOR_H
 #define QQQDETECTOR_H

include/SabreDeadChannelMap.h

@@ -1,16 +1,17 @@
-#ifndef DEADCHANNELMAP_H
+#ifndef SABREDEADCHANNELMAP_H
-#define DEADCHANNELMAP_H
+#define SABREDEADCHANNELMAP_H
 
 #include <string>
 #include <iostream>
 #include <unordered_map>
 
-class DeadChannelMap {
+class SabreDeadChannelMap 
+{
 public:
-	DeadChannelMap();
+	SabreDeadChannelMap();
-	DeadChannelMap(std::string& name);
+	SabreDeadChannelMap(const std::string& name);
-	~DeadChannelMap();
+	~SabreDeadChannelMap();
-	void LoadMapfile(std::string& name);
+	void LoadMapfile(const std::string& name);
 	bool IsDead(int det, int channel, int ringwedgeFlag);
 	bool IsValid() { return validFlag; };
 

include/SabreEfficiency.h

@@ -4,7 +4,7 @@
 #include "DetectorEfficiency.h"
 #include "SabreDetector.h"
 #include "Target.h"
-#include "DeadChannelMap.h"
+#include "SabreDeadChannelMap.h"
 #include "Nucleus.h"
 
 class SabreEfficiency : public DetectorEfficiency {
@@ -23,7 +23,7 @@ private:
     
 	Mask::Target deadlayer;
     Mask::Target sabre_eloss;
-    DeadChannelMap dmap;
+    SabreDeadChannelMap dmap;
 
 	//Sabre constants
 	const double INNER_R = 0.0326;

include/StripDetector.h

@@ -2,12 +2,12 @@
 #define __STRIP_DETECTOR_H
 
 // +z is along beam axis
-// +y is vertically "upward" in the lab frame
+// +y is vertically "downward" in the lab frame
 
 //angles must be in radians, but distances can be whatever
 //PROVIDED all input distances are the same
 
-//Front strips from lowest y to highest y
+//Front strips from largest y to smallest y
 
 //Back strips from lowest z to highest z
 
@@ -19,6 +19,13 @@
 #include "Rotation.h"
 #include "RandomGenerator.h"
 
+struct StripHit
+{
+	int front_strip_index=-1;
+	int back_strip_index=-1;
+	double front_ratio=0.0;
+};
+
 class StripDetector {
 public:
   
@@ -34,7 +41,7 @@ public:
 	inline void TurnOffRandomizedCoordinates() { rndmFlag = false; }
 
 	Mask::Vec3 GetHitCoordinates(int front_stripch, double front_strip_ratio);
-	std::pair<int,double> GetChannelRatio(double theta, double phi);
+	StripHit GetChannelRatio(double theta, double phi);
 
 private:
 	inline bool ValidChannel(int f) { return ((f >= 0 && f < num_strips) ? true : false); };

input.txt

@@ -1,12 +1,12 @@
 ----------Data Information----------
-OutputFile: /media/gordon/ANASENData/MaskData/Kinematics/7Bedp_870keV_beam_50CD2_test.mask
+OutputFile: /data1/gwm17/7BeNov2021/sim/7Bedp_1100keV_beam_50CD2.mask
 ----------Reaction Information----------
 ReactionType: TwoStepRxn
 Z	A (order is target, projectile, ejectile, break1, break3(if pure decay is target, ejectile))
-1   2
+1	2
 4	7
 1	1
-2   4
+2	4
 ----------Target Information----------
 NumberOfLayers: 1
 begin_layer
@@ -18,7 +18,7 @@ begin_layer
 end_layer
 ----------Sampling Information----------
 NumberOfSamples: 100000
-BeamMeanEnergy(MeV): 0.87 BeamEnergySigma(MeV): 0.0
+BeamMeanEnergy(MeV): 1.1 BeamEnergySigma(MeV): 0.0
 EjectileThetaType(0=Lab,1=CM): 1
 EjectileThetaMin(deg): 0.0 EjectileThetaMax(deg): 180.0
 EjectilePhiMin(deg): 0.0 EjectilePhiMax(deg): 360.0

src/Detectors/AnasenDeadChannelMap.cpp

@@ -0,0 +1,133 @@
+#include "AnasenDeadChannelMap.h"
+#include <fstream>
+#include <iostream>
+#include <cmath>
+
+AnasenDeadChannelMap::AnasenDeadChannelMap() :
+	valid_flag(false)
+{
+	InitMap();
+}
+
+AnasenDeadChannelMap::AnasenDeadChannelMap(const std::string& filename) :
+	valid_flag(false)
+{
+	InitMap();
+	LoadMapfile(filename);
+}
+
+AnasenDeadChannelMap::~AnasenDeadChannelMap() {}
+
+void AnasenDeadChannelMap::InitMap()
+{
+	for(int i=0; i<(bqqq_offset+4.0*nchannels_qqq); i++)
+		dcMap[i] = false;
+}
+
+int AnasenDeadChannelMap::ConvertStringTypeIndexToOffset(const std::string& type, const std::string& index, const std::string& side)
+{
+	int channel_offset=0;
+	if(type == "BARREL1B")
+		channel_offset += 6*nchannels_sx3;
+	else if (type == "BARREL2A")
+		channel_offset += barrel2_offset;
+	else if (type == "BARREL2B")
+		channel_offset += barrel2_offset + 6*nchannels_sx3;
+	else if (type == "FQQQ")
+		channel_offset += fqqq_offset;
+	else if (type == "BQQQ")
+		channel_offset += bqqq_offset;
+
+	if(index == "B")
+		channel_offset += nchannels_sx3;
+	else if(index == "C")
+		channel_offset += 2.0*nchannels_sx3;
+	else if(index == "D")
+		channel_offset += 3.0*nchannels_sx3;
+	else if(index == "E")
+		channel_offset += 4.0*nchannels_sx3;
+	else if(index == "F")
+		channel_offset += 5.0*nchannels_sx3;
+	else if(index == "1")
+		channel_offset += nchannels_qqq;
+	else if(index == "2")
+		channel_offset += 2.0*nchannels_qqq;
+	else if(index == "3")
+		channel_offset += 3.0*nchannels_qqq;
+
+	if(side == "BACK")
+		channel_offset += nfronts_sx3;
+	else if(side == "WEDGE")
+		channel_offset += nfronts_qqq;
+
+	return channel_offset;
+}
+
+void AnasenDeadChannelMap::LoadMapfile(const std::string& filename)
+{
+	valid_flag = false;
+	std::ifstream input(filename);
+	if(!input.is_open())
+	{
+		std::cerr<<"Unable to open input file at AnasenDeadChannelMap::LoadMapfile(). Exiting."<<std::endl;
+		return;
+	}
+
+	std::string junk;
+	std::string type, side, index;
+	int channel;
+
+	int dead_channel;
+	while(input>>junk)
+	{
+		input>>type>>index>>side>>junk>>channel;
+
+		if(side == "FRONT")
+		{
+			channel = std::floor(channel/2.0);
+		}
+
+		dead_channel = ConvertStringTypeIndexToOffset(type, index, side) + channel;
+		
+		dcMap[dead_channel] = true;
+	}
+
+	valid_flag = true;
+}
+
+const bool AnasenDeadChannelMap::IsDead(AnasenDetectorType type, int detIndex, int channel, AnasenDetectorSide side) const
+{
+	int channel_index=-1;
+	switch(type)
+	{
+		case AnasenDetectorType::Barrel1:
+		{
+			channel_index = detIndex*nchannels_sx3 + side*nfronts_sx3 + channel;
+			break;
+		}
+		case AnasenDetectorType::Barrel2:
+		{
+			channel_index = barrel2_offset + detIndex*nchannels_sx3 + side*nfronts_sx3 + channel;
+			break;
+		}
+		case AnasenDetectorType::FQQQ:
+		{
+			channel_index = fqqq_offset + detIndex*nchannels_qqq + side*nfronts_qqq + channel;
+			break;
+		}
+		case AnasenDetectorType::BQQQ:
+		{
+			channel_index = bqqq_offset + detIndex*nchannels_qqq + side*nfronts_qqq + channel;
+			break;
+		}
+	}
+
+	auto data = dcMap.find(channel_index);
+	if(data == dcMap.end())
+	{
+		std::cerr<<"Bad channel index "<<channel_index<<" at AnasenDeadChannelMap::IsDead()"<<std::endl;
+		return false;
+	}
+
+	return data->second;
+}

src/Detectors/AnasenEfficiency.cpp

@@ -158,7 +158,7 @@ double AnasenEfficiency::RunConsistencyCheck() {
 	for(auto& point : r1_points) {
 		for(auto& sx3 : m_Ring1) {
 			auto result = sx3.GetChannelRatio(point.GetTheta(), point.GetPhi());
-			coords = sx3.GetHitCoordinates(result.first, result.second);
+			coords = sx3.GetHitCoordinates(result.front_strip_index, result.front_ratio);
 			if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
 				count++;
 				break;
@@ -168,7 +168,7 @@ double AnasenEfficiency::RunConsistencyCheck() {
 	for(auto& point : r2_points) {
 		for(auto& sx3 : m_Ring2) {
 			auto result = sx3.GetChannelRatio(point.GetTheta(), point.GetPhi());
-			coords = sx3.GetHitCoordinates(result.first, result.second);
+			coords = sx3.GetHitCoordinates(result.front_strip_index, result.front_ratio);
 			if(IsDoubleEqual(point.GetX(), coords.GetX()) && IsDoubleEqual(point.GetY(), coords.GetY()) && IsDoubleEqual(point.GetZ(), coords.GetZ())) {
 				count++;
 				break;
@@ -205,14 +205,15 @@ double AnasenEfficiency::RunConsistencyCheck() {
 DetectorResult AnasenEfficiency::IsRing1(Mask::Nucleus& nucleus) {
 
 	DetectorResult observation;
-	//Mask::Vec3 coords;
 	double thetaIncident;
-	for(auto& sx3 : m_Ring1) {
+	for(int i=0; i<n_sx3_per_ring; i++) {
-		auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
+		auto result = m_Ring1[i].GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
-		if(result.first != -1) {
+		if(result.front_strip_index != -1 /*&& !dmap.IsDead(AnasenDetectorType::Barrel1, i, result.front_strip_index, AnasenDetectorSide::Front)*/
+			&& !dmap.IsDead(AnasenDetectorType::Barrel1, i, result.back_strip_index, AnasenDetectorSide::Back)) 
+		{
 			observation.detectFlag = true;
-			observation.direction = sx3.GetHitCoordinates(result.first, result.second);
+			observation.direction = m_Ring1[i].GetHitCoordinates(result.front_strip_index, result.front_ratio);
-			thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR());
+			thetaIncident = std::acos(observation.direction.Dot(m_Ring1[i].GetNormRotated())/observation.direction.GetR());
 			if(thetaIncident > M_PI/2.0)
 				thetaIncident = M_PI - thetaIncident;
 
@@ -229,12 +230,14 @@ DetectorResult AnasenEfficiency::IsRing2(Mask::Nucleus& nucleus) {
 
 	DetectorResult observation;
 	double thetaIncident;
-	for(auto& sx3 : m_Ring2) {
+	for(int i=0; i<n_sx3_per_ring; i++) {
-		auto result = sx3.GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
+		auto result = m_Ring2[i].GetChannelRatio(nucleus.GetTheta(), nucleus.GetPhi());
-		if(result.first != -1) {
+		if(result.front_strip_index != -1 /*&& !dmap.IsDead(AnasenDetectorType::Barrel2, i, result.front_strip_index, AnasenDetectorSide::Front)*/
+			&& !dmap.IsDead(AnasenDetectorType::Barrel2, i, result.back_strip_index, AnasenDetectorSide::Back)) 
+		{
 			observation.detectFlag = true;
-			observation.direction = sx3.GetHitCoordinates(result.first, result.second);
+			observation.direction = m_Ring2[i].GetHitCoordinates(result.front_strip_index, result.front_ratio);
-			thetaIncident = std::acos(observation.direction.Dot(sx3.GetNormRotated())/observation.direction.GetR());
+			thetaIncident = std::acos(observation.direction.Dot(m_Ring2[i].GetNormRotated())/observation.direction.GetR());
 			if(thetaIncident > M_PI/2.0)
 				thetaIncident = M_PI - thetaIncident;
 
@@ -251,12 +254,14 @@ DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) {
 
 	DetectorResult observation;
 	double thetaIncident;
-	for(auto& qqq : m_forwardQQQs) {
+	for(int i=0; i<n_qqq; i++) {
-		auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
+		auto result = m_forwardQQQs[i].GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
-		if(result.first != -1) {
+		if(result.first != -1 /*&& !dmap.IsDead(AnasenDetectorType::FQQQ, i, result.first, AnasenDetectorSide::Front)*/ &&
+			!dmap.IsDead(AnasenDetectorType::FQQQ, i, result.second, AnasenDetectorSide::Back)) 
+		{
 			observation.detectFlag = true;
-			observation.direction = qqq.GetHitCoordinates(result.first, result.second);
+			observation.direction = m_forwardQQQs[i].GetHitCoordinates(result.first, result.second);
-			thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR());
+			thetaIncident = std::acos(observation.direction.Dot(m_forwardQQQs[i].GetNorm())/observation.direction.GetR());
 			if(thetaIncident > M_PI/2.0)
 				thetaIncident = M_PI - thetaIncident;
 
@@ -267,12 +272,14 @@ DetectorResult AnasenEfficiency::IsQQQ(Mask::Nucleus& nucleus) {
 	}
 
 	
-	for(auto& qqq : m_backwardQQQs) {
+	for(int i=0; i<n_qqq; i++) {
-		auto result = qqq.GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
+		auto result = m_backwardQQQs[i].GetTrajectoryRingWedge(nucleus.GetTheta(), nucleus.GetPhi());
-		if(result.first != -1) {
+		if(result.first != -1 /*&& !dmap.IsDead(AnasenDetectorType::BQQQ, i, result.first, AnasenDetectorSide::Front)*/ &&
+			!dmap.IsDead(AnasenDetectorType::BQQQ, i, result.second, AnasenDetectorSide::Back)) 
+		{
 			observation.detectFlag = true;
-			observation.direction = qqq.GetHitCoordinates(result.first, result.second);
+			observation.direction = m_backwardQQQs[i].GetHitCoordinates(result.first, result.second);
-			thetaIncident = std::acos(observation.direction.Dot(qqq.GetNorm())/observation.direction.GetR());
+			thetaIncident = std::acos(observation.direction.Dot(m_backwardQQQs[i].GetNorm())/observation.direction.GetR());
 			if(thetaIncident > M_PI/2.0)
 				thetaIncident = M_PI - thetaIncident;
 
@@ -381,6 +388,14 @@ void AnasenEfficiency::CountCoincidences(const Mask::MaskFileData& data, std::ve
 }
 
 void AnasenEfficiency::CalculateEfficiency(const std::string& inputname, const std::string& outputname, const std::string& statsname) {
+
+	if(!dmap.IsValid())
+	{
+		std::cerr<<"Invalid Dead Channel Map at AnasenEfficiency::CalculateEfficiency()! If you want to run with all possible";
+		std::cerr<<"channels active, simply load an empty file."<<std::endl;
+		return;
+	}
+
 	std::cout<<"----------ANASEN Efficiency Calculation----------"<<std::endl;
 	std::cout<<"Loading in output from kinematics simulation: "<<inputname<<std::endl;
 	std::cout<<"Running efficiency calculation..."<<std::endl;

src/Detectors/DeadChannelMap.cpp

@@ -1,68 +0,0 @@
-#include "DeadChannelMap.h"
-#include <fstream>
-
-DeadChannelMap::DeadChannelMap() :
-	validFlag(false)
-{
-	int maxchan = 5*24+1*16+7;
-	for(int i=0; i<maxchan; i++) {
-		dcMap[i] = false;
-	}
-}
-
-DeadChannelMap::DeadChannelMap(std::string& name) :
-	validFlag(false)
-{
-	LoadMapfile(name);
-}
-
-DeadChannelMap::~DeadChannelMap() {}
-
-void DeadChannelMap::LoadMapfile(std::string& name) {
-	std::ifstream input(name);
-	if(!input.is_open()) {
-		std::cerr<<"Unable to load dead channels, file: "<<name<<" could not be opened"<<std::endl;
-		validFlag = false;
-		return;
-	}
-
-	
-
-	std::string junk, rw;
-	int detID, channel;
-	int this_channel;
-
-	std::getline(input, junk);
-	while(input>>detID) {
-		input>>rw>>channel;
-		if(rw == "RING") {
-			this_channel = detID*24+RING*16+channel;
-			dcMap[this_channel] = true;
-		} else if(rw == "WEDGE") {
-			this_channel = detID*24+WEDGE*16+channel;
-			dcMap[this_channel] = true;
-		} else {
-			std::cerr<<"Invalid ring/wedge designation at DeadChannelMap"<<std::endl;
-		}
-	}
-
-	validFlag = true;
-}
-
-bool DeadChannelMap::IsDead(int detID, int channel, int ringwedgeFlag) {
-	if(!IsValid() || (ringwedgeFlag != 0 && ringwedgeFlag != 1)) {
-		std::cerr<<"Error at DeadChannelMap IsDead(), bad input parameters"<<std::endl;
-		return false;
-	}
-
-	int this_channel = detID*24+ringwedgeFlag*16+channel;
-
-	auto iter = dcMap.find(this_channel);
-	if(iter != dcMap.end()) {
-		return iter->second;
-	} else {
-		std::cerr<<"Error at DeadChannelMap IsDead(), invalid channel: "<<this_channel<<std::endl;
-		std::cerr<<"detID: "<<detID<<" ringwedgeFlag: "<<ringwedgeFlag<<" channel: "<<channel<<std::endl;
-		return false;
-	}
-}

src/Detectors/DetectorEfficiency.cpp

@@ -17,7 +17,7 @@ int main(int argc, char** argv) {
 
 	/*
 	SabreEfficiency sabre;
-	std::string mapfile = "./etc/DeadChannels.txt";
+	std::string mapfile = "./etc/SabreDeadChannels.txt";
 	sabre.SetDeadChannelMap(mapfile);
 	sabre.CalculateEfficiency(inputname, outputname, statsname);
 	//std::cout<<"Running consistency check(1=success): "<<sabre.RunConsistencyCheck()<<std::endl;
@@ -26,6 +26,8 @@ int main(int argc, char** argv) {
 
 	try {
  		AnasenEfficiency anasen;
+ 		std::string mapfile = "./etc/AnasenDeadChannels.txt";
+ 		anasen.SetDeadChannelMap(mapfile);
  		anasen.CalculateEfficiency(inputname, outputname, statsname);
  	//std::cout<<"Running consistency check(1=success): "<<anasen.RunConsistencyCheck()<<std::endl;
  	//anasen.DrawDetectorSystem("/data1/gwm17/TRIUMF_7Bed/simulation/ANASENGeo_centered_target_targetGap_BackQQQ_fixedZ.txt");

src/Detectors/SabreDeadChannelMap.cpp

@@ -0,0 +1,80 @@
+#include "SabreDeadChannelMap.h"
+#include <fstream>
+
+SabreDeadChannelMap::SabreDeadChannelMap() :
+	validFlag(false)
+{
+	int maxchan = 5*24+1*16+7;
+	for(int i=0; i<maxchan; i++)
+		dcMap[i] = false;
+}
+
+SabreDeadChannelMap::SabreDeadChannelMap(const std::string& name) :
+	validFlag(false)
+{
+	LoadMapfile(name);
+}
+
+SabreDeadChannelMap::~SabreDeadChannelMap() {}
+
+void SabreDeadChannelMap::LoadMapfile(const std::string& name) 
+{
+	std::ifstream input(name);
+	if(!input.is_open()) 
+	{
+		std::cerr<<"Unable to load dead channels, file: "<<name<<" could not be opened"<<std::endl;
+		validFlag = false;
+		return;
+	}
+
+	
+
+	std::string junk, rw;
+	int detID, channel;
+	int this_channel;
+
+	std::getline(input, junk);
+	while(input>>detID) 
+	{
+		input>>rw>>channel;
+		if(rw == "RING") 
+		{
+			this_channel = detID*24+RING*16+channel;
+			dcMap[this_channel] = true;
+		}
+		else if(rw == "WEDGE")
+		{
+			this_channel = detID*24+WEDGE*16+channel;
+			dcMap[this_channel] = true;
+		}
+		else 
+		{
+			std::cerr<<"Invalid ring/wedge designation at SabreDeadChannelMap"<<std::endl;
+		}
+	}
+
+	validFlag = true;
+}
+
+bool SabreDeadChannelMap::IsDead(int detID, int channel, int ringwedgeFlag) 
+{
+	if(!IsValid() || (ringwedgeFlag != 0 && ringwedgeFlag != 1)) 
+	{
+		std::cerr<<"Error at SabreDeadChannelMap IsDead(), bad input parameters"<<std::endl;
+		return false;
+	}
+
+	int this_channel = detID*24+ringwedgeFlag*16+channel;
+
+	auto iter = dcMap.find(this_channel);
+	if(iter != dcMap.end()) 
+	{
+		return iter->second;
+	} 
+	else 
+	{
+		std::cerr<<"Error at SabreDeadChannelMap IsDead(), invalid channel: "<<this_channel<<std::endl;
+		std::cerr<<"detID: "<<detID<<" ringwedgeFlag: "<<ringwedgeFlag<<" channel: "<<channel<<std::endl;
+		return false;
+	}
+}

src/Detectors/StripDetector.cpp

@@ -53,8 +53,8 @@ StripDetector::~StripDetector() {}
 void StripDetector::CalculateCorners() {
 	double y_min, y_max, z_min, z_max; 
 	for (int s=0; s<num_strips; s++) {
-		y_max = -total_width/2.0 + front_strip_width*(s+1);
+		y_max = total_width/2.0 - front_strip_width*s;
-		y_min = -total_width/2.0 + front_strip_width*s;
+		y_min = total_width/2.0 - front_strip_width*(s+1);
 		z_max = center_z + length/2.0;
 		z_min = center_z - length/2.0;
 		front_strip_coords[s][2] = Mask::Vec3(center_rho, y_max, z_max);
@@ -104,8 +104,9 @@ Mask::Vec3 StripDetector::GetHitCoordinates(int front_stripch, double front_stri
 
 }
 
-std::pair<int,double> StripDetector::GetChannelRatio(double theta, double phi) {
+StripHit StripDetector::GetChannelRatio(double theta, double phi) {
 
+	StripHit hit;
 	while (phi < 0) phi += 2*M_PI;
 
 	//to make the math easier (and the same for each det), rotate the input phi
@@ -118,7 +119,7 @@ std::pair<int,double> StripDetector::GetChannelRatio(double theta, double phi) {
 	double det_max_phi = atan2(total_width/2,center_rho);
 	double det_min_phi = -det_max_phi;
   
-	if (phi < det_min_phi || phi > det_max_phi) return std::make_pair(-1,0);
+	if (phi < det_min_phi || phi > det_max_phi) return hit;
 
 	//for theta it's not so simple, so we have to go through the typical plane-intersect method
 	//first thing's first: we have a fixed x for the entire detector plane:
@@ -127,18 +128,26 @@ std::pair<int,double> StripDetector::GetChannelRatio(double theta, double phi) {
 	double yhit = xhit*tan(phi);
 	double zhit = sqrt(xhit*xhit+yhit*yhit)/tan(theta);
 
-	int front_ch_hit=-1;
-	double ratio=0;
 	for (int s=0; s<num_strips; s++) {
 		if (xhit >= front_strip_coords[s][0].GetX() && xhit <= front_strip_coords[s][0].GetX() && //Check min and max x (constant in flat)
 			yhit >= front_strip_coords[s][1].GetY() && yhit <= front_strip_coords[s][2].GetY() && //Check min and max y
 			zhit >= front_strip_coords[s][1].GetZ() && zhit <= front_strip_coords[s][0].GetZ()) //Check min and max z
 		{
-			front_ch_hit = s;
+			hit.front_strip_index = s;
-			ratio = (zhit-center_z)/(length/2);
+			hit.front_ratio = (zhit-center_z)/(length/2);
 			break;
 		}
 	}
 
-	return std::make_pair(front_ch_hit,ratio);
+	for (int s=0; s<num_strips; s++) {
+		if (xhit >= back_strip_coords[s][0].GetX() && xhit <= back_strip_coords[s][0].GetX() && //Check min and max x (constant in flat)
+			yhit >= back_strip_coords[s][1].GetY() && yhit <= back_strip_coords[s][2].GetY() && //Check min and max y
+			zhit >= back_strip_coords[s][1].GetZ() && zhit <= back_strip_coords[s][0].GetZ()) //Check min and max z
+		{
+			hit.back_strip_index = s;
+			break;
+		}
+	}
+
+	return hit;
 }