Fork me on GitHub

Context Navigation

source: git/modules/RunPUPPI.cc@ 21eab4f

Visit:

Last change on this file since 21eab4f was 659c7b6, checked in by Michele Selvaggi <michele.selvaggi@…>, 8 years ago
fix PUPPI includes
Property mode set to `100644`
File size: 11.3 KB

Rev	Line
[fa33983]	1	#include "modules/RunPUPPI.h"
	2
[e921a28]	3	#include "PUPPI/RecoObj2.hh"
	4	#include "PUPPI/AlgoObj.hh"
[659c7b6]	5	#include "PUPPI/PuppiContainer.hh"
	6
	7	#include "fastjet/PseudoJet.hh"
[fa33983]	8
	9	#include "classes/DelphesClasses.h"
	10	#include "classes/DelphesFactory.h"
	11	#include "classes/DelphesFormula.h"
	12
	13	#include <algorithm>
	14	#include <stdexcept>
	15	#include <iostream>
	16	#include <sstream>
	17	#include <vector>
	18
	19	using namespace std;
[659c7b6]	20	using namespace fastjet;
[fa33983]	21
	22	//------------------------------------------------------------------------------
	23	RunPUPPI::RunPUPPI() :
	24	fItTrackInputArray(0),
	25	fItNeutralInputArray(0)
	26	{}
	27
	28	//------------------------------------------------------------------------------
	29	RunPUPPI::~RunPUPPI(){}
	30
	31	//------------------------------------------------------------------------------
	32
	33	void RunPUPPI::Init(){
	34	// input collection
	35	fTrackInputArray = ImportArray(GetString("TrackInputArray", "Calorimeter/towers"));
	36	fItTrackInputArray = fTrackInputArray->MakeIterator();
	37	fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "Calorimeter/towers"));
	38	fItNeutralInputArray = fNeutralInputArray->MakeIterator();
	39	fPVInputArray = ImportArray(GetString("PVInputArray", "PV"));
	40	fPVItInputArray = fPVInputArray->MakeIterator();
[e921a28]	41	// puppi parameters
	42	fApplyNoLep = GetBool("UseNoLep", true);
[fa33983]	43	fMinPuppiWeight = GetDouble("MinPuppiWeight", 0.01);
	44	fUseExp = GetBool("UseExp", false);
[e921a28]	45	// read eta min ranges
[fa33983]	46	ExRootConfParam param = GetParam("EtaMinBin");
	47	fEtaMinBin.clear();
	48	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMinBin.push_back(param[iMap].GetDouble());
[e921a28]	49	// read eta max ranges
[fa33983]	50	param = GetParam("EtaMaxBin");
	51	fEtaMaxBin.clear();
	52	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fEtaMaxBin.push_back(param[iMap].GetDouble());
[e921a28]	53	// read pt min value
[fa33983]	54	param = GetParam("PtMinBin");
	55	fPtMinBin.clear();
	56	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fPtMinBin.push_back(param[iMap].GetDouble());
	57	// read cone size
	58	param = GetParam("ConeSizeBin");
	59	fConeSizeBin.clear();
	60	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fConeSizeBin.push_back(param[iMap].GetDouble());
	61	// read RMS min pt
	62	param = GetParam("RMSPtMinBin");
	63	fRMSPtMinBin.clear();
	64	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSPtMinBin.push_back(param[iMap].GetDouble());
[e921a28]	65	// read RMS scale factor
[fa33983]	66	param = GetParam("RMSScaleFactorBin");
	67	fRMSScaleFactorBin.clear();
	68	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fRMSScaleFactorBin.push_back(param[iMap].GetDouble());
[e921a28]	69	// read neutral pt min cut
[fa33983]	70	param = GetParam("NeutralMinEBin");
	71	fNeutralMinEBin.clear();
	72	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralMinEBin.push_back(param[iMap].GetDouble());
[e921a28]	73	// read neutral pt min slope
[fa33983]	74	param = GetParam("NeutralPtSlope");
	75	fNeutralPtSlope.clear();
	76	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fNeutralPtSlope.push_back(param[iMap].GetDouble());
	77	// read apply chs
[e921a28]	78	//param = GetParam("ApplyCHS");
	79	//fApplyCHS.clear();
	80	//for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyCHS.push_back(param[iMap].GetBool());
	81	// read use charged
[fa33983]	82	param = GetParam("UseCharged");
	83	fUseCharged.clear();
	84	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fUseCharged.push_back(param[iMap].GetBool());
[e921a28]	85	// read apply chs correction
[fa33983]	86	param = GetParam("ApplyLowPUCorr");
	87	fApplyLowPUCorr.clear();
	88	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fApplyLowPUCorr.push_back(param[iMap].GetBool());
	89	// read metric id
	90	param = GetParam("MetricId");
	91	fMetricId.clear();
	92	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fMetricId.push_back(param[iMap].GetInt());
[e921a28]	93	// scheme for combining
	94	param = GetParam("CombId");
	95	fCombId.clear();
	96	for(int iMap = 0; iMap < param.GetSize(); ++iMap) fCombId.push_back(param[iMap].GetInt());
[fa33983]	97	// create output array
	98	fOutputArray = ExportArray(GetString("OutputArray", "puppiParticles"));
	99	fOutputTrackArray = ExportArray(GetString("OutputArrayTracks", "puppiTracks"));
	100	fOutputNeutralArray = ExportArray(GetString("OutputArrayNeutrals", "puppiNeutrals"));
[e921a28]	101	// Create algorithm list for puppi
	102	std::vector<AlgoObj> puppiAlgo;
	103	if(puppiAlgo.empty()){
	104	if(!(fEtaMinBin.size() == fEtaMaxBin.size() and fEtaMinBin.size() == fPtMinBin.size() and fEtaMinBin.size() == fConeSizeBin.size() and fEtaMinBin.size() == fRMSPtMinBin.size()
	105	and fEtaMinBin.size() == fRMSScaleFactorBin.size() and fEtaMinBin.size() == fNeutralMinEBin.size() and fEtaMinBin.size() == fNeutralPtSlope.size()
	106	and fEtaMinBin.size() == fUseCharged.size()
	107	and fEtaMinBin.size() == fApplyLowPUCorr.size() and fEtaMinBin.size() == fMetricId.size())) {
	108	std::cerr<<" Error in PUPPI configuration, algo info should have the same size --> exit from the code"<<std::endl;
	109	std::exit(EXIT_FAILURE);
	110	}
	111	}
	112	for( size_t iAlgo = 0 ; iAlgo < fEtaMinBin.size() ; iAlgo++){
	113	AlgoObj algoTmp ;
	114	algoTmp.etaMin = fEtaMinBin.at(iAlgo);
	115	algoTmp.etaMax = fEtaMaxBin.at(iAlgo);
	116	algoTmp.ptMin = fPtMinBin.at(iAlgo);
	117	algoTmp.minNeutralPt = fNeutralMinEBin.at(iAlgo);
	118	algoTmp.minNeutralPtSlope = fNeutralPtSlope.at(iAlgo);
	119	//Eta Extrapolation stuff is missing
	120	//Loop through file requiring algos for same bins to be adjacent
	121	while(iAlgo < fEtaMinBin.size() and algoTmp.etaMin == fEtaMinBin.at(iAlgo) and algoTmp.etaMax == fEtaMaxBin.at(iAlgo)) {
	122	AlgoSubObj algoSubTmp;
	123	algoSubTmp.metricId = fMetricId.at(iAlgo);
	124	algoSubTmp.useCharged = fUseCharged.at(iAlgo);
	125	algoSubTmp.applyLowPUCorr = fApplyLowPUCorr.at(iAlgo);
	126	algoSubTmp.combId = fCombId.at(iAlgo);
	127	algoSubTmp.coneSize = fConeSizeBin.at(iAlgo);
	128	algoSubTmp.rmsPtMin = fRMSPtMinBin.at(iAlgo);
	129	algoSubTmp.rmsScaleFactor = fRMSScaleFactorBin.at(iAlgo);
	130	algoTmp.subAlgos.push_back(algoSubTmp);
	131	iAlgo++;
	132	}
	133	iAlgo--;
	134	//if(std::find(puppiAlgo.begin(),puppiAlgo.end(),algoTmp) != puppiAlgo.end()) continue;
	135	puppiAlgo.push_back(algoTmp);
	136	}
	137	fPuppi = new PuppiContainer(true,fUseExp,fMinPuppiWeight,puppiAlgo);
[fa33983]	138	}
	139
	140	//------------------------------------------------------------------------------
	141
	142	void RunPUPPI::Finish(){
	143	if(fItTrackInputArray) delete fItTrackInputArray;
	144	if(fItNeutralInputArray) delete fItNeutralInputArray;
	145	}
	146
	147	//------------------------------------------------------------------------------
	148
	149	void RunPUPPI::Process(){
	150
	151	Candidate candidate, particle;
	152	TLorentzVector momentum;
	153
[e921a28]	154	//DelphesFactory *factory = GetFactory();
[fa33983]	155
	156	// loop over input objects
[e921a28]	157	fItTrackInputArray ->Reset();
	158	fItNeutralInputArray ->Reset();
	159	fPVItInputArray ->Reset();
[fa33983]	160
[fb4337d]	161	std::vector<Candidate *> InputParticles;
[fa33983]	162	InputParticles.clear();
	163
	164	// take the leading vertex
	165	float PVZ = 0.;
	166	Candidate pv = static_cast<Candidate>(fPVItInputArray->Next());
	167	if (pv) PVZ = pv->Position.Z();
	168
	169	// Fill input particles for puppi
	170	std::vector<RecoObj> puppiInputVector;
[fb4337d]	171	puppiInputVector.clear();
[e921a28]	172	int lNBad = 0;
[fa33983]	173	// Loop on charge track candidate
	174	while((candidate = static_cast<Candidate*>(fItTrackInputArray->Next()))){
	175	momentum = candidate->Momentum;
	176	RecoObj curRecoObj;
	177	curRecoObj.pt = momentum.Pt();
	178	curRecoObj.eta = momentum.Eta();
	179	curRecoObj.phi = momentum.Phi();
	180	curRecoObj.m = momentum.M();
	181	particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
[e921a28]	182	//if(fApplyNoLep && TMath::Abs(candidate->PID) == 11) continue; //Dumb cut to minimize the nolepton on electron
	183	//if(fApplyNoLep && TMath::Abs(candidate->PID) == 13) continue;
[fa33983]	184	if (candidate->IsRecoPU and candidate->Charge !=0) { // if it comes fromPU vertexes after the resolution smearing and the dZ matching within resolution
[e921a28]	185	lNBad++;
[fa33983]	186	curRecoObj.id = 2;
[e921a28]	187	curRecoObj.vtxId = 0.7*(fPVInputArray->GetEntries()); //Hack apply reco vtx efficiency of 70% for calibration
[fa33983]	188	if(TMath::Abs(candidate->PID) == 11) curRecoObj.pfType = 2;
	189	else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
	190	else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
	191	else curRecoObj.pfType = 1;
	192	curRecoObj.dZ = particle->Position.Z()-PVZ;
	193	}
	194	else if(!candidate->IsRecoPU and candidate->Charge !=0) {
	195	curRecoObj.id = 1; // charge from LV
	196	curRecoObj.vtxId = 1; // from PV
	197	if(TMath::Abs(candidate->PID) == 11) curRecoObj.pfType = 2;
	198	else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
	199	else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
	200	else curRecoObj.pfType = 1;
	201	curRecoObj.dZ = particle->Position.Z()-PVZ;
	202	}
	203	else {
	204	std::cerr<<" RunPUPPI: problem with a charged track --> it has charge 0 "<<std::endl;
	205	continue;
	206	}
	207
	208	puppiInputVector.push_back(curRecoObj);
[fb4337d]	209	InputParticles.push_back(candidate);
[fa33983]	210	}
	211
	212	// Loop on neutral calo cells
	213	while((candidate = static_cast<Candidate*>(fItNeutralInputArray->Next()))){
	214	momentum = candidate->Momentum;
	215	RecoObj curRecoObj;
	216	curRecoObj.pt = momentum.Pt();
	217	curRecoObj.eta = momentum.Eta();
	218	curRecoObj.phi = momentum.Phi();
	219	curRecoObj.m = momentum.M();
[e921a28]	220	curRecoObj.charge = 0;
[fa33983]	221	particle = static_cast<Candidate*>(candidate->GetCandidates()->Last());
	222
	223	if(candidate->Charge == 0){
	224	curRecoObj.id = 0; // neutrals have id==0
	225	curRecoObj.vtxId = 0; // neutrals have vtxId==0
	226	if(TMath::Abs(candidate->PID) == 11) curRecoObj.pfType = 2;
	227	else if(TMath::Abs(candidate->PID) == 13) curRecoObj.pfType = 3;
	228	else if(TMath::Abs(candidate->PID) == 22) curRecoObj.pfType = 4;
	229	else curRecoObj.pfType = 5;
	230	curRecoObj.dZ = particle->Position.Z()-PVZ;
	231	}
	232	else{
	233	std::cerr<<" RunPUPPI: problem with a neutrals cells --> it has charge !=0 "<<std::endl;
	234	continue;
	235	}
	236	puppiInputVector.push_back(curRecoObj);
[fb4337d]	237	InputParticles.push_back(candidate);
[fa33983]	238	}
	239	// Create PUPPI container
[e921a28]	240	fPuppi->initialize(puppiInputVector);
	241	fPuppi->puppiWeights();
[659c7b6]	242	std::vector<PseudoJet> puppiParticles = fPuppi->puppiParticles();
[fa33983]	243
	244	// Loop on final particles
[659c7b6]	245	for (std::vector<PseudoJet>::iterator it = puppiParticles.begin() ; it != puppiParticles.end() ; it++) {
[fa33983]	246	if(it->user_index() <= int(InputParticles.size())){
[fb4337d]	247	candidate = static_cast<Candidate *>(InputParticles.at(it->user_index())->Clone());
	248	candidate->Momentum.SetPxPyPzE(it->px(),it->py(),it->pz(),it->e());
[fa33983]	249	fOutputArray->Add(candidate);
	250	if( puppiInputVector.at(it->user_index()).id == 1 or puppiInputVector.at(it->user_index()).id == 2) fOutputTrackArray->Add(candidate);
	251	else if (puppiInputVector.at(it->user_index()).id == 0) fOutputNeutralArray->Add(candidate);
	252	}
	253	else{
	254	std::cerr<<" particle not found in the input Array --> skip "<<std::endl;
	255	continue;
	256	}
	257	}
	258
	259	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: