source: trunk/modules/MadGraphJetParticleSelector.cc@ 18

Last change on this file since 18 was 17, checked in by Pavel Demin, 16 years ago

fix InsertExcludedParticleID

File size: 7.0 KB
Line 
1
2#include "modules/MadGraphJetParticleSelector.h"
3
4
5#include "ExRootAnalysis/ExRootResult.h"
6#include "ExRootAnalysis/ExRootClasses.h"
7
8#include "ExRootAnalysis/ExRootFilter.h"
9#include "ExRootAnalysis/ExRootClassifier.h"
10
11#include "ExRootAnalysis/ExRootFactory.h"
12#include "ExRootAnalysis/ExRootCandidate.h"
13
14#include "TMath.h"
15#include "TString.h"
16#include "TLorentzVector.h"
17#include "TClonesArray.h"
18
19#include <iostream>
20#include <set>
21
22using namespace std;
23
24
25//------------------------------------------------------------------------------
26
27class MadGraphJetParticleClassifier : public ExRootClassifier
28{
29public:
30
31 MadGraphJetParticleClassifier(TClonesArray *branch);
32
33 Int_t GetCategory(TObject *object);
34
35 void SetEtaMax(Double_t eta);
36 void InsertSpecialParticleID(Int_t pid);
37 void InsertExcludedAncestorID(Int_t pid);
38 void InsertExcludedParticleID(Int_t pid);
39 void SetHadronizationInfo(Bool_t info);
40
41private:
42
43 Bool_t hasBadAncestor(ExRootGenParticle *object);
44
45 Double_t fEtaMax;
46
47 TClonesArray *fBranchParticle;
48
49 set< Int_t > fSpecialParticleIDSet;
50 set< Int_t > fExcludedAncestorIDSet;
51 set< Int_t > fExcludedParticleIDSet;
52};
53
54//------------------------------------------------------------------------------
55
56MadGraphJetParticleClassifier::MadGraphJetParticleClassifier(TClonesArray *branch) :
57 fBranchParticle(branch)
58{
59}
60
61//------------------------------------------------------------------------------
62
63void MadGraphJetParticleClassifier::SetEtaMax(Double_t eta)
64{
65 fEtaMax = eta;
66}
67
68//------------------------------------------------------------------------------
69
70void MadGraphJetParticleClassifier::InsertSpecialParticleID(Int_t pid)
71{
72 fSpecialParticleIDSet.insert(pid);
73}
74
75//------------------------------------------------------------------------------
76
77void MadGraphJetParticleClassifier::InsertExcludedAncestorID(Int_t pid)
78{
79 fExcludedAncestorIDSet.insert(pid);
80}
81
82//------------------------------------------------------------------------------
83
84void MadGraphJetParticleClassifier::InsertExcludedParticleID(Int_t pid)
85{
86 fExcludedParticleIDSet.insert(pid);
87}
88
89//------------------------------------------------------------------------------
90
91Bool_t MadGraphJetParticleClassifier::hasBadAncestor(ExRootGenParticle *object)
92{
93 const int kMaxAncestors = 10;
94 Int_t i, pidAbs;
95 ExRootGenParticle *particle = object;
96 set< Int_t >::const_iterator itAncestorIDSet;
97 set< Int_t >::const_iterator itParticleIDSet;
98
99 for(i = 0; i < kMaxAncestors; ++i)
100 {
101 if(particle->M1 < 0) return kFALSE;
102
103 particle = static_cast<ExRootGenParticle*>(fBranchParticle->At(particle->M1));
104
105 pidAbs = TMath::Abs(particle->PID);
106
107 // skip particles with pid included in list
108 itAncestorIDSet = fExcludedAncestorIDSet.find(pidAbs);
109
110 if(itAncestorIDSet != fExcludedAncestorIDSet.end()) return kTRUE;
111 }
112
113 return kFALSE;
114}
115
116//------------------------------------------------------------------------------
117
118Int_t MadGraphJetParticleClassifier::GetCategory(TObject *object)
119{
120 ExRootGenParticle *particle = static_cast<ExRootGenParticle*>(object);
121
122 set< Int_t >::const_iterator itParticleIDSet;
123
124 Int_t pidAbs = TMath::Abs(particle->PID);
125 Double_t etaAbs = TMath::Abs(particle->Eta);
126
127 // skip beam particles and initial state partons
128 if(particle->M1 < 2) return -1;
129
130 if(etaAbs > fEtaMax) return -1;
131
132 if(particle->Status != 1) return -1;
133
134 // skip particles with pid included in list
135 itParticleIDSet = fExcludedParticleIDSet.find(pidAbs);
136
137 if(itParticleIDSet != fExcludedParticleIDSet.end()) return -1;
138
139 // check ancestors for particles with pid included in list
140 itParticleIDSet = fSpecialParticleIDSet.find(pidAbs);
141
142 if(itParticleIDSet != fSpecialParticleIDSet.end())
143 {
144 if(hasBadAncestor(particle)) return -1;
145 }
146
147 return 0;
148}
149
150//------------------------------------------------------------------------------
151
152MadGraphJetParticleSelector::MadGraphJetParticleSelector() :
153 fFilter(0), fClassifier(0)
154{
155}
156
157//------------------------------------------------------------------------------
158
159MadGraphJetParticleSelector::~MadGraphJetParticleSelector()
160{
161}
162
163//------------------------------------------------------------------------------
164
165void MadGraphJetParticleSelector::Init()
166{
167 ExRootConfParam param;
168
169 Int_t i, pid, sizeParam;
170
171 // import ROOT tree branch
172
173 fBranchParticle = UseBranch("GenParticle");
174
175 // create classifier and filter
176
177 fClassifier = new MadGraphJetParticleClassifier(fBranchParticle);
178 fFilter = new ExRootFilter(fBranchParticle);
179
180 fEtaMax = GetDouble("EtaMax", 5.0);
181 fClassifier->SetEtaMax(fEtaMax);
182
183 // read particle IDs from configuration file and setup classifier
184
185 param = GetParam("SpecialParticleIDs");
186 sizeParam = param.GetSize();
187
188 for(i = 0; i < sizeParam; ++i)
189 {
190 pid = param[i].GetInt();
191 fClassifier->InsertSpecialParticleID(pid);
192 }
193
194 // read ancestor IDs from configuration file and setup classifier
195
196 param = GetParam("ExcludedAncestorIDs");
197 sizeParam = param.GetSize();
198
199 for(i = 0; i < sizeParam; ++i)
200 {
201 pid = param[i].GetInt();
202 fClassifier->InsertExcludedAncestorID(pid);
203 }
204
205 // read particle IDs from configuration file and setup classifier
206
207 param = GetParam("ExcludedParticleIDs");
208 sizeParam = param.GetSize();
209
210 for(i = 0; i < sizeParam; ++i)
211 {
212 pid = param[i].GetInt();
213 fClassifier->InsertExcludedParticleID(pid);
214 }
215
216 // create output arrays
217
218 fOutputArray = ExportArray("candidates");
219
220}
221
222//------------------------------------------------------------------------------
223
224void MadGraphJetParticleSelector::Finish()
225{
226 if(fFilter) delete fFilter;
227 if(fClassifier) delete fClassifier;
228}
229
230//------------------------------------------------------------------------------
231
232void MadGraphJetParticleSelector::Process()
233{
234 TObjArray *array = 0;
235 ExRootGenParticle *particle = 0;
236 ExRootCandidate *candidate = 0;
237 ExRootFactory *factory = GetFactory();
238
239 TLorentzVector momentum;
240
241 fFilter->Reset();
242 array = fFilter->GetSubArray(fClassifier, 0);
243
244 if(array == 0) return;
245
246 TIter itArray(array);
247
248 while(particle = static_cast<ExRootGenParticle*>(itArray.Next()))
249 {
250 momentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
251
252 candidate = factory->NewCandidate();
253
254 candidate->SetP4(momentum);
255 candidate->SetType(particle->PID);
256
257 fOutputArray->Add(candidate);
258 }
259
260/*
261 cout << "==============================" << endl;
262 Int_t indexParticle = -1;
263 itArray.Reset();
264 while(particle = static_cast<ExRootGenParticle*>(itArray.Next()))
265 {
266 ++indexParticle;
267 cout << "--->\t" << particle->Status << "\t" << particle->PID << "\t";
268 cout << particle->M1 << "\t" << particle->M2 << "\t";
269 cout << particle->Px << "\t" << particle->Py << "\t" << particle->Pz << endl;
270 }
271*/
272}
273
274//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.