Reworking of the ObserverTimeEvolution

include/crpropa/module/Observer.h

@@ -240,8 +240,22 @@ class ObserverParticleIdVeto: public ObserverFeature {
  This observer is very useful if the time evolution of the particle density is needed. It detects all candidates in lin-spaced, log-spaced, or user-defined time intervals and limits the nextStep of candidates to prevent overshooting of detection intervals.
  */
 class ObserverTimeEvolution: public ObserverFeature {
-private:
+protected:
+	int numb;
+	bool islog = false, doDetListConstruction = true;
+	double min, max;
+	/** Vector containing all used times. 
+	 It is only constructed by the user manually.
+	 If it is not empty, the vector will be used instead of the getTime function.
+	 (leave empty if you want to rather use functions)
+	*/
 	std::vector<double> detList;
+	/**
+	 A temporary storage for detList, this enables the return of a List in getTimes
+	 without risking to modify detList
+	 */
+	mutable std::vector<double> tempDetList;
+
 public:
 	/** Default constructor
 	 */
@@ -250,22 +264,105 @@ class ObserverTimeEvolution: public ObserverFeature {
 	 @param min		minimum time
 	 @param dist	time interval for detection
 	 @param numb	number of time intervals
+
+	 This constructor calculates the maximum from max = min + (numb - 1) * dist
 	 */
 	ObserverTimeEvolution(double min, double dist, double numb);
 	/** Constructor
 	 @param min		minimum time
 	 @param max	    maximum time
 	 @param numb	number of time intervals
 	 @param log     log (input: true) or lin (input: false) scaling between min and max with numb steps
+
+	 This constructor sets the maximum directly and gets numb automatically.
+	 You need to set the log parameter, since an overload for the first three doubles exist.
 	 */
 	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
+	/** Constructor
+	 @param detList	user defined vector<double> with times to check
+
+	 This constructor uses a predefined vector containing the times that should be observed.
+	 The so created detList can then be modified via addTime, addTimeRange and setTimes.
+	 */
+	ObserverTimeEvolution(const std::vector<double> &detList);
+	/** Destructor
+	 */
+	~ObserverTimeEvolution(){}
+
+	/** Function
+	 Generates the detList if it is empty when for example the 
+	 ObserverTimeEvolution(const std::vector<double> &detList) constructor
+	 was not used.
+	 Use this function to create a detList with can then be modified.
+	 When detList is not empty its entries are used instead of a runtime calculation of the times.	 
+	 */
+	void constructDetListIfEmpty();
+	/** Function
+	 @param candidate	Candidate usally given by a module list
+
+	 Checks whether to make a detection at the current step of candidate or not.
+	 This function is called in Observer.process with the simulated Candidate.
+	 */
+	DetectionState checkDetection(Candidate *candidate) const;
+	/** Function
+	 @param enableConstruction	if true, constructs detList from range of min, max, numb
+	 when calling addTime
+	 Clears the content of detList
+	 */
+	void clear();
+	/** Function
+	 Checks if detList is empty
+	 */
+	bool empty(){return detList.empty();}
+
+	// setter functions:
+	/** Function
+	 @param time	Time that should be appended to detList
+
+	 Makes a push_back on detList with the given time.	 
+	 */
+	void addTime(const double &time);
+	/** Function
+	 @param min		minimum time
+	 @param max	    maximum time
+	 @param numb	number of time intervals
+	 @param log     log (input: true) or lin (input: false) scaling between min and max with numb steps
+
+	 Appends a linear or logarithmic time range to detList via repeatedly calling addTime
+	 */
 	void addTimeRange(double min, double max, double numb, bool log = false);
+	/** Function
+	 @param detList	vector<double> containing times when to observe
+
+	 Sets this->detList to detList, with this it is possible to fully modify detList
+	 */
+	void setTimes(const std::vector<double> &detList);
+	void setMin(double min){this->min = min;}
+	void setMax(double max){this->max = max;}
+	void setNumb(int numb){this->numb = numb;}
+	void setIslog(bool log){this->islog = islog;}
+
+	// getter functions:
+	/** Function
+	 @param index	index of the required time
+
+	 Replaces the previous preconstructed detList and return the time for a specific index
+	 */
+	virtual double getTime(std::size_t index) const;
+	double getMin() const {return min;}
+	double getMax() const {return max;}
+	int getNumb() const {return numb;}
+	bool getIsLog() const {return islog;}
+	/** Function 
+	 Returns a vector<double> containing all times between min and max generated by getTime.
+	 This function does not return detList directly, it rather appends a new vector with
+	 getTime numb times.
+	 */
 	const std::vector<double>& getTimes() const;
-	DetectionState checkDetection(Candidate *candidate) const;
+	/** Function
+	 Returns a string containing a representation of all times.
+	 This function does not create a detList.
+	 */
 	std::string getDescription() const;
 };
 /** @} */

src/module/Observer.cpp

@@ -246,23 +246,29 @@ 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 - 1) * 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){
+	setTimes(detList);
+}
 
 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 +278,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 +294,8 @@ DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
 		}
 		else {
 
-			if (index < detList.size()-1) {
-				c->limitNextStep(detList[index+1]-length);
+			if (index < numb-2) {
+				c->limitNextStep(getTime(index+1)-length);
 			}
 			c->setProperty(DI, Variant::fromUInt64(index+1));
 
@@ -299,29 +305,81 @@ DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
 	return NOTHING;
 }
 
-void ObserverTimeEvolution::addTime(const double& t) {
-	detList.push_back(t);
+void ObserverTimeEvolution::clear(){
+	doDetListConstruction = false;
+	detList.clear();
+	numb = 0;
+}
+
+void ObserverTimeEvolution::constructDetListIfEmpty(){
+	if (detList.empty() && doDetListConstruction) {
+		std::vector<double> detListTemp;
+		size_t counter = 0;
+		while (getTime(counter)<=max) {
+			detListTemp.push_back(getTime(counter));
+			counter++;
+		}
+		detList.assign(detListTemp.begin(), detListTemp.end());
+	}
+}
+
+void ObserverTimeEvolution::addTime(const double &time){
+	constructDetListIfEmpty();
+	detList.push_back(time);
+	numb += 1;  // increase number of entries by one
 }
 
 void ObserverTimeEvolution::addTimeRange(double min, double max, double numb, bool log) {
 	for (size_t i = 0; i < numb; i++) {
-		if (log == true) {
+		if (log) {
+			if ( min == 0 ){
+				std::cout << "min can not be 0 if log=true" << std::endl;
+				throw new std::invalid_argument("min can not be 0 if log=true");
+			}
 			addTime(min * pow(max / min, i / (numb - 1.0)));
 		} else {
-			addTime(min + i * (max - min) / numb);
+			addTime(min + i * (max - min) / (numb - 1.0));
 		}
 	}
+	// allready corrected by addTime, just here to be safe
+	this->numb = detList.size();
+}
+
+void ObserverTimeEvolution::setTimes(const std::vector<double> &detList){
+	this->detList.assign(detList.begin(), detList.end());
+	this->numb = detList.size();
+	this->min = detList.front();
+	this->max = detList.back();
+	this->doDetListConstruction = false;
+}
+
+double ObserverTimeEvolution::getTime(size_t index) const {
+	if (!detList.empty()) {
+		return detList.at(index);
+	} else if (islog) {
+		if ( min == 0 ){
+			std::cout << "min can not be 0 if islog=true" << std::endl;
+			throw new std::invalid_argument("min can not be 0 if islog=true");
+		}
+		return min * pow(max / min, index / (numb - 1.0));
+	} else {
+		return min + index * (max - min) / (numb - 1.0);
+	}
 }
 
 const std::vector<double>& ObserverTimeEvolution::getTimes() const {
-	return detList;
+	tempDetList.clear();
+	for (size_t i = 0; i < numb; i++){
+		tempDetList.push_back(getTime(i));
+	}
+	return tempDetList;
 }
 
 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();
 }