Fork me on GitHub

source: git/modules/StatusPidFilter.cc@ 6ba7847

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 6ba7847 was f3598d0, checked in by Fabio Monti <fabio.monti@…>, 7 years ago

add stable gen photon info

  • Property mode set to 100644
File size: 5.8 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/ExRootResult.h"
35#include "ExRootAnalysis/ExRootFilter.h"
36#include "ExRootAnalysis/ExRootClassifier.h"
37
38#include "TMath.h"
39#include "TString.h"
40#include "TFormula.h"
41#include "TRandom3.h"
42#include "TObjArray.h"
43#include "TDatabasePDG.h"
44#include "TLorentzVector.h"
45
46#include <algorithm>
47#include <stdexcept>
48#include <iostream>
49#include <sstream>
50
51using namespace std;
52
53namespace {
54 // integer power (faster than TMath::Pow() + cast to integer)
55 int ipow(int base, int exp)
56 {
57 int result = 1;
58 while (exp)
59 {
60 if (exp & 1)
61 result *= base;
62 exp >>= 1;
63 base *= base;
64 }
65
66 return result;
67 }
68
69 // standalone function to extract the i-th digit from a number (counting from 0 = rightmost, etc..)
70 int digit(int val, int i)
71 {
72 int y = ipow(10,i);
73 int z = val/y;
74 int val2 = val / (y * 10);
75 return (z - val2*10);
76 }
77
78 // return the first two digits if this is a "fundamental" particle
79 // ID = 100 is a special case (internal generator ID's are 81-100)
80 // also, 101 and 102 are now used (by HepPID) for geantinos
81 int fundamentalID(int pdgCode)
82 {
83 pdgCode = abs(pdgCode);
84 if( ( digit(pdgCode, 9) == 1 ) && ( digit(pdgCode, 8) == 0 ) ) { return 0; }
85 if( digit(pdgCode, 2) == 0 && digit(pdgCode, 3) == 0) {
86 return pdgCode%10000;
87 } else if( pdgCode <= 102 ) {
88 return pdgCode;
89 } else {
90 return 0;
91 }
92 }
93
94 bool hasBottom(int pdgCode)
95 {
96 if ( (pdgCode/10000000) > 0)
97 return false;
98 if (pdgCode <= 100)
99 return false;
100 if (fundamentalID(pdgCode) <= 100 && fundamentalID(pdgCode) > 0)
101 return false;
102 if( digit(pdgCode, 3) == 5 || digit(pdgCode, 2) == 5 || digit(pdgCode, 1) == 5 )
103 return true;
104 return false;
105 }
106
107 bool isTauDaughter(int pdgCode, int M1, const TObjArray *fInputArray){
108 //not needed, just to speed up the code - can be further refined but gives only negligible improvement:
109 if ( pdgCode==15 || pdgCode<11 || (pdgCode > 22 && pdgCode < 100) || pdgCode>1000 )
110 return false;
111
112 if ( M1 < 0 )
113 return false;
114
115 Candidate *mother;
116 mother = static_cast<Candidate*>(fInputArray->At( M1 ));
117 if ( TMath::Abs(mother->PID) == 15 )
118 return true;
119
120 return false;
121 }
122
123}
124
125
126//------------------------------------------------------------------------------
127
128StatusPidFilter::StatusPidFilter() :
129 fItInputArray(0)
130{
131}
132
133//------------------------------------------------------------------------------
134
135StatusPidFilter::~StatusPidFilter()
136{
137}
138
139//------------------------------------------------------------------------------
140
141void StatusPidFilter::Init()
142{
143 // PT threshold
144
145 fPTMin = GetDouble("PTMin", 0.5);
146
147 // import input array
148 fInputArray = ImportArray(GetString("InputArray", "Delphes/allParticles"));
149 fItInputArray = fInputArray->MakeIterator();
150
151 // create output array
152
153 fOutputArray = ExportArray(GetString("OutputArray", "filteredParticles"));
154}
155
156//------------------------------------------------------------------------------
157
158void StatusPidFilter::Finish()
159{
160 if(fItInputArray) delete fItInputArray;
161}
162
163//------------------------------------------------------------------------------
164
165void StatusPidFilter::Process()
166{
167 Candidate *candidate;
168 Int_t status, pdgCode;
169 Bool_t pass;
170
171 fItInputArray->Reset();
172 while((candidate = static_cast<Candidate*>(fItInputArray->Next())))
173 {
174 status = candidate->Status;
175 pdgCode = TMath::Abs(candidate->PID);
176
177 pass = kFALSE;
178
179 // hard scattering particles (first condition for Py6, second for Py8)
180 if(status == 3) pass = kTRUE;
181 if(status > 20 && status < 30 ) pass = kTRUE;
182
183 // electrons, muons, taus and neutrinos
184 if(pdgCode > 10 && pdgCode < 17) pass = kTRUE;
185
186 // heavy quarks
187 if(pdgCode == 4 ||pdgCode == 5 || pdgCode == 6) pass = kTRUE;
188
189 // Gauge bosons and other fundamental bosons
190 if(pdgCode > 22 && pdgCode < 43) pass = kTRUE;
191
192 //Stable photons
193 if(pdgCode == 22 && status==1) pass = kTRUE;
194
195 // logic ported from HepPDF: http://lcgapp.cern.ch/project/simu/HepPDT/HepPDT.2.05.02/html/ParticleID_8cc-source.html#l00081
196 bool is_b_hadron = hasBottom(pdgCode);
197 bool is_b_quark = (pdgCode == 5);
198
199 bool is_tau_daughter = isTauDaughter(pdgCode, candidate->M1, fInputArray);
200
201 if (is_b_hadron)
202 pass = kTRUE;
203
204 if (is_tau_daughter)
205 pass = kTRUE;
206
207 // fPTMin not applied to b_hadrons / b_quarks to allow for b-enriched sample stitching
208 // fPTMin not applied to tau decay products to allow visible-tau four momentum determination
209 if(!pass || (candidate->Momentum.Pt() < fPTMin && !(is_b_hadron || is_b_quark || is_tau_daughter)) ) continue;
210
211 fOutputArray->Add(candidate);
212 }
213}
214
Note: See TracBrowser for help on using the repository browser.