Merge pull request #20 from gcowan/rapidpid

Rapidpid

.gitlab-ci.yml

@@ -83,3 +83,15 @@ run-D02Kpi:
     artifacts:
         paths:
             - build/D02Kpi_hists.root
+
+run-B2DplusD0:
+    script:
+        - export RAPIDSIM_ROOT=$PWD
+        - cd build
+        - src/RapidSim.exe ../validation/B2DplusD0 10000
+        - mkdir plots
+        - root -b -q -l "../validation/compareHistograms.C(\"B2DplusD0\")"
+    stage: test
+    artifacts:
+        paths:
+            - build/B2DplusD0_hists.root

CMakeLists.txt

@@ -9,7 +9,7 @@ include("GNUInstallDirs")
 
 
 SET( ${PROJECT_NAME}_MAJOR_VERSION 1 )
-SET( ${PROJECT_NAME}_MINOR_VERSION 3 )
+SET( ${PROJECT_NAME}_MINOR_VERSION 4 )
 SET( ${PROJECT_NAME}_PATCH_LEVEL 0 )
 
 SET( CMAKE_BUILD_TYPE Release )

README.md

@@ -189,19 +189,27 @@ Particle settings should be defined after the corresponding `@#` tag using the s
 
 * `useEvtGen` :
   * Perform decays using the external EvtGen generator
-  * Syntax is `useEvtGen : TRUE
+  * Syntax is `useEvtGen : TRUE`
   * Note any value for this parameter will turn EvtGen ON (even FALSE)
   * To turn EvtGen OFF do not define this parameter
 
+* `pid` :
+  * Sets the particle identification scheme to use
+  * Syntax is `pid : <scheme>`, where 
+    * `<scheme>` is the name of the file that defines the scheme (default `LHCbGenericPID`)
+  * More types may be defined in $RAPIDSIM_ROOT/config/pid or $RAPIDSIM_CONFIG/config/pid
+
 ### Particle settings
 
 * `name`:
   * A user-friendly name for this particle to be used in variable names
   * Default: Automatically generated unique name of particle
 
 * `smear`:
-  * The type of momentum smearing to apply to this particle
-  * Supported: `LHCbGeneric`, `LHCbElectron`, `AtlasMuon`, or `AtlasHadron`
+  * The type of momentum or IP smearing to apply to this particle
+    * May be used multiple times for the same particle to define different types of smearing
+  * Supported momentum types: `LHCbGeneric`, `LHCbElectron`, `AtlasMuon`, or `AtlasHadron`
+  * Supported IP types: `LHCbGenericIP`
     * More types may be defined in $RAPIDSIM_ROOT/config/smear or $RAPIDSIM_CONFIG/config/smear
   * Default: `LHCbElectron` (for electrons/positrons), otherwise `LHCbGeneric`
 
@@ -238,6 +246,8 @@ Particle settings should be defined after the corresponding `@#` tag using the s
 * `PY`: The Y momentum of the combination of the given particles
 * `PZ`: The Z momentum of the combination of the given particles
 * `PT`: The transverse momentum of the combination of the given particles
+* `IP` : Impact parameter to own primary vertex
+* `SIGMAIP` : Error on impact parameter to own primary vertex
 * `eta`: The pseudorapidity of the combination
 * `phi`: The azimuthal angle of the combination
 * `y`: The rapidity of the combination
@@ -250,6 +260,11 @@ Particle settings should be defined after the corresponding `@#` tag using the s
 * `costheta`: The cosine of theta
 * `Mcorr`: The corrected mass of the combination of the given particles
   correcting for any invisible particles
+* `ProbNNmu` : Particle ID feature for the muon hypothesis
+* `ProbNNe` : Particle ID feature for the electron hypothesis
+* `ProbNNpi` : Particle ID feature for the pion hypothesis
+* `ProbNNk` : Particle ID feature for the kaon hypothesis
+* `ProbNNp` : Particle ID feature for the proton hypothesis
 
 ## External Generators
 

config/pid/LHCbGenericPID

@@ -0,0 +1,4 @@
+LHCb_K_PID.root 321 ProbNNk ProbNNmu ProbNNpi ProbNNp
+LHCb_Pi_PID.root 211 ProbNNk ProbNNmu ProbNNpi ProbNNp
+LHCb_Mu_PID.root 13 ProbNNk ProbNNmu ProbNNpi ProbNNp
+LHCb_P_PID.root 2212 ProbNNk ProbNNmu ProbNNpi ProbNNp

src/RapidConfig.cc

@@ -33,6 +33,16 @@ RapidConfig::~RapidConfig() {
 		delete itr2->second;
 		ipSmearCategories_.erase(itr2++);
 	}
+	std::map<RapidParam::ParamType, std::map<unsigned int, TH3D*>*>::iterator itr3 = pidHists_.begin();
+	while (itr3 != pidHists_.end()) {
+		std::map<unsigned int, TH3D*>::iterator itr4 = (*itr3).second->begin();
+		while (itr4 != (*itr3).second->end()) {
+			delete itr4->second;
+			(*itr3).second->erase(itr4++);
+		}
+		delete itr3->second;
+		pidHists_.erase(itr3++);
+	}
 	while(!parts_.empty()) {
 		delete parts_[parts_.size()-1];
 		parts_.pop_back();
@@ -534,6 +544,15 @@ bool RapidConfig::configGlobal(TString command, TString value) {
 	} else if(command=="useEvtGen") {
 		external_ = new RapidExternalEvtGen();
 		std::cout << "INFO in RapidConfig::configGlobal : will use external EvtGen generator to decay particles." << std::endl;
+	} else if (command=="pid") {
+		std::cout << "INFO in RapidConfig::configGlobal : setting pid type to " << value << "." << std::endl;
+		int from(0);
+		TString histFile;
+		value.Tokenize(histFile,from," ");
+		histFile = histFile.Strip(TString::kBoth);
+
+		pidLoaded_ = loadPID(histFile);
+		if(!pidLoaded_) return false;
 	}
 
 	return true;
@@ -592,7 +611,7 @@ RapidParam* RapidConfig::loadParam(TString paramStr) {
 		return 0;
 	}
 
-	RapidParam* param = new RapidParam(name,type,partlist,truth);
+	RapidParam * param = new RapidParam(name, type, partlist, truth);
 	return param;
 }
 
@@ -831,6 +850,88 @@ void RapidConfig::setSmearing(unsigned int particle, TString category) {
 	return;
 }
 
+bool RapidConfig::loadPID(TString category) {
+	TString path;
+	std::ifstream fin;
+	bool found(false);
+
+	path = getenv("RAPIDSIM_CONFIG");
+	std::cout << "INFO in RapidConfig::loadPID " << std::endl;
+
+	if(path!="") {
+		fin.open(path+"/config/pid/"+category, std::ifstream::in);
+		if( fin.good()) {
+			std::cout << "INFO in RapidConfig::loadPID : found pid category " << category << " in RAPIDSIM_CONFIG." << std::endl;
+			std::cout << "                                    this version will be used." << std::endl;
+			found = true;
+		} else {
+			std::cout << "INFO in RapidConfig::loadPID : pid category " << category << " not found in RAPIDSIM_CONFIG." << std::endl;
+			std::cout << "                                    checking RAPIDSIM_ROOT." << std::endl;
+			fin.close();
+		}
+	}
+
+	if(!found) {
+		path = getenv("RAPIDSIM_ROOT");
+		fin.open(path+"/config/pid/"+category, std::ifstream::in);
+		if( ! fin.good()) {
+			std::cout << "WARNING in RapidConfig::loadPID : failed to load pid category" << category << std::endl;
+			fin.close();
+			return false;
+		}
+	}
+
+	TString line("");
+	TString buffer("");
+	while (fin.good()) {
+		line.ReadLine(fin);
+		int from(0);
+		bool fileLoaded(false);
+		bool idLoaded(false);
+		unsigned int id(0);
+		TFile* file = NULL;
+		while (line.Tokenize(buffer, from)) {
+			if (buffer.Contains(".root") && !fileLoaded) {
+				std::cout << "INFO in RapidConfig::loadPID : loading root file " << buffer << std::endl;
+				file = TFile::Open(path+"/rootfiles/pid/"+buffer);
+				fileLoaded = true;
+				if(!file) {
+					std::cout << "WARNING in RapidConfig::loadPID : failed to load root file " << buffer << std::endl;
+					fin.close();
+					return false;
+				}
+				continue;
+			}
+			if ( fileLoaded && !idLoaded) {
+				id = buffer.Atoi();
+				idLoaded = true;
+				continue;
+			}
+			if ( fileLoaded && idLoaded && buffer.Contains("Prob")) {
+				std::cout << "INFO in RapidConfig::loadPID : loading histogram " << buffer << std::endl;
+				TH3D * hist = dynamic_cast<TH3D*>(file->Get(buffer));
+				//hist->Print();
+				if(!hist) {
+					std::cout << "WARNING in RapidConfig::loadPID : failed to load histogram " << buffer << std::endl;
+				}
+				RapidParam::ParamType type = RapidParam::typeFromString(buffer);
+				if(pidHists_.find(type)==pidHists_.end()) {
+					pidHists_[type] = new std::map<unsigned int, TH3D*>();
+				}
+				pidHists_[type]->insert(std::pair<unsigned int,TH3D*>(id, hist));
+			}
+			if(pidHists_.empty()) {
+				std::cout << "WARNING in RapidConfig::loadPID : failed to load any histograms for PID category " << category << std::endl;
+				file->Close();
+				fin.close();
+				return false;
+			}
+		}
+	}
+	fin.close();
+	return true;
+}
+
 bool RapidConfig::loadAcceptRejectHist(TString histFile, TString histName, RapidParam* paramX, RapidParam* paramY) {
 	if(accRejHisto_) {
 		std::cout << "WARNING in RapidConfig::loadAcceptRejectHist : accept/reject histogram already set." << std::endl
@@ -963,20 +1064,25 @@ void RapidConfig::setupDefaultParams() {
 	while(TString(paramStrStable_).Tokenize(buffer,from," ")) {
 		buffer = buffer.Strip(TString::kBoth,',');
 		RapidParam::ParamType type = RapidParam::typeFromString(buffer);
+
+		if ( buffer.Contains("ProbNN") && !pidLoaded_) loadPID("LHCbGenericPID");
+
 		if(type==RapidParam::UNKNOWN) {
 			std::cout << "WARNING in RapidConfig::setDefaultParams : Unknown parameter type " << buffer << "ignored." << std::endl;
 			continue;
 		} else {
 			for(unsigned int i=0; i<parts_.size(); ++i) {
 				RapidParticle* part = parts_[i];
 				if(part->nDaughters() == 0) {
-					std::vector<RapidParticle*> partlist;
-					partlist.push_back(part);
-					RapidParam* param = new RapidParam("", type, partlist, false);
+					std::map<unsigned int, TH3D*>* pidHists=0;
+					if (pidHists_.find(type)!=pidHists_.end() && pidHists_[type] && !pidHists_[type]->empty()) {
+						pidHists = pidHists_[type];
+					}
+					RapidParam* param = new RapidParam("", type, part, false, pidHists);
 					param->name();
 					paramsStable_.push_back(param);
 					if ( param->canBeTrue() ) {
-						param = new RapidParam("", type, partlist, true);
+						param = new RapidParam("", type, part, true, pidHists);
 						param->name();
 						paramsStable_.push_back(param);
 					}
@@ -997,13 +1103,11 @@ void RapidConfig::setupDefaultParams() {
 			for(unsigned int i=0; i<parts_.size(); ++i) {
 				RapidParticle* part = parts_[i];
 				if(part->nDaughters() > 0) {
-					std::vector<RapidParticle*> partlist;
-					partlist.push_back(part);
-					RapidParam* param = new RapidParam("", type, partlist, false);
+					RapidParam* param = new RapidParam("", type, part, false);
 					param->name();
 					paramsDecaying_.push_back(param);
 					if ( param->canBeTrue() ) {
-						param = new RapidParam("", type, partlist, true);
+						param = new RapidParam("", type, part, true);
 						param->name();
 						paramsDecaying_.push_back(param);
 					}

src/RapidConfig.h

@@ -6,8 +6,10 @@
 
 #include "TH1.h"
 #include "TH2.h"
+#include "TH3.h"
 #include "TString.h"
 #include "RapidAcceptance.h"
+#include "RapidParam.h"
 
 class RapidCut;
 class RapidDecay;
@@ -54,6 +56,7 @@ class RapidConfig {
 
 		void setSmearing(unsigned int particle, TString category);
 		bool loadSmearing(TString category);
+		bool loadPID(TString category);
 
 		bool loadAcceptRejectHist(TString histFile, TString histName, RapidParam* paramX, RapidParam* paramY);
 		bool loadParentKinematics();
@@ -98,6 +101,10 @@ class RapidConfig {
 		RapidParam* accRejParameterX_;
 		RapidParam* accRejParameterY_;
 
+		// PID histogram file
+		bool pidLoaded_;
+		std::map<RapidParam::ParamType, std::map<unsigned int, TH3D*>*> pidHists_;
+
 		//type of geometric acceptance to apply
 		RapidAcceptance::AcceptanceType acceptanceType_;