From ff1d205ab0aebccf6f80d135075e5016e8252aa5 Mon Sep 17 00:00:00 2001
From: Adam Davis <adamdddave@googlemail.com>
Date: Thu, 30 Aug 2018 17:10:43 +0200
Subject: [PATCH 1/2] Add Vertex Smearing capabilities
---
src/RapidDecay.cc | 2 +-
src/RapidParticle.h | 10 ++-
src/RapidVertex.cc | 12 +++-
src/RapidVertex.h | 5 +-
src/RapidVertexSmear.h | 14 ++++
src/RapidVertexSmearHisto.cc | 129 +++++++++++++++++++++++++++++++++++
src/RapidVertexSmearHisto.h | 40 +++++++++++
7 files changed, 205 insertions(+), 7 deletions(-)
create mode 100644 src/RapidVertexSmear.h
create mode 100644 src/RapidVertexSmearHisto.cc
create mode 100644 src/RapidVertexSmearHisto.h
diff --git a/src/RapidDecay.cc b/src/RapidDecay.cc
index e72188a..3f75a5a 100644
--- a/src/RapidDecay.cc
+++ b/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;
diff --git a/src/RapidParticle.h b/src/RapidParticle.h
index 577d834..91b0244 100644
--- a/src/RapidParticle.h
+++ b/src/RapidParticle.h
@@ -12,7 +12,7 @@
class RapidMomentumSmear;
class RapidParticleData;
class RapidIPSmear;
-//class RapidVtxSmear;
+class RapidVerterxSmear;
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
diff --git a/src/RapidVertex.cc b/src/RapidVertex.cc
index de963e1..b4ee08d 100644
--- a/src/RapidVertex.cc
+++ b/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);
+}
diff --git a/src/RapidVertex.h b/src/RapidVertex.h
index 41f1659..f733e96 100644
--- a/src/RapidVertex.h
+++ b/src/RapidVertex.h
@@ -6,6 +6,8 @@
#include "Math/Point3D.h"
#include "Math/Vector3D.h"
+#include "RapidVertexSmear.h"
+
class RapidVertex {
public:
RapidVertex(double x, double y, double z)
@@ -13,10 +15,11 @@ class RapidVertex {
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_;
diff --git a/src/RapidVertexSmear.h b/src/RapidVertexSmear.h
new file mode 100644
index 0000000..2206d18
--- /dev/null
+++ b/src/RapidVertexSmear.h
@@ -0,0 +1,14 @@
+#ifndef RAPIDVERTEXSMEAR_H
+#define RAPIDVERTEXSMEAR_H 1
+
+#include "Math/Point3D.h"
+
+class RapidVertexSmear {
+ public:
+ virtual ~RapidVertexSmear() {}
+
+ virtual ROOT::Math::XYZPoint smearVertex(ROOT::Math::XYZPoint vtx, int direction)=0;
+};
+
+
+#endif
diff --git a/src/RapidVertexSmearHisto.cc b/src/RapidVertexSmearHisto.cc
new file mode 100644
index 0000000..8b2ceb3
--- /dev/null
+++ b/src/RapidVertexSmearHisto.cc
@@ -0,0 +1,129 @@
+#include "RapidVertexSmearHisto.h"
+
+#include <iostream>
+
+#include "TMath.h"
+#include "TRandom.h"
+
+RapidVertexSmearHisto::~RapidVertexSmearHisto() {
+ while(!histos1D_.empty()) {
+ delete histos1D_[histos1D_.size()-1];
+ histos1D_.pop_back();
+ }
+
+ while(!histos2D_.empty()) {
+ delete histos2D_[histos2D_.size()-1];
+ histos2D_.pop_back();
+ }
+
+ while(!histos3D_.empty()) {
+ delete histos3D_[histos3D_.size()-1];
+ histos3D_.pop_back();
+ }
+
+}
+
+
+ROOT::Math::XYZPoint RapidVertexSmearHisto::smearVertex(ROOT::Math::XYZPoint vtx, int direction){
+ double smearx{0}, smeary{0}, smearz{0};
+ unsigned int iHist=0;
+ if(histos1D_.size()!=0){
+ while(true){//select the right histogram
+ if( iHist==thresholds_.size()-1 ) break;
+ if( threshVal_ < thresholds_[iHist+1] ) break;
+ }
+ switch(direction){
+ case 0 : smearx = histos1D_.at(iHist)->GetRandom();break;
+ case 1 : smeary = histos1D_.at(iHist)->GetRandom();break;
+ case 2 : smeary = histos1D_.at(iHist)->GetRandom();break;
+ default :smearx = histos1D_.at(iHist)->GetRandom();break;
+ }
+
+ }
+ else if(histos2D_.size()!=0){
+ while(true){
+ if( iHist==thresholds_.size()-1 ) break;
+ if( threshVal_ < thresholds_[iHist+1] ) break;
+ }
+ switch(direction){
+ case 0: histos2D_.at(iHist)->GetRandom2(smearx,smeary);break;
+ case 1: histos2D_.at(iHist)->GetRandom2(smearx,smearz);break;
+ case 2: histos2D_.at(iHist)->GetRandom2(smeary,smearz);break;
+ default:histos2D_.at(iHist)->GetRandom2(smearx,smeary);break;
+ }
+ }
+ else if(histos3D_.size()!=0){
+ while(true){
+ if( iHist==thresholds_.size()-1 ) break;
+ if( threshVal_ < thresholds_[iHist+1] ) break;
+ }
+ histos3D_.at(iHist)->GetRandom3(smearx,smeary,smearz);
+ }
+ else{
+ std::cout<<"Non-configured smearing. Should never come here. BREAK!"<<std::endl;
+ assert (0);
+ }
+ return ROOT::Math::XYZPoint(vtx.X()+smearx,vtx.Y()+smeary,vtx.Z()+smearz);
+}
+
+
+void RapidVertexSmearHisto::init(std::vector<double> thresholds, std::vector<TH1*>histos){
+ if(thresholds.size() < histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH1 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess histograms ignored." << std::endl;
+
+ while(thresholds.size() < histos.size()) {
+ histos.pop_back();
+ }
+ } else if(thresholds.size() > histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH1 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess thresholds ignored." << std::endl;
+
+ while(thresholds.size() > histos.size()) {
+ thresholds.pop_back();
+ }
+ }
+
+ thresholds_ = thresholds;
+ histos1D_ = histos;
+}
+void RapidVertexSmearHisto::init(std::vector<double> thresholds, std::vector<TH2*>histos){
+ if(thresholds.size() < histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH2 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess histograms ignored." << std::endl;
+
+ while(thresholds.size() < histos.size()) {
+ histos.pop_back();
+ }
+ } else if(thresholds.size() > histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH2 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess thresholds ignored." << std::endl;
+
+ while(thresholds.size() > histos.size()) {
+ thresholds.pop_back();
+ }
+ }
+
+ thresholds_ = thresholds;
+ histos2D_ = histos;
+}
+void RapidVertexSmearHisto::init(std::vector<double> thresholds, std::vector<TH3*>histos){
+ if(thresholds.size() < histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH3 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess histograms ignored." << std::endl;
+
+ while(thresholds.size() < histos.size()) {
+ histos.pop_back();
+ }
+ } else if(thresholds.size() > histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH3 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess thresholds ignored." << std::endl;
+
+ while(thresholds.size() > histos.size()) {
+ thresholds.pop_back();
+ }
+ }
+
+ thresholds_ = thresholds;
+ histos3D_ = histos;
+}
diff --git a/src/RapidVertexSmearHisto.h b/src/RapidVertexSmearHisto.h
new file mode 100644
index 0000000..3b38f22
--- /dev/null
+++ b/src/RapidVertexSmearHisto.h
@@ -0,0 +1,40 @@
+#ifndef RAPIDVERTEXSMEARHISTO_H
+#define RAPIDVERTEXSMEARHISTO_H 1
+
+#include "TH1.h"
+#include "TH2.h"
+#include "TH3.h"
+
+#include "RapidVertexSmear.h"
+
+
+class RapidVertexSmearHisto : public RapidVertexSmear {
+ ///Simple class to smear Vertex position by a TH1, TH2 or TH3.
+ //variable 'direction' controls which direction to smear
+ // TH1: 0 == x, 1== y, 2== z
+ // TH2: 0 == xy, 2== xz, 3== yz
+ // TH3: only 0==xyz
+ //thresholds allow for possible extra variable dependence, e.g. decay time.
+
+ public:
+ RapidVertexSmearHisto(std::vector<double> thresholds, std::vector<TH1*> histos){ init(thresholds, histos); };
+ RapidVertexSmearHisto(std::vector<double> thresholds, std::vector<TH2*> histos){ init(thresholds, histos); };
+ RapidVertexSmearHisto(std::vector<double> thresholds, std::vector<TH3*> histos){ init(thresholds, histos); };
+ ~RapidVertexSmearHisto();
+
+ ROOT::Math::XYZPoint smearVertex(ROOT::Math::XYZPoint vtx, int direction);
+ void SetThresholdVariableValue(double val){threshVal_=val;}
+ private:
+ void init(std::vector<double> thresholds, std::vector<TH1*>histos);
+ void init(std::vector<double> thresholds, std::vector<TH2*>histos);
+ void init(std::vector<double> thresholds, std::vector<TH3*>histos);
+ std::vector<TH1*> histos1D_;
+ std::vector<TH2*> histos2D_;
+ std::vector<TH3*> histos3D_;
+ std::vector<double>thresholds_;
+ double threshVal_;
+};
+
+
+
+#endif
From ee21eedf9eff2d7dcd5c1104c25498670693e454 Mon Sep 17 00:00:00 2001
From: Adam Davis <adamdddave@googlemail.com>
Date: Thu, 30 Aug 2018 23:54:37 +0200
Subject: [PATCH 2/2] Add support for secondary vertex smearing Implemented
through base class RapidVertexSmear Extended to Histogram (1,2 and 3D)
through RapidVertexSmearHisto Configurable on decaying particles implemented
as ``` @0 name : D0_0 endVertexSmear : MyVertexSmear ``` where
`MyVertexSmear` points to a configuration file in
`$RAPIDSIM_CONFIG/config/smear` or `$RAPIDSIM_ROOT/config/smear`.
This file has the structure
```
Dstar_vertex_resolutions.root
VTX3D
0. hdst3d
```
where the first line is the rootfile in `$RAPIDSIM_{CONFIG/ROOT}/rootfiles/smear`
second line is the type of smearing (here a TH3 to be randomly sampled)
and 3rd line is first a threshold and second the histogram for that threshold.
todo:
(1) add support for threshold variable choice in configuration
(2) add support for smearind direction in config file
---
README.md | 15 +++
src/RapidConfig.cc | 207 +++++++++++++++++++++++++++++++++++
src/RapidConfig.h | 11 +-
src/RapidParticle.h | 2 +-
src/RapidVertexSmearHisto.cc | 186 +++++++++++++++----------------
5 files changed, 324 insertions(+), 97 deletions(-)
diff --git a/README.md b/README.md
index df2a75f..dea8d7f 100644
--- a/README.md
+++ b/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
diff --git a/src/RapidConfig.cc b/src/RapidConfig.cc
index 7f02d2c..36a8850 100644
--- a/src/RapidConfig.cc
+++ b/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;
+}
diff --git a/src/RapidConfig.h b/src/RapidConfig.h
index 77b4002..9860c7c 100644
--- a/src/RapidConfig.h
+++ b/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,6 +59,10 @@ 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();
@@ -66,9 +70,10 @@ class RapidConfig {
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_;
diff --git a/src/RapidParticle.h b/src/RapidParticle.h
index 91b0244..eb38c5a 100644
--- a/src/RapidParticle.h
+++ b/src/RapidParticle.h
@@ -12,7 +12,7 @@
class RapidMomentumSmear;
class RapidParticleData;
class RapidIPSmear;
-class RapidVerterxSmear;
+class RapidVertexSmear;
class RapidParticle {
public:
diff --git a/src/RapidVertexSmearHisto.cc b/src/RapidVertexSmearHisto.cc
index 8b2ceb3..a110cdc 100644
--- a/src/RapidVertexSmearHisto.cc
+++ b/src/RapidVertexSmearHisto.cc
@@ -6,124 +6,124 @@
#include "TRandom.h"
RapidVertexSmearHisto::~RapidVertexSmearHisto() {
- while(!histos1D_.empty()) {
- delete histos1D_[histos1D_.size()-1];
- histos1D_.pop_back();
- }
+ while(!histos1D_.empty()) {
+ delete histos1D_[histos1D_.size()-1];
+ histos1D_.pop_back();
+ }
- while(!histos2D_.empty()) {
- delete histos2D_[histos2D_.size()-1];
- histos2D_.pop_back();
- }
+ while(!histos2D_.empty()) {
+ delete histos2D_[histos2D_.size()-1];
+ histos2D_.pop_back();
+ }
- while(!histos3D_.empty()) {
- delete histos3D_[histos3D_.size()-1];
- histos3D_.pop_back();
- }
+ while(!histos3D_.empty()) {
+ delete histos3D_[histos3D_.size()-1];
+ histos3D_.pop_back();
+ }
}
ROOT::Math::XYZPoint RapidVertexSmearHisto::smearVertex(ROOT::Math::XYZPoint vtx, int direction){
- double smearx{0}, smeary{0}, smearz{0};
- unsigned int iHist=0;
- if(histos1D_.size()!=0){
- while(true){//select the right histogram
- if( iHist==thresholds_.size()-1 ) break;
- if( threshVal_ < thresholds_[iHist+1] ) break;
- }
- switch(direction){
- case 0 : smearx = histos1D_.at(iHist)->GetRandom();break;
- case 1 : smeary = histos1D_.at(iHist)->GetRandom();break;
- case 2 : smeary = histos1D_.at(iHist)->GetRandom();break;
- default :smearx = histos1D_.at(iHist)->GetRandom();break;
- }
+ double smearx{0}, smeary{0}, smearz{0};
+ unsigned int iHist=0;
+ if(histos1D_.size()!=0){
+ while(true){//select the right histogram
+ if( iHist==thresholds_.size()-1 ) break;
+ if( threshVal_ < thresholds_[iHist+1] ) break;
+ }
+ switch(direction){
+ case 0 : smearx = histos1D_.at(iHist)->GetRandom();break;
+ case 1 : smeary = histos1D_.at(iHist)->GetRandom();break;
+ case 2 : smeary = histos1D_.at(iHist)->GetRandom();break;
+ default :smearx = histos1D_.at(iHist)->GetRandom();break;
+ }
- }
- else if(histos2D_.size()!=0){
- while(true){
- if( iHist==thresholds_.size()-1 ) break;
- if( threshVal_ < thresholds_[iHist+1] ) break;
- }
- switch(direction){
- case 0: histos2D_.at(iHist)->GetRandom2(smearx,smeary);break;
- case 1: histos2D_.at(iHist)->GetRandom2(smearx,smearz);break;
- case 2: histos2D_.at(iHist)->GetRandom2(smeary,smearz);break;
- default:histos2D_.at(iHist)->GetRandom2(smearx,smeary);break;
- }
- }
- else if(histos3D_.size()!=0){
- while(true){
- if( iHist==thresholds_.size()-1 ) break;
- if( threshVal_ < thresholds_[iHist+1] ) break;
- }
- histos3D_.at(iHist)->GetRandom3(smearx,smeary,smearz);
- }
- else{
- std::cout<<"Non-configured smearing. Should never come here. BREAK!"<<std::endl;
- assert (0);
- }
- return ROOT::Math::XYZPoint(vtx.X()+smearx,vtx.Y()+smeary,vtx.Z()+smearz);
+ }
+ else if(histos2D_.size()!=0){
+ while(true){
+ if( iHist==thresholds_.size()-1 ) break;
+ if( threshVal_ < thresholds_[iHist+1] ) break;
+ }
+ switch(direction){
+ case 0: histos2D_.at(iHist)->GetRandom2(smearx,smeary);break;
+ case 1: histos2D_.at(iHist)->GetRandom2(smearx,smearz);break;
+ case 2: histos2D_.at(iHist)->GetRandom2(smeary,smearz);break;
+ default:histos2D_.at(iHist)->GetRandom2(smearx,smeary);break;
+ }
+ }
+ else if(histos3D_.size()!=0){
+ while(true){
+ if( iHist==thresholds_.size()-1 ) break;
+ if( threshVal_ < thresholds_[iHist+1] ) break;
+ }
+ histos3D_.at(iHist)->GetRandom3(smearx,smeary,smearz);
+ }
+ else{
+ std::cout<<"Non-configured smearing. Should never come here. BREAK!"<<std::endl;
+ assert (0);
+ }
+ return ROOT::Math::XYZPoint(vtx.X()+smearx,vtx.Y()+smeary,vtx.Z()+smearz);
}
void RapidVertexSmearHisto::init(std::vector<double> thresholds, std::vector<TH1*>histos){
- if(thresholds.size() < histos.size()) {
- std::cout << "WARNING in RapidVertexSmearHisto::init : TH1 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
- std::cout << " excess histograms ignored." << std::endl;
+ if(thresholds.size() < histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH1 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess histograms ignored." << std::endl;
- while(thresholds.size() < histos.size()) {
- histos.pop_back();
- }
- } else if(thresholds.size() > histos.size()) {
- std::cout << "WARNING in RapidVertexSmearHisto::init : TH1 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
- std::cout << " excess thresholds ignored." << std::endl;
+ while(thresholds.size() < histos.size()) {
+ histos.pop_back();
+ }
+ } else if(thresholds.size() > histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH1 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess thresholds ignored." << std::endl;
- while(thresholds.size() > histos.size()) {
- thresholds.pop_back();
- }
- }
+ while(thresholds.size() > histos.size()) {
+ thresholds.pop_back();
+ }
+ }
- thresholds_ = thresholds;
- histos1D_ = histos;
+ thresholds_ = thresholds;
+ histos1D_ = histos;
}
void RapidVertexSmearHisto::init(std::vector<double> thresholds, std::vector<TH2*>histos){
- if(thresholds.size() < histos.size()) {
- std::cout << "WARNING in RapidVertexSmearHisto::init : TH2 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
- std::cout << " excess histograms ignored." << std::endl;
+ if(thresholds.size() < histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH2 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess histograms ignored." << std::endl;
- while(thresholds.size() < histos.size()) {
- histos.pop_back();
- }
- } else if(thresholds.size() > histos.size()) {
- std::cout << "WARNING in RapidVertexSmearHisto::init : TH2 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
- std::cout << " excess thresholds ignored." << std::endl;
+ while(thresholds.size() < histos.size()) {
+ histos.pop_back();
+ }
+ } else if(thresholds.size() > histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH2 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess thresholds ignored." << std::endl;
- while(thresholds.size() > histos.size()) {
- thresholds.pop_back();
- }
- }
+ while(thresholds.size() > histos.size()) {
+ thresholds.pop_back();
+ }
+ }
- thresholds_ = thresholds;
+ thresholds_ = thresholds;
histos2D_ = histos;
}
void RapidVertexSmearHisto::init(std::vector<double> thresholds, std::vector<TH3*>histos){
- if(thresholds.size() < histos.size()) {
- std::cout << "WARNING in RapidVertexSmearHisto::init : TH3 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
- std::cout << " excess histograms ignored." << std::endl;
+ if(thresholds.size() < histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH3 : too many histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess histograms ignored." << std::endl;
- while(thresholds.size() < histos.size()) {
+ while(thresholds.size() < histos.size()) {
histos.pop_back();
- }
- } else if(thresholds.size() > histos.size()) {
- std::cout << "WARNING in RapidVertexSmearHisto::init : TH3 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
- std::cout << " excess thresholds ignored." << std::endl;
+ }
+ } else if(thresholds.size() > histos.size()) {
+ std::cout << "WARNING in RapidVertexSmearHisto::init : TH3 : too few histograms provided. Number of histograms should match number of thresholds." << std::endl;
+ std::cout << " excess thresholds ignored." << std::endl;
- while(thresholds.size() > histos.size()) {
- thresholds.pop_back();
- }
- }
+ while(thresholds.size() > histos.size()) {
+ thresholds.pop_back();
+ }
+ }
- thresholds_ = thresholds;
- histos3D_ = histos;
+ thresholds_ = thresholds;
+ histos3D_ = histos;
}