Fork me on GitHub

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

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

add missing fDaughterMap.clear()

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