WIP: Add support for secondary vertex smearing

README.md

@@ -255,6 +255,21 @@ Particle settings should be defined after the corresponding `@#` tag using the s
     parameters to be passed to the model
   * Default: PHSP
 
+* `endVertexSmear` :
+  * Sets the smearing of the end vertex, depending on a configuration in
+    `$RAPIDSIM_ROOT/config/smear` or `$RAPIDSIM_CONFIG/config/smear`
+  * Current implementation is for smearing by histograms, implemented in `RapidVertexSmearHisto`
+  * Syntax is `endVertexSmear : MySmear`
+    where `MySmear` is a configuration file defined as above. The file will look like
+    ```
+    rootfile.root
+    VTX3D
+    0. histname
+    ```
+    where `rootfile.root` is the rootfile in `$RAPIDSIM_CONFIG/rootfiles/smear` or `RAPIDSIM_ROOT/rootfiles/smear`
+    `VTX3D` is the smearing type, here a `TH3*` histogram
+    `0.` is the first threshold
+    `histname` is the histogram name in `rootfile.root` for this threshold.
 ## Parameters
 
 * `M`: The invariant mass of the combination of the given particles

src/RapidConfig.cc

@@ -23,6 +23,7 @@
 #include "RapidParticle.h"
 #include "RapidParticleData.h"
 #include "RapidPID.h"
+#include "RapidVertexSmearHisto.h"
 
 RapidConfig::~RapidConfig() {
 	std::map<TString, RapidMomentumSmear*>::iterator itr = momSmearCategories_.begin();
@@ -40,6 +41,11 @@ RapidConfig::~RapidConfig() {
 		delete itr3->second;
 		pidHists_.erase(itr3++);
 	}
+	std::map<TString, RapidVertexSmear*>::iterator itr4 = vtxSmearCategories_.begin();
+	while(itr4!=vtxSmearCategories_.end()) {
+	        delete itr4->second;
+		vtxSmearCategories_.erase(itr4++);
+	}
 	while(!parts_.empty()) {
 		delete parts_[parts_.size()-1];
 		parts_.pop_back();
@@ -429,6 +435,8 @@ bool RapidConfig::configParticle(unsigned int part, TString command, TString val
 		parts_[part]->setEvtGenDecayModel(value);
 		std::cout << "INFO in RapidConfig::configParticle : set EvtGen decay model for particle " << parts_[part]->name() << std::endl
 			  << "                                    : " << value << std::endl;
+	} else if(command=="endVertexSmear"){
+	       setVertexSmearing(part, value);
 	}
 
 	return true;
@@ -1296,3 +1304,202 @@ TH1* RapidConfig::reduceHistogram(TH1* histo, double min, double max) {
 
 	return histoOut;
 }
+
+
+bool RapidConfig::loadVertexSmearing(TString category){
+	//this "category" is actually the name of a file in /rootfiles/smear/category
+	//the category is itself a textfile including the format
+	// rootfname 
+	TString path;
+	std::ifstream fin;
+	bool found(false);
+	found = true;
+	path = getenv("RAPIDSIM_CONFIG");
+
+	if(path!=""){
+		fin.open(path+"/config/smear/"+category,std::ifstream::in);
+		if( fin.good() ){
+			std::cout << "INFO in RapidConfig::loadVertexSmearing : found smearing category " << category << "in RAPIDSIM_CONFIG." << std::endl
+				  << "                                          this version will be used." <<std::endl;
+		} else {
+			std::cout << "INFO in RapidConfig::loadVertexSmearing : smearing category " << category << " not found in RAPIDSIM_CONFIG." << std::endl
+				  << "                                          checking RAPIDSIM_ROOT. "<< std::endl;
+			fin.close();
+		}
+	}else {
+		std::cout <<" INFO in RapidConfig::loadVertexSmearing : $RAPIDSIM_CONFIG not set. trying RAPIDSIM_ROOT." << std::endl;
+		found = false;
+	}
+	if(!found){
+
+		path = getenv("RAPIDSIM_ROOT");
+
+		fin.open(path+"/config/smear/"+category,std::ifstream::in);
+		if( ! fin.good() ){
+			std::cout << "WARNING in RapidConfig::loadVertexSmearing : failed to load smearing from $RAPIDSIM_ROOT for category " << category << std::endl;
+		}
+		else{
+			std::cout<< "INFO in RapidConfig::loadVertexSmearing : found file! continuing " <<std::endl;
+		}
+	}
+
+
+	TString filename("");
+	filename.ReadToken(fin);
+
+	TFile* file = NULL;
+	if( filename !="NULL" ) {
+		file  =  TFile::Open(path+"/rootfiles/smear/"+filename);
+		if(!file) {
+			std::cout << "WARNING in RapidConfig::loadVertexSmearing : failed to load root file " << filename << std::endl;
+			fin.close();
+			return false;
+		}
+	}
+	// format from reading: VTXND thres1 hist1
+	TString type("");
+	type.ReadToken(fin);
+	if(type=="VTX1D") { //todo. add direction as third fin quantity
+		TString histname("");
+		std::vector<TH1*> hists;
+		std::vector<double> thresholds;
+		double threshold(0.);
+		while( true ){
+			fin >> threshold;
+			histname.ReadToken(fin);
+
+			//add check if user doesn't put a return before the end of the configuration.
+			if(""==histname){
+				std::cout<<" INFO in RapidConfig::loadVertexSmearing : out of histograms !"<<std::endl;
+				break;
+			}
+			if(fin.good()) {
+				TH1* hist = dynamic_cast<TH1*>(file->Get(histname));
+				if(!hist || !check1D(hist)) {
+					std::cout << "WARNING in RapidConfig::loadVertexSmearing : failed to load histgram " << histname << std::endl
+						  << "                                             threshold will be ignored." << std::endl;
+				} else {
+					hists.push_back(hist);
+					thresholds.push_back(threshold);
+					break;
+				}	
+			}      
+		}
+		if(hists.size() == 0) {
+			std::cout << "WARNING in RapidConfig::loadVertexSmearing : failed to load any histograms for category " << category << std::endl;
+			file->Close();
+			fin.close();
+			return false;
+		}
+		vtxSmearCategories_[category] = new RapidVertexSmearHisto(thresholds,hists);
+	} else if(type=="VTX2D") {
+		TString histname("");
+		std::vector<TH2*> hists;
+		std::vector<double> thresholds;
+		double threshold(0.);
+		while( true ){
+			fin >> threshold;
+			histname.ReadToken(fin);
+
+			//add check if user doesn't put a return before the end of the configuration.
+			if(""==histname){
+				std::cout<<" INFO in RapidConfig::loadVertexSmearing : out of histograms !"<<std::endl;
+				break;
+			}
+			if(fin.good()) {
+				TH2* hist = dynamic_cast<TH2*>(file->Get(histname));
+				if(!hist || !check2D(hist)) {
+					std::cout << "WARNING in RapidConfig::loadVertexSmearing : failed to load histgram " << histname << std::endl
+						  << "                                             threshold will be ignored." << std::endl;
+				} else {
+					hists.push_back(hist);
+					thresholds.push_back(threshold);
+					break;
+				}	
+			}      
+		}
+		if(hists.size() == 0) {
+			std::cout << "WARNING in RapidConfig::loadVertexSmearing : failed to load any histograms for category " << category << std::endl;
+			file->Close();
+			fin.close();
+			return false;
+		}
+		vtxSmearCategories_[category] = new RapidVertexSmearHisto(thresholds,hists);
+	} else if(type=="VTX3D") {
+		TString histname("");
+		std::vector<TH3*> hists;
+		std::vector<double> thresholds;
+		double threshold(0.);
+		while( true ){
+			fin >> threshold;
+			histname.ReadToken(fin);
+
+			//add check if user doesn't put a return before the end of the configuration.
+			if(""==histname){
+				std::cout<<" INFO in RapidConfig::loadVertexSmearing : out of histograms !"<<std::endl;
+				break;
+			}
+			if(fin.good()) {
+				TH3* hist = dynamic_cast<TH3*>(file->Get(histname));
+				if(!hist || !check3D(hist)) {
+					std::cout << "WARNING in RapidConfig::loadVertexSmearing : failed to load histgram " << histname << std::endl
+						  << "                                             threshold will be ignored." << std::endl;
+				} else {
+					hists.push_back(hist);
+					thresholds.push_back(threshold);
+					break;
+				}	
+			}      
+		}
+		if(hists.size() == 0) {
+			std::cout << "WARNING in RapidConfig::loadVertexSmearing : failed to load any histograms for category " << category << std::endl;
+			file->Close();
+			fin.close();
+			return false;
+		}
+		vtxSmearCategories_[category] = new RapidVertexSmearHisto(thresholds,hists);
+	} 
+	fin.close();
+	return true;
+}
+
+void RapidConfig::setVertexSmearing(unsigned int particle, TString category){
+	//expect that category is a space separated list of type thresh hist thresh hist...
+	TString buffer ;
+	int from = 0;
+	std::vector<TString> things;
+	while(category.Tokenize(buffer,from, " ")){
+		things.push_back(buffer);
+	}
+
+	for(auto thing : things){
+		std::cout<<"got thing "<<thing<<std::endl;
+	}
+	std::cout << "in setVertexSmearing : particle =  "<<particle<<", category  = "<<category<<std::endl;
+	if(particle>parts_.size()){
+		std::cout << "WARNING in RapidConfig::setVertexSmearing : particle " << particle << " does not exist - smearing functions not set." <<std::endl;
+		return;
+	}
+	if( parts_[particle]->nDaughters() ==0 ){
+		std::cout << "WARNING in RapidConfig::setVertexSmearing : particle " << particle << " is not composite - cannot set end vertex of non-composite particle" << std::endl;
+		return;
+	}
+
+	std::cout<<" INFO in RapidConfig::setVertexSmearing : setting vertex smearing function for particle " << particle << "(category : " << category << ")"<<std::endl;
+	bool loadedsmear = true;
+	if(!vtxSmearCategories_.count(category)) {
+		if(!loadVertexSmearing(category)){
+			std::cout << "WARNING in RapidConfig::setVertexSmearing : failed to load end vertex smearing for category " << category << "." <<std::endl
+				  << "                                            End vertex smearing functions not set for particle " << particle << "." << std::endl;
+			loadedsmear = false;
+		}
+		else if(!vtxSmearCategories_.count(category)) {
+			loadedsmear = false;
+		}
+		if(loadedsmear) {
+			parts_[particle]->setVertexSmear(vtxSmearCategories_[category]);
+		}
+	}
+
+	return;
+}

src/RapidConfig.h

@@ -17,7 +17,7 @@ class RapidExternalGenerator;
 class RapidHistWriter;
 class RapidMomentumSmear;
 class RapidIPSmear;
-class RapidVtxSmear;
+class RapidVertexSmear;
 class RapidParam;
 class RapidParticle;
 class RapidPID;
@@ -59,16 +59,21 @@ class RapidConfig {
 		bool loadSmearing(TString category);
 		bool loadPID(TString category);
 
+		//vertex smearing
+		void setVertexSmearing(unsigned int particle, TString category);
+		bool loadVertexSmearing(TString category);
+		//
 		bool loadAcceptRejectHist(TString histFile, TString histName, RapidParam* paramX, RapidParam* paramY);
 		bool loadParentKinematics();
 
 		void setupDefaultParams();
 
 		bool check1D(TH1* hist) { return (dynamic_cast<TH1F*>(hist) || dynamic_cast<TH1D*>(hist)); }
 		bool check2D(TH1* hist) { return (dynamic_cast<TH2F*>(hist) || dynamic_cast<TH2D*>(hist)); }
+		bool check3D(TH1* hist) { return (dynamic_cast<TH3F*>(hist) || dynamic_cast<TH3D*>(hist)); }
 
 		TH1* reduceHistogram(TH1* histo, double min, double max);
-
+		
 		TString fileName_;
 
 		std::vector<RapidParticle*> parts_;
@@ -95,7 +100,7 @@ class RapidConfig {
 		//IP smearing lookup for each smearing category
 		std::map<TString, RapidIPSmear*> ipSmearCategories_;
 		//Vtx smearing lookup for each smearing category, placeholder for now
-		//std::map<TString, RapidVtxSmear*> vtxSmearCategories_;
+		std::map<TString, RapidVertexSmear*> vtxSmearCategories_;
 
 		//accept reject hist to sculpt kinematics
 		TH1* accRejHisto_;

src/RapidDecay.cc

@@ -243,7 +243,7 @@ bool RapidDecay::genDecay(bool acceptAny) {
 				double dvx  = part->getOriginVertex()->getVertex(true).X() + part->getP().Vect().Unit().X()*dist;
 				double dvy  = part->getOriginVertex()->getVertex(true).Y() + part->getP().Vect().Unit().Y()*dist;
 				double dvz  = part->getOriginVertex()->getVertex(true).Z() + part->getP().Vect().Unit().Z()*dist;
-				part->getDecayVertex()->setXYZ(dvx,dvy,dvz);
+				part->getDecayVertex()->setXYZ(dvx,dvy,dvz,part->getDecayVertexSmear());
 			}
 
 			int j=0;

src/RapidParticle.h

@@ -12,7 +12,7 @@
 class RapidMomentumSmear;
 class RapidParticleData;
 class RapidIPSmear;
-//class RapidVtxSmear;
+class RapidVertexSmear;
 
 class RapidParticle {
 	public:
@@ -23,7 +23,8 @@ class RapidParticle {
 			  evtGenModel_("PHSP"),
 			  currentHypothesis_(0),
 			  originVertex_(0),
-			  decayVertex_(0)
+                          decayVertex_(0),
+                          decayVertexSmear_(0)
 			{setPtEtaPhi(0,0,0); setupVertices();}
 
 		~RapidParticle() {}
@@ -100,7 +101,9 @@ class RapidParticle {
 
 		TString evtGenDecayModel() { return evtGenModel_; }
 		void setEvtGenDecayModel(TString value) { evtGenModel_ = value; }
-
+		// added methods for verex smearing
+		void setVertexSmear(RapidVertexSmear* smear){decayVertexSmear_ = smear;}
+		RapidVertexSmear* getDecayVertexSmear(){return decayVertexSmear_;}
 	private:
 		bool hasFlavour(int flavour);
 
@@ -154,5 +157,6 @@ class RapidParticle {
 
 		RapidVertex * originVertex_;
 		RapidVertex * decayVertex_;
+		RapidVertexSmear* decayVertexSmear_;
 };
 #endif

src/RapidVertex.cc

@@ -10,9 +10,14 @@ ROOT::Math::XYZPoint RapidVertex::getVertex(bool truth) {
 	else return vertexSmeared_;
 }
 
-void RapidVertex::setXYZ(double x, double y, double z) {
+void RapidVertex::setXYZ(double x, double y, double z, RapidVertexSmear* smear) {
 	vertexTrue_ = ROOT::Math::XYZPoint(x,y,z);
-	smearVertex();
+	if(!smear){
+	  smearVertex();
+	}
+	else{
+	  smearVertex(smear);
+	}
 }
 
 void RapidVertex::smearVertex() {
@@ -28,3 +33,6 @@ void RapidVertex::smearVertex() {
 			vertexTrue_.Z() + gRandom->Gaus(0,zS)*1000.);
 }
 
+void RapidVertex::smearVertex(RapidVertexSmear* smear){
+  vertexSmeared_ = smear->smearVertex(vertexTrue_,0);
+}

src/RapidVertex.h

@@ -6,17 +6,20 @@
 #include "Math/Point3D.h"
 #include "Math/Vector3D.h"
 
+#include "RapidVertexSmear.h"
+
 class RapidVertex {
 	public:
 		RapidVertex(double x, double y, double z)
 			: nPVTracks_(5), vertexTrue_(x, y, z) {smearVertex();}
 
 		ROOT::Math::XYZPoint getVertex(bool truth);
 
-		void setXYZ(double x, double y, double z);
+		void setXYZ(double x, double y, double z, RapidVertexSmear* smear = nullptr);
 
 	private:
 		void smearVertex();
+		void smearVertex(RapidVertexSmear* smear);
 		unsigned int nPVTracks_;
 
 		ROOT::Math::XYZPoint vertexTrue_;