Fork me on GitHub

source: git/modules/StatusPidFilter.cc

Last change on this file was 341014c, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

apply .clang-format to all .h, .cc and .cpp files

  • Property mode set to 100644
File size: 6.3 KB
Line 
1/*
2 * Delphes: a framework for fast simulation of a generic collider experiment
3 * Copyright (C) 2012-2014 Universite catholique de Louvain (UCL), Belgium
4 *
5 * This program is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19/** \class StatusPidFilter
20 *
21 * Removes all generated particles except electrons, muons, taus,
22 * and particles with status == 3.
23 *
24 * \author J. Hirschauer - FNAL
25 *
26 */
27
28#include "modules/StatusPidFilter.h"
29
30#include "classes/DelphesClasses.h"
31#include "classes/DelphesFactory.h"
32#include "classes/DelphesFormula.h"
33
34#include "ExRootAnalysis/ExRootClassifier.h"
35#include "ExRootAnalysis/ExRootFilter.h"
36#include "ExRootAnalysis/ExRootResult.h"
37
38#include "TDatabasePDG.h"
39#include "TFormula.h"
40#include "TLorentzVector.h"
41#include "TMath.h"
42#include "TObjArray.h"
43#include "TRandom3.h"
44#include "TString.h"
45
46#include <algorithm>
47#include <iostream>
48#include <sstream>
49#include <stdexcept>
50
51using namespace std;
52
53namespace
54{
55// integer power (faster than TMath::Pow() + cast to integer)
56int ipow(int base, int exp)
57{
58 int result = 1;
59 while(exp)
60 {
61 if(exp & 1)
62 result *= base;
63 exp >>= 1;
64 base *= base;
65 }
66
67 return result;
68}
69
70// standalone function to extract the i-th digit from a number (counting from 0 = rightmost, etc..)
71int digit(int val, int i)
72{
73 int y = ipow(10, i);
74 int z = val / y;
75 int val2 = val / (y * 10);
76 return (z - val2 * 10);
77}
78
79// return the first two digits if this is a "fundamental" particle
80// ID = 100 is a special case (internal generator ID's are 81-100)
81// also, 101 and 102 are now used (by HepPID) for geantinos
82int fundamentalID(int pdgCode)
83{
84 pdgCode = abs(pdgCode);
85 if((digit(pdgCode, 9) == 1) && (digit(pdgCode, 8) == 0))
86 {
87 return 0;
88 }
89 if(digit(pdgCode, 2) == 0 && digit(pdgCode, 3) == 0)
90 {
91 return pdgCode % 10000;
92 }
93 else if(pdgCode <= 102)
94 {
95 return pdgCode;
96 }
97 else
98 {
99 return 0;
100 }
101}
102
103bool hasBottom(int pdgCode)
104{
105 if((pdgCode / 10000000) > 0)
106 return false;
107 if(pdgCode <= 100)
108 return false;
109 if(fundamentalID(pdgCode) <= 100 && fundamentalID(pdgCode) > 0)
110 return false;
111 if(digit(pdgCode, 3) == 5 || digit(pdgCode, 2) == 5 || digit(pdgCode, 1) == 5)
112 return true;
113 return false;
114}
115
116bool isTauDaughter(int pdgCode, int M1, const TObjArray *fInputArray)
117{
118 //not needed, just to speed up the code - can be further refined but gives only negligible improvement:
119 if(pdgCode == 15 || pdgCode < 11 || (pdgCode > 22 && pdgCode < 100) || pdgCode > 1000)
120 return false;
121
122 if(M1 < 0)
123 return false;
124
125 Candidate *mother;
126 mother = static_cast<Candidate *>(fInputArray->At(M1));
127 if(TMath::Abs(mother->PID) == 15)
128 return true;
129
130 return false;
131}
132
133bool isWDaughter(int M1, const TObjArray *fInputArray)
134{
135 if(M1 < 0) return false;
136
137 Candidate *mother;
138 mother = static_cast<Candidate *>(fInputArray->At(M1));
139 if(TMath::Abs(mother->PID) == 24) return true;
140
141 return false;
142}
143
144} // namespace
145
146//------------------------------------------------------------------------------
147
148StatusPidFilter::StatusPidFilter() :
149 fItInputArray(0)
150{
151}
152
153//------------------------------------------------------------------------------
154
155StatusPidFilter::~StatusPidFilter()
156{
157}
158
159//------------------------------------------------------------------------------
160
161void StatusPidFilter::Init()
162{
163 // PT threshold
164 fPTMin = GetDouble("PTMin", 0.5);
165
166 // keep or remove pileup particles
167 fRequireNotPileup = GetBool("RequireNotPileup", false);
168
169 // import input array
170 fInputArray = ImportArray(GetString("InputArray", "Delphes/allParticles"));
171 fItInputArray = fInputArray->MakeIterator();
172
173 // create output array
174
175 fOutputArray = ExportArray(GetString("OutputArray", "filteredParticles"));
176}
177
178//------------------------------------------------------------------------------
179
180void StatusPidFilter::Finish()
181{
182 if(fItInputArray) delete fItInputArray;
183}
184
185//------------------------------------------------------------------------------
186
187void StatusPidFilter::Process()
188{
189 Candidate *candidate;
190 Int_t status, pdgCode;
191 Bool_t pass;
192
193 fItInputArray->Reset();
194 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
195 {
196 status = candidate->Status;
197 pdgCode = TMath::Abs(candidate->PID);
198
199 pass = kFALSE;
200
201 // Store all SUSY particles
202 if(pdgCode >= 1000001 && pdgCode <= 1000039) pass = kTRUE;
203
204 // hard scattering particles (first condition for Py6, second for Py8)
205 if(status == 3) pass = kTRUE;
206 if(status > 20 && status < 30) pass = kTRUE;
207
208 // electrons, muons, taus and neutrinos
209 if(pdgCode > 10 && pdgCode < 17) pass = kTRUE;
210
211 // heavy quarks
212 if(pdgCode == 4 || pdgCode == 5 || pdgCode == 6) pass = kTRUE;
213
214 // Gauge bosons and other fundamental bosons
215 if(pdgCode > 22 && pdgCode < 43) pass = kTRUE;
216
217 //Stable photons
218 if(pdgCode == 22 && status == 1) pass = kTRUE;
219
220 // logic ported from HepPDF: http://lcgapp.cern.ch/project/simu/HepPDT/HepPDT.2.05.02/html/ParticleID_8cc-source.html#l00081
221 bool is_b_hadron = hasBottom(pdgCode);
222 bool is_b_quark = (pdgCode == 5);
223
224 bool is_tau_daughter = isTauDaughter(pdgCode, candidate->M1, fInputArray);
225
226 if(is_b_hadron)
227 pass = kTRUE;
228
229 if(is_tau_daughter)
230 pass = kTRUE;
231
232 bool is_W_daughter = isWDaughter(candidate->M1, fInputArray);
233 if(is_W_daughter)
234 pass = kTRUE;
235
236 // fPTMin not applied to b_hadrons / b_quarks to allow for b-enriched sample stitching
237 // fPTMin not applied to tau decay products to allow visible-tau four momentum determination
238 if(!pass || (candidate->Momentum.Pt() < fPTMin && !(is_b_hadron || is_b_quark || is_tau_daughter || is_W_daughter))) continue;
239
240 // not pileup particles
241 if(fRequireNotPileup && (candidate->IsPU > 0)) continue;
242
243 fOutputArray->Add(candidate);
244 }
245}
Note: See TracBrowser for help on using the repository browser.