Fork me on GitHub

source: svn/trunk/classes/DelphesHepMCReader.cc@ 1105

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

add weight for HepMCEvent
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 9.2 KB
Line 
1
2/** \class DelphesHepMCReader
3 *
4 * Reads HepMC file
5 *
6 *
7 * $Date: 2013-05-01 17:29:33 +0000 (Wed, 01 May 2013) $
8 * $Revision: 1105 $
9 *
10 *
11 * \author P. Demin - UCL, Louvain-la-Neuve
12 *
13 */
14
15#include "classes/DelphesHepMCReader.h"
16
17#include <stdexcept>
18#include <iostream>
19#include <sstream>
20
21#include <map>
22#include <vector>
23
24#include <stdio.h>
25
26#include "TObjArray.h"
27#include "TStopwatch.h"
28#include "TDatabasePDG.h"
29#include "TParticlePDG.h"
30#include "TLorentzVector.h"
31
32#include "classes/DelphesClasses.h"
33#include "classes/DelphesFactory.h"
34#include "classes/DelphesStream.h"
35
36#include "ExRootAnalysis/ExRootTreeBranch.h"
37
38using namespace std;
39
40static const int kBufferSize = 1024;
41
42//---------------------------------------------------------------------------
43
44DelphesHepMCReader::DelphesHepMCReader() :
45 fInputFile(0), fBuffer(0), fPDG(0),
46 fVertexCounter(-1), fInCounter(-1), fOutCounter(-1),
47 fParticleCounter(0)
48{
49 fBuffer = new char[kBufferSize];
50
51 fPDG = TDatabasePDG::Instance();
52}
53
54//---------------------------------------------------------------------------
55
56DelphesHepMCReader::~DelphesHepMCReader()
57{
58 if(fBuffer) delete[] fBuffer;
59}
60
61//---------------------------------------------------------------------------
62
63void DelphesHepMCReader::SetInputFile(FILE *inputFile)
64{
65 fInputFile = inputFile;
66}
67
68//---------------------------------------------------------------------------
69
70void DelphesHepMCReader::Clear()
71{
72 fStateSize = 0;
73 fState.clear();
74 fWeightSize = 0;
75 fWeight.clear();
76 fVertexCounter = -1;
77 fInCounter = -1;
78 fOutCounter = -1;
79 fMotherMap.clear();
80 fDaughterMap.clear();
81 fParticleCounter = 0;
82}
83
84//---------------------------------------------------------------------------
85
86bool DelphesHepMCReader::EventReady()
87{
88 return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
89}
90
91//---------------------------------------------------------------------------
92
93bool DelphesHepMCReader::ReadBlock(DelphesFactory *factory,
94 TObjArray *allParticleOutputArray,
95 TObjArray *stableParticleOutputArray,
96 TObjArray *partonOutputArray)
97{
98 map< int, pair< int, int > >::iterator itMotherMap;
99 map< int, pair< int, int > >::iterator itDaughterMap;
100 char key;
101 int i, rc, state;
102 double weight;
103
104 if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
105
106 DelphesStream bufferStream(fBuffer + 1);
107
108 key = fBuffer[0];
109
110 if(key == 'E')
111 {
112 Clear();
113
114 rc = bufferStream.ReadInt(fEventNumber)
115 && bufferStream.ReadInt(fMPI)
116 && bufferStream.ReadDbl(fScale)
117 && bufferStream.ReadDbl(fAlphaQCD)
118 && bufferStream.ReadDbl(fAlphaQED)
119 && bufferStream.ReadInt(fProcessID)
120 && bufferStream.ReadInt(fSignalCode)
121 && bufferStream.ReadInt(fVertexCounter)
122 && bufferStream.ReadInt(fBeamCode[0])
123 && bufferStream.ReadInt(fBeamCode[1])
124 && bufferStream.ReadInt(fStateSize);
125
126 if(!rc)
127 {
128 cerr << "** ERROR: " << "invalid event format" << endl;
129 return kFALSE;
130 }
131
132 for(i = 0; i < fStateSize; ++i)
133 {
134 rc = rc && bufferStream.ReadInt(state);
135 fState.push_back(state);
136 }
137
138 rc = rc && bufferStream.ReadInt(fWeightSize);
139
140 if(!rc)
141 {
142 cerr << "** ERROR: " << "invalid event format" << endl;
143 return kFALSE;
144 }
145
146 for(i = 0; i < fWeightSize; ++i)
147 {
148 rc = rc && bufferStream.ReadDbl(weight);
149 fWeight.push_back(weight);
150 }
151
152 if(!rc)
153 {
154 cerr << "** ERROR: " << "invalid event format" << endl;
155 return kFALSE;
156 }
157 }
158 else if(key == 'F')
159 {
160 rc = bufferStream.ReadInt(fID1)
161 && bufferStream.ReadInt(fID2)
162 && bufferStream.ReadDbl(fX1)
163 && bufferStream.ReadDbl(fX2)
164 && bufferStream.ReadDbl(fScalePDF)
165 && bufferStream.ReadDbl(fPDF1)
166 && bufferStream.ReadDbl(fPDF2);
167
168 if(!rc)
169 {
170 cerr << "** ERROR: " << "invalid PDF format" << endl;
171 return kFALSE;
172 }
173 }
174 else if(key == 'V' && fVertexCounter > 0)
175 {
176 rc = bufferStream.ReadInt(fOutVertexCode)
177 && bufferStream.ReadInt(fVertexID)
178 && bufferStream.ReadDbl(fX)
179 && bufferStream.ReadDbl(fY)
180 && bufferStream.ReadDbl(fZ)
181 && bufferStream.ReadDbl(fT)
182 && bufferStream.ReadInt(fInCounter)
183 && bufferStream.ReadInt(fOutCounter);
184
185 if(!rc)
186 {
187 cerr << "** ERROR: " << "invalid vertex format" << endl;
188 return kFALSE;
189 }
190 --fVertexCounter;
191 }
192 else if(key == 'P' && fOutCounter > 0)
193 {
194 rc = bufferStream.ReadInt(fParticleCode)
195 && bufferStream.ReadInt(fPID)
196 && bufferStream.ReadDbl(fPx)
197 && bufferStream.ReadDbl(fPy)
198 && bufferStream.ReadDbl(fPz)
199 && bufferStream.ReadDbl(fE)
200 && bufferStream.ReadDbl(fMass)
201 && bufferStream.ReadInt(fStatus)
202 && bufferStream.ReadDbl(fTheta)
203 && bufferStream.ReadDbl(fPhi)
204 && bufferStream.ReadInt(fInVertexCode);
205
206 if(!rc)
207 {
208 cerr << "** ERROR: " << "invalid particle format" << endl;
209 return kFALSE;
210 }
211
212 if(fInVertexCode < 0)
213 {
214 itMotherMap = fMotherMap.find(fInVertexCode);
215 if(itMotherMap == fMotherMap.end())
216 {
217 fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
218 }
219 else
220 {
221 itMotherMap->second.second = fParticleCounter;
222 }
223 }
224
225 if(fInCounter <= 0)
226 {
227 itDaughterMap = fDaughterMap.find(fOutVertexCode);
228 if(itDaughterMap == fDaughterMap.end())
229 {
230 fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
231 }
232 else
233 {
234 itDaughterMap->second.second = fParticleCounter;
235 }
236 }
237
238 AnalyzeParticle(factory, allParticleOutputArray,
239 stableParticleOutputArray, partonOutputArray);
240
241 if(fInCounter > 0)
242 {
243 --fInCounter;
244 }
245 else
246 {
247 --fOutCounter;
248 }
249
250 ++fParticleCounter;
251 }
252
253 if(EventReady())
254 {
255 FinalizeParticles(allParticleOutputArray);
256 }
257
258 return kTRUE;
259}
260
261//---------------------------------------------------------------------------
262
263void DelphesHepMCReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
264 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
265{
266 HepMCEvent *element;
267
268 element = static_cast<HepMCEvent *>(branch->NewEntry());
269 element->Number = fEventNumber;
270
271 element->ProcessID = fProcessID;
272 element->MPI = fMPI;
273 element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
274 element->Scale = fScale;
275 element->AlphaQED = fAlphaQED;
276 element->AlphaQCD = fAlphaQCD;
277
278 element->ID1 = fID1;
279 element->ID2 = fID2;
280 element->X1 = fX1;
281 element->X2 = fX2;
282 element->ScalePDF = fScalePDF;
283 element->PDF1 = fPDF1;
284 element->PDF2 = fPDF2;
285
286 element->ReadTime = readStopWatch->RealTime();
287 element->ProcTime = procStopWatch->RealTime();
288}
289
290//---------------------------------------------------------------------------
291
292void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
293 TObjArray *allParticleOutputArray,
294 TObjArray *stableParticleOutputArray,
295 TObjArray *partonOutputArray)
296{
297 Candidate *candidate;
298 TParticlePDG *pdgParticle;
299 int pdgCode;
300
301 candidate = factory->NewCandidate();
302
303 candidate->PID = fPID;
304 pdgCode = TMath::Abs(candidate->PID);
305
306 candidate->Status = fStatus;
307
308 pdgParticle = fPDG->GetParticle(fPID);
309 candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
310 candidate->Mass = pdgParticle ? pdgParticle->Mass() : -999.9;
311
312 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
313
314 candidate->M2 = 1;
315 candidate->D2 = 1;
316 if(fInCounter > 0)
317 {
318 candidate->M1 = 1;
319 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
320 }
321 else
322 {
323 candidate->M1 = fOutVertexCode;
324 candidate->Position.SetXYZT(fX, fY, fZ, fT);
325 }
326 if(fInVertexCode < 0)
327 {
328 candidate->D1 = fInVertexCode;
329 }
330 else
331 {
332 candidate->D1 = 1;
333 }
334
335 allParticleOutputArray->Add(candidate);
336
337 if(!pdgParticle) return;
338
339 if(fStatus == 1 && pdgParticle->Stable())
340 {
341 stableParticleOutputArray->Add(candidate);
342 }
343 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
344 {
345 partonOutputArray->Add(candidate);
346 }
347}
348
349//---------------------------------------------------------------------------
350
351void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
352{
353 Candidate *candidate;
354 map< int, pair< int, int > >::iterator itMotherMap;
355 map< int, pair< int, int > >::iterator itDaughterMap;
356 int i;
357
358 for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
359 {
360 candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
361
362 if(candidate->M1 > 0)
363 {
364 candidate->M1 = -1;
365 candidate->M2 = -1;
366 }
367 else
368 {
369 itMotherMap = fMotherMap.find(candidate->M1);
370 if(itMotherMap == fMotherMap.end())
371 {
372 candidate->M1 = -1;
373 candidate->M2 = -1;
374 }
375 else
376 {
377 candidate->M1 = itMotherMap->second.first;
378 candidate->M2 = itMotherMap->second.second;
379 }
380 }
381 if(candidate->D1 > 0)
382 {
383 candidate->D1 = -1;
384 candidate->D2 = -1;
385 }
386 else
387 {
388 itDaughterMap = fDaughterMap.find(candidate->D1);
389 if(itDaughterMap == fDaughterMap.end())
390 {
391 candidate->D1 = -1;
392 candidate->D2 = -1;
393 }
394 else
395 {
396 candidate->D1 = itDaughterMap->second.first;
397 candidate->D2 = itDaughterMap->second.second;
398 }
399 }
400 }
401}
402
403//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.