source: trunk/modules/MadGraphShowerPartonSelector.cc@ 12

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

add SISCone jet algorithm and update names for ConeJetFinder modules

File size: 9.2 KB
Line 
1
2#include "modules/MadGraphShowerPartonSelector.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 MadGraphShowerPartonClassifier : public ExRootClassifier
28{
29public:
30
31 MadGraphShowerPartonClassifier(TClonesArray *branch);
32
33 Int_t GetCategory(TObject *object);
34
35 void SetEtaMax(Double_t eta);
36 void InsertParticleID(Int_t pid);
37 void InsertExclAncestorID(Int_t pid);
38 void InsertInclAncestorID(Int_t pid);
39 void SetHadronizationInfo(Bool_t info);
40
41private:
42
43 Bool_t hasBadAncestor(ExRootGenParticle *object);
44 Bool_t hasGoodAncestor(ExRootGenParticle *object);
45
46 Double_t fEtaMax;
47
48 TClonesArray *fBranchParticle;
49
50 set< Int_t > fParticleIDSet;
51 set< Int_t > fExclAncestorIDSet;
52 set< Int_t > fInclAncestorIDSet;
53};
54
55//------------------------------------------------------------------------------
56
57MadGraphShowerPartonClassifier::MadGraphShowerPartonClassifier(TClonesArray *branch) :
58 fBranchParticle(branch)
59{
60}
61
62//------------------------------------------------------------------------------
63
64void MadGraphShowerPartonClassifier::SetEtaMax(Double_t eta)
65{
66 fEtaMax = eta;
67}
68
69//------------------------------------------------------------------------------
70
71void MadGraphShowerPartonClassifier::InsertParticleID(Int_t pid)
72{
73 fParticleIDSet.insert(pid);
74}
75
76//------------------------------------------------------------------------------
77
78void MadGraphShowerPartonClassifier::InsertExclAncestorID(Int_t pid)
79{
80 fExclAncestorIDSet.insert(pid);
81}
82
83//------------------------------------------------------------------------------
84
85void MadGraphShowerPartonClassifier::InsertInclAncestorID(Int_t pid)
86{
87 fInclAncestorIDSet.insert(pid);
88}
89
90//------------------------------------------------------------------------------
91
92Bool_t MadGraphShowerPartonClassifier::hasBadAncestor(ExRootGenParticle *object)
93{
94 const int kMaxAncestors = 10;
95 Int_t i, pidAbs;
96 ExRootGenParticle *particle = object;
97 set< Int_t >::const_iterator itAncestorIDSet;
98
99 for(i = 0; i < kMaxAncestors && particle->Status != 3; ++i)
100 {
101 if(particle->M1 < 0) return kFALSE;
102
103 particle = static_cast<ExRootGenParticle*>(fBranchParticle->At(particle->M1));
104 }
105
106 if(particle->PID == 21) return kFALSE;
107
108 // keep all particles if there is no pid in the list
109 if(fExclAncestorIDSet.empty())
110 {
111 return kFALSE;
112 }
113
114 pidAbs = TMath::Abs(particle->PID);
115
116 // skip particles with pid included in list
117 itAncestorIDSet = fExclAncestorIDSet.find(pidAbs);
118
119 if(itAncestorIDSet != fExclAncestorIDSet.end()) return kTRUE;
120 if(particle->M2 > -1) return kFALSE;
121
122 for(i = 0; i < kMaxAncestors; ++i)
123 {
124 if(particle->M1 < 0) return kFALSE;
125
126 if(particle->PID == 21) return kFALSE;
127
128 particle = static_cast<ExRootGenParticle*>(fBranchParticle->At(particle->M1));
129
130 pidAbs = TMath::Abs(particle->PID);
131
132 // skip particles with pid included in list
133 itAncestorIDSet = fExclAncestorIDSet.find(pidAbs);
134
135 if(itAncestorIDSet != fExclAncestorIDSet.end()) return kTRUE;
136 if(particle->M2 > -1) return kFALSE;
137 }
138
139 return kFALSE;
140}
141
142//------------------------------------------------------------------------------
143
144Bool_t MadGraphShowerPartonClassifier::hasGoodAncestor(ExRootGenParticle *object)
145{
146 const int kMaxAncestors = 10;
147 Int_t i, pidAbs;
148 ExRootGenParticle *particle = object;
149 set< Int_t >::const_iterator itAncestorIDSet;
150
151 for(i = 0; i < kMaxAncestors && particle->Status != 3; ++i)
152 {
153 if(particle->M1 < 0) return kFALSE;
154
155 particle = static_cast<ExRootGenParticle*>(fBranchParticle->At(particle->M1));
156 }
157
158 if(particle->PID == 21) return kFALSE;
159
160 // keep all particles if there is no pid in the list
161 if(fInclAncestorIDSet.empty())
162 {
163 return kTRUE;
164 }
165
166 pidAbs = TMath::Abs(particle->PID);
167
168 // keep particles with pid included in list
169 itAncestorIDSet = fInclAncestorIDSet.find(pidAbs);
170
171 if(itAncestorIDSet != fInclAncestorIDSet.end()) return kTRUE;
172 if(particle->M2 > -1) return kFALSE;
173
174 for(i = 0; i < kMaxAncestors; ++i)
175 {
176 if(particle->M1 < 0) return kFALSE;
177
178 if(particle->PID == 21) return kFALSE;
179
180 particle = static_cast<ExRootGenParticle*>(fBranchParticle->At(particle->M1));
181
182 pidAbs = TMath::Abs(particle->PID);
183
184 // keep particles with pid included in list
185 itAncestorIDSet = fInclAncestorIDSet.find(pidAbs);
186
187 if(itAncestorIDSet != fInclAncestorIDSet.end()) return kTRUE;
188 if(particle->M2 > -1) return kFALSE;
189 }
190
191 return kFALSE;
192}
193
194//------------------------------------------------------------------------------
195
196Int_t MadGraphShowerPartonClassifier::GetCategory(TObject *object)
197{
198 ExRootGenParticle *particle = static_cast<ExRootGenParticle*>(object);
199 ExRootGenParticle *daughter;
200
201 set< Int_t >::const_iterator itParticleIDSet;
202
203 Int_t pidAbs = TMath::Abs(particle->PID);
204 Double_t etaAbs = TMath::Abs(particle->Eta);
205
206 // skip beam particles and initial state partons
207 if(particle->M1 < 2) return -1;
208
209 // skip particles with pid not included in list
210 itParticleIDSet = fParticleIDSet.find(pidAbs);
211
212 if(itParticleIDSet == fParticleIDSet.end() || etaAbs > fEtaMax) return -1;
213
214 // with hadronization
215 if(particle->Status == 2)
216 {
217 // skip particles if they do not form a string
218 if(particle->D1 > -1)
219 {
220 daughter = static_cast<ExRootGenParticle*>(fBranchParticle->At(particle->D1));
221 if(daughter->PID != 92) return -1;
222 }
223 }
224 // without hadronization
225 else if(particle->Status != 1) return -1;
226
227 if(hasBadAncestor(particle)) return -1;
228
229 if(!hasGoodAncestor(particle)) return -1;
230
231 return 0;
232}
233
234//------------------------------------------------------------------------------
235
236MadGraphShowerPartonSelector::MadGraphShowerPartonSelector() :
237 fFilter(0), fClassifier(0)
238{
239}
240
241//------------------------------------------------------------------------------
242
243MadGraphShowerPartonSelector::~MadGraphShowerPartonSelector()
244{
245}
246
247//------------------------------------------------------------------------------
248
249void MadGraphShowerPartonSelector::Init()
250{
251 ExRootConfParam param;
252
253 Int_t i, pid, sizeParam;
254
255 // import ROOT tree branch
256
257 fBranchParticle = UseBranch("GenParticle");
258
259 // create classifier and filter
260
261 fClassifier = new MadGraphShowerPartonClassifier(fBranchParticle);
262 fFilter = new ExRootFilter(fBranchParticle);
263
264 fEtaMax = GetDouble("EtaMax", 5.0);
265 fClassifier->SetEtaMax(fEtaMax);
266
267 // read particle IDs from configuration file and setup classifier
268
269 param = GetParam("PartonIDs");
270 sizeParam = param.GetSize();
271
272 for(i = 0; i < sizeParam; ++i)
273 {
274 pid = param[i].GetInt();
275 fClassifier->InsertParticleID(pid);
276 }
277
278 // read ancestor IDs from configuration file and setup classifier
279
280 param = GetParam("ExcludedAncestorIDs");
281 sizeParam = param.GetSize();
282
283 for(i = 0; i < sizeParam; ++i)
284 {
285 pid = param[i].GetInt();
286 fClassifier->InsertExclAncestorID(pid);
287 }
288
289 // read ancestor IDs from configuration file and setup classifier
290
291 param = GetParam("IncludedAncestorIDs");
292 sizeParam = param.GetSize();
293
294 for(i = 0; i < sizeParam; ++i)
295 {
296 pid = param[i].GetInt();
297 fClassifier->InsertInclAncestorID(pid);
298 }
299
300 // create output arrays
301
302 fOutputArray = ExportArray("candidates");
303
304}
305
306//------------------------------------------------------------------------------
307
308void MadGraphShowerPartonSelector::Finish()
309{
310 if(fFilter) delete fFilter;
311 if(fClassifier) delete fClassifier;
312}
313
314//------------------------------------------------------------------------------
315
316void MadGraphShowerPartonSelector::Process()
317{
318 TObjArray *array = 0;
319 ExRootGenParticle *particle = 0;
320 ExRootCandidate *candidate = 0;
321 ExRootFactory *factory = GetFactory();
322
323 TLorentzVector momentum;
324
325 fFilter->Reset();
326 array = fFilter->GetSubArray(fClassifier, 0);
327
328 if(array == 0) return;
329
330 TIter itArray(array);
331
332 while(particle = static_cast<ExRootGenParticle*>(itArray.Next()))
333 {
334 momentum.SetPxPyPzE(particle->Px, particle->Py, particle->Pz, particle->E);
335
336 candidate = factory->NewCandidate();
337
338 candidate->SetP4(momentum);
339 candidate->SetType(particle->PID);
340
341 fOutputArray->Add(candidate);
342 }
343
344/*
345 cout << "==============================" << endl;
346 Int_t indexParticle = -1;
347 itArray.Reset();
348 while(particle = static_cast<ExRootGenParticle*>(itArray.Next()))
349 {
350 ++indexParticle;
351 cout << "--->\t" << particle->Status << "\t" << particle->PID << "\t";
352 cout << particle->M1 << "\t" << particle->M2 << "\t";
353 cout << particle->Px << "\t" << particle->Py << "\t" << particle->Pz << endl;
354 }
355*/
356}
357
358//------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.