Fork me on GitHub

Context Navigation

source: git/modules/PileUpJetID.cc@ de6d698

Visit:

ImprovedOutputFile Timing dual_readout llp

Last change on this file since de6d698 was 2862770, checked in by Michele Selvaggi <michele.selvaggi@…>, 10 years ago
included last PileUpJetID version from Seth
Property mode set to `100644`
File size: 10.2 KB

Line
1
2	/** \class PileUpJetID
3	*
4	* CMS PileUp Jet ID Variables
5	*
6	* \author S. Zenz
7	*
8	*/
9
10	#include "modules/PileUpJetID.h"
11
12	#include "classes/DelphesClasses.h"
13	#include "classes/DelphesFactory.h"
14	#include "classes/DelphesFormula.h"
15
16	#include "ExRootAnalysis/ExRootResult.h"
17	#include "ExRootAnalysis/ExRootFilter.h"
18	#include "ExRootAnalysis/ExRootClassifier.h"
19
20	#include "TMath.h"
21	#include "TString.h"
22	#include "TFormula.h"
23	#include "TRandom3.h"
24	#include "TObjArray.h"
25	//#include "TDatabasePDG.h"
26	#include "TLorentzVector.h"
27
28	#include <algorithm>
29	#include <stdexcept>
30	#include <iostream>
31	#include <sstream>
32
33	using namespace std;
34
35	//------------------------------------------------------------------------------
36
37	PileUpJetID::PileUpJetID() :
38	fItJetInputArray(0),fTrackInputArray(0),fNeutralInputArray(0)
39	{
40
41	}
42
43	//------------------------------------------------------------------------------
44
45	PileUpJetID::~PileUpJetID()
46	{
47
48	}
49
50	//------------------------------------------------------------------------------
51
52	void PileUpJetID::Init()
53	{
54
55	fJetPTMin = GetDouble("JetPTMin", 20.0);
56	fParameterR = GetDouble("ParameterR", 0.5);
57	fUseConstituents = GetInt("UseConstituents", 0);
58
59
60	/*
61	Double_t fMeanSqDeltaRMaxBarrel; // \|eta\| < 1.5
62	Double_t fBetaMinBarrel; // \|eta\| < 2.5
63	Double_t fMeanSqDeltaRMaxEndcap; // 1.5 < \|eta\| < 4.0
64	Double_t fBetaMinEndcap; // 1.5 < \|eta\| < 4.0
65	Double_t fMeanSqDeltaRMaxForward; // \|eta\| > 4.0
66	*/
67
68	fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
69	fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
70	fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
71	fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
72	fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
73	fJetPTMinForNeutrals = GetDouble("JetPTMinForNeutrals", 20.0);
74	fNeutralPTMin = GetDouble("NeutralPTMin", 2.0);
75
76
77
78	cout << " set MeanSqDeltaRMaxBarrel " << fMeanSqDeltaRMaxBarrel << endl;
79	cout << " set BetaMinBarrel " << fBetaMinBarrel << endl;
80	cout << " set MeanSqDeltaRMaxEndcap " << fMeanSqDeltaRMaxEndcap << endl;
81	cout << " set BetaMinEndcap " << fBetaMinEndcap << endl;
82	cout << " set MeanSqDeltaRMaxForward " << fMeanSqDeltaRMaxForward << endl;
83
84
85
86	fAverageEachTower = false; // for timing
87
88	// import input array(s)
89
90	fJetInputArray = ImportArray(GetString("JetInputArray", "FastJetFinder/jets"));
91	fItJetInputArray = fJetInputArray->MakeIterator();
92
93
94	// cout << "BeforE SCZ additions in init" << endl;
95	// cout << GetString("TrackInputArray", "ParticlePropagator/tracks") << endl;
96	// cout << GetString("EFlowTrackInputArray", "ParticlePropagator/tracks") << endl;
97
98	fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
99	fItTrackInputArray = fTrackInputArray->MakeIterator();
100
101	fNeutralInputArray = ImportArray(GetString("NeutralInputArray", "ParticlePropagator/tracks"));
102	fItNeutralInputArray = fNeutralInputArray->MakeIterator();
103
104
105	// create output array(s)
106
107	fOutputArray = ExportArray(GetString("OutputArray", "jets"));
108
109	fNeutralsInPassingJets = ExportArray(GetString("NeutralsInPassingJets","eflowtowers"));
110
111
112	// cout << " end of INIT " << endl;
113
114	}
115
116	//------------------------------------------------------------------------------
117
118	void PileUpJetID::Finish()
119	{
120	// cout << "In finish" << endl;
121
122	if(fItJetInputArray) delete fItJetInputArray;
123	if(fItTrackInputArray) delete fItTrackInputArray;
124	if(fItNeutralInputArray) delete fItNeutralInputArray;
125
126	}
127
128	//------------------------------------------------------------------------------
129
130	void PileUpJetID::Process()
131	{
132	// cout << "start of process" << endl;
133
134	Candidate candidate, constituent;
135	TLorentzVector momentum, area;
136
137	// cout << "BeforE SCZ additions in process" << endl;
138
139	// SCZ
140	Candidate *trk;
141
142	// loop over all input candidates
143	fItJetInputArray->Reset();
144	while((candidate = static_cast<Candidate*>(fItJetInputArray->Next())))
145	{
146	momentum = candidate->Momentum;
147	area = candidate->Area;
148
149	float sumT0 = 0.;
150	float sumT1 = 0.;
151	float sumT10 = 0.;
152	float sumT20 = 0.;
153	float sumT30 = 0.;
154	float sumT40 = 0.;
155	float sumWeightsForT = 0.;
156	candidate->Ntimes = 0;
157
158	float sumpt = 0.;
159	float sumptch = 0.;
160	float sumptchpv = 0.;
161	float sumptchpu = 0.;
162	float sumdrsqptsq = 0.;
163	float sumptsq = 0.;
164	int nc = 0;
165	int nn = 0;
166	float pt_ann[5];
167
168	for (int i = 0 ; i < 5 ; i++) {
169	pt_ann[i] = 0.;
170	}
171
172	if (fUseConstituents) {
173	TIter itConstituents(candidate->GetCandidates());
174	while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
175	float pt = constituent->Momentum.Pt();
176	float dr = candidate->Momentum.DeltaR(constituent->Momentum);
177	// cout << " There exists a constituent with dr=" << dr << endl;
178	sumpt += pt;
179	sumdrsqptsq += drdrpt*pt;
180	sumptsq += pt*pt;
181	if (constituent->Charge == 0) {
182	nn++;
183	} else {
184	if (constituent->IsRecoPU) {
185	sumptchpu += pt;
186	} else {
187	sumptchpv += pt;
188	}
189	sumptch += pt;
190	nc++;
191	}
192	for (int i = 0 ; i < 5 ; i++) {
193	if (dr > 0.1i && dr < 0.1(i+1)) {
194	pt_ann[i] += pt;
195	}
196	}
197	float tow_sumT = 0;
198	float tow_sumW = 0;
199	for (int i = 0 ; i < constituent->Ecal_E_t.size() ; i++) {
200	float w = TMath::Sqrt(constituent->Ecal_E_t[i].first);
201	if (fAverageEachTower) {
202	tow_sumT += w*constituent->Ecal_E_t[i].second;
203	tow_sumW += w;
204	} else {
205	sumT0 += w*constituent->Ecal_E_t[i].second;
206	sumT1 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.001);
207	sumT10 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.010);
208	sumT20 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.020);
209	sumT30 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.030);
210	sumT40 += w*gRandom->Gaus(constituent->Ecal_E_t[i].second,0.040);
211	sumWeightsForT += w;
212	candidate->Ntimes++;
213	}
214	}
215	if (fAverageEachTower && tow_sumW > 0.) {
216	sumT0 += tow_sumT;
217	sumT1 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.001);
218	sumT10 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0010);
219	sumT20 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0020);
220	sumT30 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0030);
221	sumT40 += tow_sumW*gRandom->Gaus(tow_sumT/tow_sumW,0.0040);
222	sumWeightsForT += tow_sumW;
223	candidate->Ntimes++;
224	}
225	}
226	} else {
227	// Not using constituents, using dr
228	fItTrackInputArray->Reset();
229	while ((trk = static_cast<Candidate*>(fItTrackInputArray->Next()))) {
230	if (trk->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
231	float pt = trk->Momentum.Pt();
232	sumpt += pt;
233	sumptch += pt;
234	if (trk->IsRecoPU) {
235	sumptchpu += pt;
236	} else {
237	sumptchpv += pt;
238	}
239	float dr = candidate->Momentum.DeltaR(trk->Momentum);
240	sumdrsqptsq += drdrpt*pt;
241	sumptsq += pt*pt;
242	nc++;
243	for (int i = 0 ; i < 5 ; i++) {
244	if (dr > 0.1i && dr < 0.1(i+1)) {
245	pt_ann[i] += pt;
246	}
247	}
248	}
249	}
250	fItNeutralInputArray->Reset();
251	while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
252	if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR) {
253	float pt = constituent->Momentum.Pt();
254	sumpt += pt;
255	float dr = candidate->Momentum.DeltaR(constituent->Momentum);
256	sumdrsqptsq += drdrpt*pt;
257	sumptsq += pt*pt;
258	nn++;
259	for (int i = 0 ; i < 5 ; i++) {
260	if (dr > 0.1i && dr < 0.1(i+1)) {
261	pt_ann[i] += pt;
262	}
263	}
264	}
265	}
266	}
267
268	if (sumptch > 0.) {
269	candidate->Beta = sumptchpv/sumptch;
270	candidate->BetaStar = sumptchpu/sumptch;
271	} else {
272	candidate->Beta = -999.;
273	candidate->BetaStar = -999.;
274	}
275	if (sumptsq > 0.) {
276	candidate->MeanSqDeltaR = sumdrsqptsq/sumptsq;
277	} else {
278	candidate->MeanSqDeltaR = -999.;
279	}
280	candidate->NCharged = nc;
281	candidate->NNeutrals = nn;
282	if (sumpt > 0.) {
283	candidate->PTD = TMath::Sqrt(sumptsq) / sumpt;
284	for (int i = 0 ; i < 5 ; i++) {
285	candidate->FracPt[i] = pt_ann[i]/sumpt;
286	}
287	} else {
288	candidate->PTD = -999.;
289	for (int i = 0 ; i < 5 ; i++) {
290	candidate->FracPt[i] = -999.;
291	}
292	}
293
294	fOutputArray->Add(candidate);
295
296	// New stuff
297	/*
298	fMeanSqDeltaRMaxBarrel = GetDouble("MeanSqDeltaRMaxBarrel",0.1);
299	fBetaMinBarrel = GetDouble("BetaMinBarrel",0.1);
300	fMeanSqDeltaRMaxEndcap = GetDouble("MeanSqDeltaRMaxEndcap",0.1);
301	fBetaMinEndcap = GetDouble("BetaMinEndcap",0.1);
302	fMeanSqDeltaRMaxForward = GetDouble("MeanSqDeltaRMaxForward",0.1);
303	*/
304
305	bool passId = false;
306	if (candidate->Momentum.Pt() > fJetPTMinForNeutrals && candidate->MeanSqDeltaR > -0.1) {
307	if (fabs(candidate->Momentum.Eta())<1.5) {
308	passId = ((candidate->Beta > fBetaMinBarrel) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxBarrel));
309	} else if (fabs(candidate->Momentum.Eta())<4.0) {
310	passId = ((candidate->Beta > fBetaMinEndcap) && (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxEndcap));
311	} else {
312	passId = (candidate->MeanSqDeltaR < fMeanSqDeltaRMaxForward);
313	}
314	}
315
316	// cout << " Pt Eta MeanSqDeltaR Beta PassId " << candidate->Momentum.Pt()
317	// << " " << candidate->Momentum.Eta() << " " << candidate->MeanSqDeltaR << " " << candidate->Beta << " " << passId << endl;
318
319	if (passId) {
320	if (fUseConstituents) {
321	TIter itConstituents(candidate->GetCandidates());
322	while((constituent = static_cast<Candidate*>(itConstituents.Next()))) {
323	if (constituent->Charge == 0 && constituent->Momentum.Pt() > fNeutralPTMin) {
324	fNeutralsInPassingJets->Add(constituent);
325	// cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
326	}
327	}
328	} else { // use DeltaR
329	fItNeutralInputArray->Reset();
330	while ((constituent = static_cast<Candidate*>(fItNeutralInputArray->Next()))) {
331	if (constituent->Momentum.DeltaR(candidate->Momentum) < fParameterR && constituent->Momentum.Pt() > fNeutralPTMin) {
332	fNeutralsInPassingJets->Add(constituent);
333	// cout << " Constitutent added Pt Eta Charge " << constituent->Momentum.Pt() << " " << constituent->Momentum.Eta() << " " << constituent->Charge << endl;
334	}
335	}
336	}
337	}
338
339
340	}
341	}
342
343	//------------------------------------------------------------------------------

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

Download in other formats: