Merge pull request #511 from JanNiklasB/RenewObserverTimeEvolution

Reworking of the ObserverTimeEvolution

Author: JulienDoerner

CHANGELOG.md

@@ -1,14 +1,17 @@
 ## CRPropa vNext
 
 ### Bug fixes:
- * Fixed sign for exponential decay of magn. field strength with Galactic height in LogarithmicSpiralField
- * Fixed r term in source distribution for SNR and Pulsar
- * Fixed wrong mass inheritance for secondaries other than nuclei or electron/positron
+* Fixed sign for exponential decay of magn. field strength with Galactic height in LogarithmicSpiralField
+* Fixed r term in source distribution for SNR and Pulsar
+* Fixed wrong mass inheritance for secondaries other than nuclei or electron/positron
+* Fixed wrong generation of interval ranges in ObserverTimeEvolution
 
 ### New features:
- * Added new backwards-compatible function particleMass that returns particle mass also for non-nuclei
- * Added the new Galactic magnetic field models from Unger&Farrar arXiv:2311.12120
- * Added EBL model from Finke et al. 2022
+* Added new backwards-compatible function particleMass that returns particle mass also for non-nuclei
+* Added the new Galactic magnetic field models from Unger&Farrar arXiv:2311.12120
+* Added EBL model from Finke et al. 2022
+* Added overridable getTime function to ObserverTimeEvolution which is called instead of detList. 
+  Old functionalities are preserved by adding old functions and setter/getter functions.
 
 ### Interface changes:
 

doc/pages/example_notebooks/Diffusion/MomentumDiffusion.ipynb

@@ -240,8 +240,24 @@ 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 nIntervals;  // number of time invervals
+	bool isLogarithmicScaling = false;  // enables or disables logarithmic scaling for the intervals
+	bool doDetListConstruction = true;  // enables the construction of detList in the relevant functions (addTime, addTimeRange)
+	double minimum;  // the minimum time
+	double maximum;  // the maximum time
+	/** 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 +266,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 setMinimum(double min);
+	void setMaximum(double max){this->maximum = max;}
+	void setNIntervals(int numb){this->nIntervals = numb;}
+	void setIsLogarithmicScaling(bool log){this->isLogarithmicScaling = log;}
+
+	// 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 getMinimum() const {return minimum;}
+	double getMaximum() const {return maximum;}
+	int getNIntervals() const {return nIntervals;}
+	bool getIsLogarithmicScaling() const {return isLogarithmicScaling;}
+	/** 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;
 };
 

include/crpropa/module/Observer.h

@@ -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);
+	setMaximum(min + (numb - 1) * dist);
+	setMinimum(min);
+	setNIntervals(numb);
+	setIsLogarithmicScaling(false);
 }
 
 ObserverTimeEvolution::ObserverTimeEvolution(double min, double max, double numb, bool log) {
-	addTimeRange(min, max, numb, log);
+	setMinimum(min);
+	setMaximum(max);
+	setNIntervals(numb);
+	setIsLogarithmicScaling(log);
 }
 
+ObserverTimeEvolution::ObserverTimeEvolution(const std::vector<double> &detList){
+	setTimes(detList);
+}
 
 DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
 
-	if (detList.size()) {
+	if (nIntervals) {
 		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 >= nIntervals) {
 			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 < nIntervals-2) {
+				c->limitNextStep(getTime(index+1)-length);
 			}
 			c->setProperty(DI, Variant::fromUInt64(index+1));
 
@@ -299,29 +305,86 @@ 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();
+	detList.resize(0);
+	setNIntervals(0);
+}
+
+void ObserverTimeEvolution::constructDetListIfEmpty(){
+	if (detList.empty() && doDetListConstruction) {
+		std::vector<double> detListTemp;
+		size_t counter = 0;
+		while (getTime(counter)<=maximum) {
+			detListTemp.push_back(getTime(counter));
+			counter++;
+		}
+		detList.assign(detListTemp.begin(), detListTemp.end());
+	}
+}
+
+void ObserverTimeEvolution::addTime(const double &time){
+	constructDetListIfEmpty();
+	detList.push_back(time);
+	setNIntervals(nIntervals + 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
+	setNIntervals(detList.size());
+}
+
+void ObserverTimeEvolution::setTimes(const std::vector<double> &detList){
+	this->detList.assign(detList.begin(), detList.end());
+	setNIntervals(detList.size());
+	setMinimum(detList.front());
+	setMaximum(detList.back());
+	doDetListConstruction = false;
+}
+
+void ObserverTimeEvolution::setMinimum(double min){
+	if ( min <= 0 ){
+		std::cout << "minimum can not be <= 0 if isLogarithmicScaling=true" << std::endl;
+		throw new std::invalid_argument("minimum can not be <= 0 if isLogarithmicScaling=true");
+	}
+	this->minimum = min;
+}
+
+double ObserverTimeEvolution::getTime(size_t index) const {
+	if (!detList.empty()) {
+		return detList.at(index);
+	} else if (isLogarithmicScaling) {
+		return minimum * pow(maximum / minimum, index / (nIntervals - 1.0));
+	} else {
+		return minimum + index * (maximum - minimum) / (nIntervals - 1.0);
+	}
 }
 
 const std::vector<double>& ObserverTimeEvolution::getTimes() const {
-	return detList;
+	tempDetList.resize(nIntervals);
+	for (size_t i = 0; i < nIntervals; i++){
+		tempDetList[i] = 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 < nIntervals; i++)
+	  s << "  - " << getTime(i) / kpc;
 	return s.str();
 }