Fork me on GitHub

Context Navigation

source: git/modules/RunPUPPI.cc@ af2c385

Visit:

ImprovedOutputFile Timing dual_readout llp

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

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

Download in other formats: