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 |
|
---|
50 | using namespace std;
|
---|
51 |
|
---|
52 | //------------------------------------------------------------------------------
|
---|
53 |
|
---|
54 | PdgCodeFilter::PdgCodeFilter() :
|
---|
55 | fItInputArray(0)
|
---|
56 | {
|
---|
57 | }
|
---|
58 |
|
---|
59 | //------------------------------------------------------------------------------
|
---|
60 |
|
---|
61 | PdgCodeFilter::~PdgCodeFilter()
|
---|
62 | {
|
---|
63 | }
|
---|
64 |
|
---|
65 | //------------------------------------------------------------------------------
|
---|
66 |
|
---|
67 | void 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 |
|
---|
111 | void PdgCodeFilter::Finish()
|
---|
112 | {
|
---|
113 | if(fItInputArray) delete fItInputArray;
|
---|
114 | }
|
---|
115 |
|
---|
116 | //------------------------------------------------------------------------------
|
---|
117 |
|
---|
118 | void 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 | //std::cout << "--" << pdgCode << std::endl;
|
---|
130 |
|
---|
131 | if (fRequireKeepGhostBHadron) {
|
---|
132 | if (isBHadron(abs(pdgCode)) ){
|
---|
133 | candidate->PT = candidate->PT * 1e-18;
|
---|
134 | if (candidate->PT ==0) candidate->PT = 1e-18;
|
---|
135 | candidate->Momentum.SetPtEtaPhiM(candidate->PT, candidate->Momentum.Eta(), candidate->Momentum.Phi(), candidate->Momentum.M());
|
---|
136 | fOutputArray->Add(candidate);
|
---|
137 | continue;
|
---|
138 | }
|
---|
139 | }
|
---|
140 |
|
---|
141 | const TLorentzVector &candidateMomentum = candidate->Momentum;
|
---|
142 | pt = candidateMomentum.Pt();
|
---|
143 |
|
---|
144 | if(pt < fPTMin) continue;
|
---|
145 | if(fRequireStatus && (candidate->Status != fStatus)) continue;
|
---|
146 | if(fRequireCharge && (candidate->Charge != fCharge)) continue;
|
---|
147 | if(fRequireNotPileup && (candidate->IsPU > 0)) continue;
|
---|
148 |
|
---|
149 | pass = kTRUE;
|
---|
150 | if(find(fPdgCodes.begin(), fPdgCodes.end(), pdgCode) != fPdgCodes.end()) pass = kFALSE;
|
---|
151 |
|
---|
152 | if(fInvert) pass = !pass;
|
---|
153 | if(pass) fOutputArray->Add(candidate);
|
---|
154 | }
|
---|
155 | }
|
---|
156 |
|
---|
157 | Bool_t PdgCodeFilter::isBHadron(const unsigned int absPdgId) {
|
---|
158 | if (absPdgId <= 100)
|
---|
159 | return false; // Fundamental particles and MC internals
|
---|
160 | if (absPdgId >= 1000000000)
|
---|
161 | return false; // Nuclei, +-10LZZZAAAI
|
---|
162 |
|
---|
163 | // General form of PDG ID is 7 digit form
|
---|
164 | // +- n nr nL nq1 nq2 nq3 nJ
|
---|
165 | //const int nJ = absPdgId % 10; // Spin
|
---|
166 | const int nq3 = (absPdgId / 10) % 10;
|
---|
167 | const int nq2 = (absPdgId / 100) % 10;
|
---|
168 | const int nq1 = (absPdgId / 1000) % 10;
|
---|
169 |
|
---|
170 | if (nq3 == 0)
|
---|
171 | return false; // Diquarks
|
---|
172 | if (nq1 == 0 and nq2 == 5)
|
---|
173 | return true; // B mesons
|
---|
174 | if (nq1 == 5)
|
---|
175 | return true; // B baryons
|
---|
176 |
|
---|
177 | return false;
|
---|
178 | }
|
---|
179 |
|
---|