Fork me on GitHub

source: git/modules/PdgCodeFilter.cc

Last change on this file was 54ec182, checked in by youngkwon jo <cccpy98@…>, 3 years ago

remove a comment

  • Property mode set to 100644
File size: 4.7 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 PdgCodeFilter
20 *
21 * Removes particles with specific PDG codes
22 *
23 * \author M. Selvaggi
24 *
25 */
26
27#include "modules/PdgCodeFilter.h"
28
29#include "classes/DelphesClasses.h"
30#include "classes/DelphesFactory.h"
31#include "classes/DelphesFormula.h"
32
33#include "ExRootAnalysis/ExRootClassifier.h"
34#include "ExRootAnalysis/ExRootFilter.h"
35#include "ExRootAnalysis/ExRootResult.h"
36
37#include "TDatabasePDG.h"
38#include "TFormula.h"
39#include "TLorentzVector.h"
40#include "TMath.h"
41#include "TObjArray.h"
42#include "TRandom3.h"
43#include "TString.h"
44
45#include <algorithm>
46#include <iostream>
47#include <sstream>
48#include <stdexcept>
49
50using namespace std;
51
52//------------------------------------------------------------------------------
53
54PdgCodeFilter::PdgCodeFilter() :
55 fItInputArray(0)
56{
57}
58
59//------------------------------------------------------------------------------
60
61PdgCodeFilter::~PdgCodeFilter()
62{
63}
64
65//------------------------------------------------------------------------------
66
67void PdgCodeFilter::Init()
68{
69
70 ExRootConfParam param;
71 Size_t i, size;
72
73 // PT threshold
74 fPTMin = GetDouble("PTMin", 0.0);
75
76 fInvert = GetBool("Invert", false);
77
78 // no pileup
79 fRequireNotPileup = GetBool("RequireNotPileup", false);
80
81 fRequireStatus = GetBool("RequireStatus", false);
82 fStatus = GetInt("Status", 1);
83
84 fRequireCharge = GetBool("RequireCharge", false);
85 fCharge = GetInt("Charge", 1);
86
87 // keep bhadron
88 fRequireKeepGhostBHadron = GetBool("RequireKeepGhostBHadron", false);
89
90 // import input array
91 fInputArray = ImportArray(GetString("InputArray", "Delphes/allParticles"));
92 fItInputArray = fInputArray->MakeIterator();
93
94 param = GetParam("PdgCode");
95 size = param.GetSize();
96
97 // read PdgCodes to be filtered out from the data card
98
99 fPdgCodes.clear();
100 for(i = 0; i < size; ++i)
101 {
102 fPdgCodes.push_back(param[i].GetInt());
103 }
104
105 // create output array
106 fOutputArray = ExportArray(GetString("OutputArray", "filteredParticles"));
107}
108
109//------------------------------------------------------------------------------
110
111void PdgCodeFilter::Finish()
112{
113 if(fItInputArray) delete fItInputArray;
114}
115
116//------------------------------------------------------------------------------
117
118void PdgCodeFilter::Process()
119{
120 Candidate *candidate;
121 Int_t pdgCode;
122 Bool_t pass;
123 Double_t pt;
124
125 fItInputArray->Reset();
126 while((candidate = static_cast<Candidate *>(fItInputArray->Next())))
127 {
128 pdgCode = candidate->PID;
129
130 if (fRequireKeepGhostBHadron) {
131 if (isBHadron(abs(pdgCode)) ){
132 candidate->PT = candidate->PT * 1e-18;
133 if (candidate->PT ==0) candidate->PT = 1e-18;
134 candidate->Momentum.SetPtEtaPhiM(candidate->PT, candidate->Momentum.Eta(), candidate->Momentum.Phi(), candidate->Momentum.M());
135 fOutputArray->Add(candidate);
136 continue;
137 }
138 }
139
140 const TLorentzVector &candidateMomentum = candidate->Momentum;
141 pt = candidateMomentum.Pt();
142
143 if(pt < fPTMin) continue;
144 if(fRequireStatus && (candidate->Status != fStatus)) continue;
145 if(fRequireCharge && (candidate->Charge != fCharge)) continue;
146 if(fRequireNotPileup && (candidate->IsPU > 0)) continue;
147
148 pass = kTRUE;
149 if(find(fPdgCodes.begin(), fPdgCodes.end(), pdgCode) != fPdgCodes.end()) pass = kFALSE;
150
151 if(fInvert) pass = !pass;
152 if(pass) fOutputArray->Add(candidate);
153 }
154}
155
156Bool_t PdgCodeFilter::isBHadron(const unsigned int absPdgId) {
157 if (absPdgId <= 100)
158 return false; // Fundamental particles and MC internals
159 if (absPdgId >= 1000000000)
160 return false; // Nuclei, +-10LZZZAAAI
161
162 // General form of PDG ID is 7 digit form
163 // +- n nr nL nq1 nq2 nq3 nJ
164 //const int nJ = absPdgId % 10; // Spin
165 const int nq3 = (absPdgId / 10) % 10;
166 const int nq2 = (absPdgId / 100) % 10;
167 const int nq1 = (absPdgId / 1000) % 10;
168
169 if (nq3 == 0)
170 return false; // Diquarks
171 if (nq1 == 0 and nq2 == 5)
172 return true; // B mesons
173 if (nq1 == 5)
174 return true; // B baryons
175
176 return false;
177}
178
Note: See TracBrowser for help on using the repository browser.