Merge pull request #103 from marc1uk/LightMap2

MCParticleProperties bugfixes

DataModel/Detector.cpp

@@ -0,0 +1,11 @@
+/* vim:set noexpandtab tabstop=4 wrap */
+#include "Detector.h"
+#include "Geometry.h"
+
+Position Detector::GetPositionInTank(){
+	if((DetectorElement!="Tank")||(GeometryPtr==nullptr)){
+		return DetectorPosition;
+	} else {
+		return (DetectorPosition-GeometryPtr->GetTankCentre());
+	}
+}

DataModel/Detector.h

@@ -2,10 +2,14 @@
 #ifndef DETECTORCLASS_H
 #define DETECTORCLASS_H
 
+#include <math.h> // defines M_PI
+//constexpr double PI=4.*atan(1.); // can also do this
+
 #include <SerialisableObject.h>
 #include "Position.h"
 #include "Direction.h"
 #include "Channel.h"
+class Geometry;
 
 enum class detectorstatus : uint8_t { OFF, ON, UNSTABLE };
 
@@ -14,27 +18,34 @@ class Detector : public SerialisableObject{
 	friend class boost::serialization::access;
 
 	public:
-	Detector() : DetectorID(0), DetectorElement(), DetectorPosition(), DetectorDirection(), DetectorType(""),
+	Detector() : DetectorID(0), DetectorElement(), TankLocation(), DetectorPosition(), DetectorDirection(), DetectorType(""),
 		Status(detectorstatus::OFF), Channels() {serialise=true;}
-	Detector(int detid, std::string DetEle, Position posin, Direction dirin, std::string detype, detectorstatus stat, double avgrate, map<unsigned long,Channel> channelsin={}) : DetectorID(detid), DetectorElement(DetEle), DetectorPosition(posin), DetectorDirection(dirin), DetectorType(detype), Status(stat), Channels(channelsin) {serialise=true;}
+	Detector(int detid, std::string DetEle, std::string CylLoc, Position posin, Direction dirin, std::string detype, detectorstatus stat, double avgrate, map<unsigned long,Channel> channelsin={}) : DetectorID(detid), DetectorElement(DetEle), TankLocation(CylLoc), DetectorPosition(posin), DetectorDirection(dirin), DetectorType(detype), Status(stat), Channels(channelsin) { serialise=true; }
+
 	std::string GetDetectorElement(){return DetectorElement;}
 	Position GetDetectorPosition(){return DetectorPosition;}
+	Position GetPositionInTank();
 	Direction GetDetectorDirection(){return DetectorDirection;}
 	int GetDetectorID(){return static_cast<int>(DetectorID);}
 	std::string GetDetectorType(){return DetectorType;}
 	detectorstatus GetStatus(){return Status;}
 	std::map<unsigned long,Channel>* GetChannels() {return &Channels;}
 	void AddChannel(Channel chanin){ Channels.emplace(chanin.GetChannelID(),chanin); }
+	std::string GetTankLocation(){ return TankLocation; }
+	Geometry* GetGeometryPtr(){ return GeometryPtr; }
 
 	void SetDetectorElement(std::string DetEleIn){DetectorElement=DetEleIn;}
 	void SetDetectorPosition(Position DetectorPositionIn){DetectorPosition=DetectorPositionIn;}
 	void SetDetectorDirection(Direction DetectorDirectionIn){DetectorDirection=DetectorDirectionIn;}
 	void SetDetectorID(int DetectorIDIn){DetectorID=DetectorIDIn;}
 	void SetDetectorType(std::string DetectorTypeIn){DetectorType=DetectorTypeIn;}
 	void SetStatus(detectorstatus StatusIn){Status=StatusIn;}
+	void SetTankLocation(std::string locin){TankLocation=locin;}
+	void SetGeometryPtr(Geometry* geomin){ GeometryPtr=geomin; }
 	bool Print(){
 		std::cout<<"DetectorPosition  : "; DetectorPosition.Print();
 		std::cout<<"DetectorDirection : "; DetectorDirection.Print();
+		std::cout<<"Location          : "<<TankLocation<<std::endl;
 		std::cout<<"DetectorElement   : "<<DetectorElement<<std::endl;
 		std::cout<<"DetectorID        : "<<DetectorID<<std::endl;
 		std::cout<<"DetectorType      : "<<DetectorType<<std::endl;
@@ -64,13 +75,16 @@ class Detector : public SerialisableObject{
 	std::string DetectorType;      // e.g. "Hamamatsu R7081"
 	detectorstatus Status;         // on, off, unstable....
 	std::map<unsigned long,Channel> Channels;
+	std::string TankLocation;      // "Barrel", "TopCap", "BottomCap", "MRD"*, "FACC"*, or "NA". *may change
+	Geometry* GeometryPtr=nullptr; // a pointer to the parent geometry to which this Detector belongs
 
 	template<class Archive> void serialize(Archive & ar, const unsigned int version){
 		if(serialise){
 		ar & DetectorElement;
 		ar & Channels;
 		ar & DetectorPosition;
 		ar & DetectorDirection;
+		ar & TankLocation;
 		ar & DetectorID;
 		ar & DetectorType;
 		ar & Status;

DataModel/Direction.h

@@ -20,26 +20,26 @@ class Direction : public SerialisableObject{
 		double mag=sqrt(pow(xin,2.)+pow(yin,2.)+pow(zin,2.));
 		if(mag==1.){ x=xin; y=yin; z=zin; }
 		else { x=xin/mag; y=yin/mag; z=zin/mag; }
-		theta = (x == 0.0 && y == 0.0 && z == 0.0) ? 0.0 : atan2(sqrt(x*x + y*y),z);
-		phi = (x == 0.0 && y == 0.0) ? 0. : atan2(y, x);
+		phi = (x == 0.0 && y == 0.0 && z == 0.0) ? 0.0 : atan2(sqrt(x*x + y*y),z);
+		theta = (x == 0.0 && y == 0.0) ? 0. : atan2(y, x);
 	};
 
 	Direction(double phiin, double thetain){
 		serialise=true;
 		theta=thetain;
 		phi=phiin;
-		x=cos(theta)*cos(phi);
-		y=cos(theta)*sin(phi);
-		z=sin(theta);
+		x=cos(phi)*cos(theta);
+		y=cos(phi)*sin(theta);
+		z=sin(phi);
 	};
 
 	inline double X() const {return x;}
 	inline double Y() const {return y;}
 	inline double Z() const {return z;}
-	inline double Phi() const {return phi;}
-	inline double PhiDeg() const {return phi*(180./M_PI);}
-	inline double Theta() const {return theta;}
-	inline double ThetaDeg() const {return theta*(180./M_PI);}
+	inline double GetPhi() const {return phi;}
+	inline double GetPhiDeg() const {return phi*(180./M_PI);}
+	inline double GetTheta() const {return theta;}
+	inline double GetThetaDeg() const {return theta*(180./M_PI);}
 
 	inline void SetX(double xx){x=xx;}
 	inline void SetY(double yy){y=yy;}
@@ -58,11 +58,11 @@ class Direction : public SerialisableObject{
 	}
 
 	private:
-	double x;      // meters
+	double x;        // meters
 	double y;
 	double z;
-	double theta;  // relative to beam direction (z), rads
-	double phi;    // relative to x axis, rads
+	double phi;      // clockwise looking down, 0 pointing downstream along beam, rads
+	double theta;    // pitch angle, relative to the x-z (beamline) plane, rads
 
 	template<class Archive> void serialize(Archive & ar, const unsigned int version){
 		if(serialise){

DataModel/Geometry.cpp

@@ -96,3 +96,18 @@ void Geometry::PrintChannels(){
 		}
 	}
 }
+
+void Geometry::CartesianToPolar(Position posin, double& R, double& Phi, double& Theta, bool tankcentered){
+	// Calculate angle from beam axis, measured clockwise while looking down
+	// first shift to place relative to the tank origin if needed
+	if(not tankcentered){ posin -= tank_centre; }
+	// calculate the angle from the beam axis
+	double thethetaval = atan(posin.X()/abs(posin.Z()));
+	if(posin.Z()<0.){ (posin.X()<0.) ? thethetaval=(-M_PI+thethetaval) : thethetaval=(M_PI-thethetaval); }
+	Phi = thethetaval;
+	// calculate angle from the x-z plane
+	Theta = atan(posin.Y() / sqrt(pow(posin.X(),2.)+pow(posin.Z(),2.)));
+	// calculate the radial distance from the tank centre
+	R = sqrt(pow(posin.X(),2.)+pow(posin.Z(),2.));
+	return;
+}

DataModel/Geometry.h

@@ -72,6 +72,10 @@ class Geometry : public SerialisableObject{
 	}
 
 	bool AddDetector(Detector detin){
+		// Pass a pointer to it's owning geometry to this Detector
+		// we need to do this before calling `emplace` as that must do a copy-construction
+		detin.SetGeometryPtr(this);
+
 		std::string thedetel = detin.GetDetectorElement();
 		int detectorsetindex=-1;
 		if(DetectorElements.count(thedetel)==0){
@@ -249,6 +253,10 @@ class Geometry : public SerialisableObject{
 
 	void PrintChannels();
 
+	// helper functions
+	Position GlobalToTankCentered(Position posin){ return (posin - tank_centre); }
+	void CartesianToPolar(Position posin, double& R, double& Phi, double& Theta, bool tankcentered=false);
+
 	private:
 	unsigned long NextFreeDetectorKey;
 	unsigned long NextFreeChannelKey;

DataModel/Hit.h

@@ -7,33 +7,33 @@
 using namespace std;
 
 class Hit : public SerialisableObject{
-
+	
 	friend class boost::serialization::access;
-
+	
 	public:
-  Hit() : TubeId(0), Time(0), Charge(0){serialise=true;}
-  Hit(int tubeid, double thetime, double charge) : TubeId(tubeid), Time(thetime), Charge(charge){serialise=true;}
-
+	Hit() : TubeId(0), Time(0), Charge(0){serialise=true;}
+	Hit(int thetubeid, double thetime, double thecharge) : TubeId(thetubeid), Time(thetime), Charge(thecharge){serialise=true;}
+	
 	inline int GetTubeId() const {return TubeId;}
 	inline double GetTime() const {return Time;}
 	inline double GetCharge() const {return Charge;}
-
+	
 	inline void SetTubeId(int tubeid){TubeId=tubeid;}
 	inline void SetTime(double tc){Time=tc;}
 	inline void SetCharge(double chg){Charge=chg;}
-
+	
 	bool Print() {
 		cout<<"TubeId : "<<TubeId<<endl;
 		cout<<"Time : "<<Time<<endl;
 		cout<<"Charge : "<<Charge<<endl;
 		return true;
 	}
-
+	
 	protected:
 	int TubeId;
 	double Time;
 	double Charge;
-
+	
 	template<class Archive> void serialize(Archive & ar, const unsigned int version){
 		if(serialise){
 			ar & TubeId;
@@ -43,22 +43,58 @@ class Hit : public SerialisableObject{
 	}
 };
 
-/*  Derived classes, if there's a reason to have them. So far...not really
+//  Derived classes
 
-class TDCHit : public Hit {
+class MCHit : virtual public Hit {
+
+	friend class boost::serialization::access;
+
 	public:
+	MCHit() : Hit(), Parents(std::vector<int>{}) {serialise=true;}
+	MCHit(int tubeid, double thetime, double thecharge, std::vector<int> theparents) : Hit(tubeid, thetime, thecharge), Parents(theparents) {serialise=true;}
+
+	const std::vector<int>* GetParents() const { return &Parents; }
+	void SetParents(std::vector<int> parentsin){ Parents = parentsin; }
+
+	bool Print(){
+		cout<<"TubeId : "<<TubeId<<endl;
+		cout<<"Time : "<<Time<<endl;
+		cout<<"Charge : "<<Charge<<endl;
+		if(Parents.size()){
+			cout<<"Parent MCPartice indices: {";
+			for(int parenti=0; parenti<Parents.size(); ++parenti){
+				cout<<Parents.at(parenti);
+				if((parenti+1)<Parents.size()) cout<<", ";
+			}
+			cout<<"}"<<endl;
+		} else {
+			cout<<"No recorded parents"<<endl;
+		}
+		return true;
+	}
+
+	protected:
+	std::vector<int> Parents;
+
+	template<class Archive> void serialize(Archive & ar, const unsigned int version){
+		if(serialise){
+			ar & TubeId;
+			ar & Time;
+			ar & Charge;
+			// do not serialize parents; the indices by themselves are not meaningful
+		}
+	}
+};
 
-	private:
-}
-
+/*
 class RecoHit : public Hit {
 	public:
 	RecoHit(double thetime, double thecharge) : Time(thetime), Charge(thecharge){};
 
 	inline double GetCharge(){return Charge;}
 	inline void SetCharge(double chg){Charge=chg;}
 
-	private:
+	protected:
 	double Charge;
 };
 */

DataModel/LAPPDHit.h

@@ -5,21 +5,22 @@
 #include<Hit.h>
 #include<SerialisableObject.h>
 
-using namespace std;
-
-class LAPPDHit : public Hit{
+using std::cout;
+using std::endl;
 
+class LAPPDHit : virtual public Hit{
+
 	friend class boost::serialization::access;
-
+	
 	public:
 	LAPPDHit() : Hit(), Position(0), LocalPosition(0) {serialise=true;}
-	LAPPDHit(int tubeid, double thetime, double charge, std::vector<double> Position, std::vector<double> LocalPosition) : Hit(tubeid,thetime,charge), Position(Position), LocalPosition(LocalPosition) {serialise=true;}
-
+	LAPPDHit(int thetubeid, double thetime, double thecharge, std::vector<double> theposition, std::vector<double> thelocalposition) : Hit(thetubeid,thetime,thecharge), Position(theposition), LocalPosition(thelocalposition) {serialise=true;}
+	
 	inline std::vector<double> GetPosition() const {return Position;}
 	inline std::vector<double> GetLocalPosition() const {return LocalPosition;}
 	inline void SetPosition(std::vector<double> pos){Position=pos;}
 	inline void SetLocalPosition(std::vector<double> locpos){LocalPosition=locpos;}
-
+	
 	bool Print() {
 		cout<<"TubeId : "<<TubeId<<endl;
 		cout<<"Time : "<<Time<<endl;
@@ -31,25 +32,76 @@ class LAPPDHit : public Hit{
 		cout<<"Charge : "<<Charge<<endl;
 		return true;
 	}
-
+	
 	protected:
 	std::vector<double> Position;
 	std::vector<double> LocalPosition;
+
+	template<class Archive> void serialize(Archive & ar, const unsigned int version){
+		if(serialise){
+			ar & TubeId;
+			ar & Time;
+			ar & Position;
+			ar & LocalPosition;
+			ar & Charge;
+		}
+	}
+};
 
+//  Derived classes, if there's a reason to have them
 
+class MCLAPPDHit : public MCHit, public LAPPDHit{
+
+	friend class boost::serialization::access;
+
+	public:
+	MCLAPPDHit() : LAPPDHit(), MCHit() {serialise=true;}
+	MCLAPPDHit(int thetubeid, double thetime, double thecharge, std::vector<double> theposition, std::vector<double> thelocalposition, std::vector<int> theparents){
+		TubeId = thetubeid;
+		Time=thetime;
+		Charge=thecharge;
+		Position=theposition;
+		LocalPosition=thelocalposition;
+		Parents=theparents;
+		serialise=true;
+	}
+
+	bool Print() {
+		cout<<"TubeId : "<<TubeId<<endl;
+		cout<<"Time : "<<Time<<endl;
+		cout<<"X Pos : "<<Position.at(0)<<endl;
+		cout<<"Y Pos : "<<Position.at(1)<<endl;
+		cout<<"Z Pos : "<<Position.at(2)<<endl;
+		cout<<"Parallel Pos : "<<LocalPosition.at(0)<<endl;
+		cout<<"Transverse Pos : "<<LocalPosition.at(1)<<endl;
+		cout<<"Charge : "<<Charge<<endl;
+		if(Parents.size()){
+			cout<<"Parent MCPartice indices: {";
+			for(int parenti=0; parenti<Parents.size(); ++parenti){
+				cout<<Parents.at(parenti);
+				if((parenti+1)<Parents.size()) cout<<", ";
+			}
+			cout<<"}"<<endl;
+		} else {
+			cout<<"No recorded parents"<<endl;
+		}
+		return true;
+	}
+
 	template<class Archive> void serialize(Archive & ar, const unsigned int version){
 		if(serialise){
 			ar & TubeId;
 			ar & Time;
 			ar & Position;
 			ar & LocalPosition;
 			ar & Charge;
+			// n.b. at time of writing MCHit stores no additional persistent members
+			// - it only adds parent MCParticle indices, and these aren't saved... 
 		}
 	}
 };
 
-/*  Derived classes, if there's a reason to have them. So far...not really
-
+/*
 class TDCHit : public Hit {
 	public: