Reworking of the ObserverTimeEvolution

include/crpropa/module/Observer.h

@@ -241,12 +241,21 @@ class ObserverParticleIdVeto: public ObserverFeature {
  */
 class ObserverTimeEvolution: public ObserverFeature {
 private:
-	std::vector<double> detList;
+	int numb;
+	bool islog = false, useCustomGetTime = false;
+	double min, max;
 public:
+	// List containing all used times, is constructed in constructor or by user manually
+	// (leave empty if you want to rather use functions)
+	std::vector<double> detList;
+
+	// Pointer to user defined version of getTime (should return time for corresponding index)
+	double (*customGetTime)(int index, double min, double max, int numb);
+
 	/** Default constructor
 	 */
 	ObserverTimeEvolution();
-	/** Constructor
+	/** Constructor (calculates the maximum from min + dist * numb)
 	 @param min		minimum time
 	 @param dist	time interval for detection
 	 @param numb	number of time intervals
@@ -259,12 +268,23 @@ class ObserverTimeEvolution: public ObserverFeature {
 	 @param log     log (input: true) or lin (input: false) scaling between min and max with numb steps
 	 */
 	ObserverTimeEvolution(double min, double max, double numb, bool log);
-	// Add a new time step to the detection time list of the observer
-	void addTime(const double &position);
-	// Using log or lin spacing of times in the range between min and
-	// max for observing particles
-	void addTimeRange(double min, double max, double numb, bool log = false);
-	const std::vector<double>& getTimes() const;
+	/** Constructor
+	 @param detList	user defined vector<double> with times to check
+	 */
+	ObserverTimeEvolution(const std::vector<double> &detList);
+	// Destructor
+	~ObserverTimeEvolution(){}
+
+	/** Function to set a user defined function, which takes index, min, max and numb to calculate the time at index
+	 @param userDefinedFunction	function which will be called instead of getTime, should return time for corresponding index
+	 	min time, max time, and number of times
+	 */
+	void setUserDefinedGetTime(double (*userDefinedFunction)(int index, double min, double max, int numb));
+	// Get time at specified index, either with log or lin spacing of times
+	// in the range between min and max for obeserving particles.
+	double getTime(std::size_t index) const;
+	// Return a vector containing all times between min and max generated by getTime
+	const std::vector<double>& getTimes();
 	DetectionState checkDetection(Candidate *candidate) const;
 	std::string getDescription() const;
 };

src/module/Observer.cpp

@@ -246,23 +246,37 @@ std::string ObserverParticleIdVeto::getDescription() const {
 ObserverTimeEvolution::ObserverTimeEvolution() {}
 
 ObserverTimeEvolution::ObserverTimeEvolution(double min, double dist, double numb) {
-	double max = min + numb * dist;
-	bool log = false;
-	addTimeRange(min, max, numb, log);
+	this->max = min + numb * dist;
+	this->min = min;
+	this->numb = numb;
+	this->islog = false;
 }
 
 ObserverTimeEvolution::ObserverTimeEvolution(double min, double max, double numb, bool log) {
-	addTimeRange(min, max, numb, log);
+	this->min = min;
+	this->max = max;
+	this->numb = numb;
+	this->islog = log;
 }
 
+ObserverTimeEvolution::ObserverTimeEvolution(const std::vector<double> &detList){
+	this->detList = detList;
+	this->numb = detList.size();
+	this->min = detList.front();
+	this->max = detList.back();
+}
+
+void ObserverTimeEvolution::setUserDefinedGetTime(double (*userDefinedFunction)(int index, double min, double max, int numb)){
+	this->customGetTime = userDefinedFunction;
+	this->useCustomGetTime = true;
+}
 
 DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
 
-	if (detList.size()) {
+	if (numb) {
 		double length = c->getTrajectoryLength();
 		size_t index;
 		const std::string DI = "DetectionIndex";
-		std::string value;
 
 		// Load the last detection index
 		if (c->hasProperty(DI)) {
@@ -272,13 +286,13 @@ DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
 			index = 0;
 		}
 
-		// Break if the particle has been detected once for all detList entries.
-		if (index > detList.size()) {
+		// Break if the particle has been detected once for all possible times.
+		if (index > numb) {
 			return NOTHING;
 		}
 
 		// Calculate the distance to next detection
-		double distance = length - detList[index];
+		double distance = length - getTime(index);
 
 		// Limit next step and detect candidate.
 		// Increase the index by one in case of detection
@@ -288,8 +302,8 @@ DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
 		}
 		else {
 
-			if (index < detList.size()-1) {
-				c->limitNextStep(detList[index+1]-length);
+			if (index < numb-1) {
+				c->limitNextStep(getTime(index+1)-length);
 			}
 			c->setProperty(DI, Variant::fromUInt64(index+1));
 
@@ -299,29 +313,34 @@ DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
 	return NOTHING;
 }
 
-void ObserverTimeEvolution::addTime(const double& t) {
-	detList.push_back(t);
+double ObserverTimeEvolution::getTime(size_t index) const {
+	if (!detList.empty()) {
+		return detList[index];
+	} else if (useCustomGetTime) {
+		return customGetTime(index, min, max, numb);
+	} else if (log) {
+		return min * pow(max / min, index / (numb - 1.0));
+	} else {
+		return min + index * (max - min) / (numb - 1.0);
+	}
 }
 
-void ObserverTimeEvolution::addTimeRange(double min, double max, double numb, bool log) {
-	for (size_t i = 0; i < numb; i++) {
-		if (log == true) {
-			addTime(min * pow(max / min, i / (numb - 1.0)));
-		} else {
-			addTime(min + i * (max - min) / numb);
+const std::vector<double>& ObserverTimeEvolution::getTimes() {
+	if (detList.empty()) {
+		size_t counter = 0;
+		while (getTime(counter)<=max) {
+			detList.push_back(getTime(counter));
+			counter++;
 		}
 	}
-}
-
-const std::vector<double>& ObserverTimeEvolution::getTimes() const {
 	return detList;
 }
 
 std::string ObserverTimeEvolution::getDescription() const {
 	std::stringstream s;
 	s << "List of Detection lengths in kpc";
-	for (size_t i = 0; i < detList.size(); i++)
-	  s << "  - " << detList[i] / kpc;
+	for (size_t i = 0; i < numb; i++)
+	  s << "  - " << getTime(i) / kpc;
 	return s.str();
 }
 

test/testBreakCondition.cpp

@@ -236,6 +236,7 @@ TEST(ObserverFeature, TimeEvolution) {
   Observer obs;
   obs.setDeactivateOnDetection(false);
   obs.setFlag("Detected", "Detected");
+  //min = 5, max = min + numb*dist = 5 + 2*5 = 15, detection can happen at [5, 15]
   obs.add(new ObserverTimeEvolution(5, 5, 2));
   Candidate c;
   c.setNextStep(10);
@@ -261,7 +262,7 @@ TEST(ObserverFeature, TimeEvolution) {
 
   // detection two
   c.setCurrentStep(0.1);
-  c.setTrajectoryLength(10.05);
+  c.setTrajectoryLength(15.05);
   obs.process(&c);
   EXPECT_TRUE(c.isActive());
   EXPECT_TRUE(c.hasProperty("Detected"));