Fork me on GitHub

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

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

use mass from input file
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 10.0 KB
Line 
1
2/** \class DelphesHepMCReader
3 *
4 * Reads HepMC file
5 *
6 *
7 * $Date: 2013-07-08 19:27:21 +0000 (Mon, 08 Jul 2013) $
8 * $Revision: 1189 $
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 fMomentumCoefficient = 1.0;
77 fPositionCoefficient = 1.0;
78 fVertexCounter = -1;
79 fInCounter = -1;
80 fOutCounter = -1;
81 fMotherMap.clear();
82 fDaughterMap.clear();
83 fParticleCounter = 0;
84}
85
86//---------------------------------------------------------------------------
87
88bool DelphesHepMCReader::EventReady()
89{
90 return (fVertexCounter == 0) && (fInCounter == 0) && (fOutCounter == 0);
91}
92
93//---------------------------------------------------------------------------
94
95bool DelphesHepMCReader::ReadBlock(DelphesFactory *factory,
96 TObjArray *allParticleOutputArray,
97 TObjArray *stableParticleOutputArray,
98 TObjArray *partonOutputArray)
99{
100 map< int, pair< int, int > >::iterator itMotherMap;
101 map< int, pair< int, int > >::iterator itDaughterMap;
102 char key, momentumUnit[4], positionUnit[3];
103 int i, rc, state;
104 double weight;
105
106 if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
107
108 DelphesStream bufferStream(fBuffer + 1);
109
110 key = fBuffer[0];
111
112 if(key == 'E')
113 {
114 Clear();
115
116 rc = bufferStream.ReadInt(fEventNumber)
117 && bufferStream.ReadInt(fMPI)
118 && bufferStream.ReadDbl(fScale)
119 && bufferStream.ReadDbl(fAlphaQCD)
120 && bufferStream.ReadDbl(fAlphaQED)
121 && bufferStream.ReadInt(fProcessID)
122 && bufferStream.ReadInt(fSignalCode)
123 && bufferStream.ReadInt(fVertexCounter)
124 && bufferStream.ReadInt(fBeamCode[0])
125 && bufferStream.ReadInt(fBeamCode[1])
126 && bufferStream.ReadInt(fStateSize);
127
128 if(!rc)
129 {
130 cerr << "** ERROR: " << "invalid event format" << endl;
131 return kFALSE;
132 }
133
134 for(i = 0; i < fStateSize; ++i)
135 {
136 rc = rc && bufferStream.ReadInt(state);
137 fState.push_back(state);
138 }
139
140 rc = rc && bufferStream.ReadInt(fWeightSize);
141
142 if(!rc)
143 {
144 cerr << "** ERROR: " << "invalid event format" << endl;
145 return kFALSE;
146 }
147
148 for(i = 0; i < fWeightSize; ++i)
149 {
150 rc = rc && bufferStream.ReadDbl(weight);
151 fWeight.push_back(weight);
152 }
153
154 if(!rc)
155 {
156 cerr << "** ERROR: " << "invalid event format" << endl;
157 return kFALSE;
158 }
159 }
160 else if(key == 'U')
161 {
162 rc = sscanf(fBuffer + 1, "%3s %2s", momentumUnit, positionUnit);
163
164 if(rc != 2)
165 {
166 cerr << "** ERROR: " << "invalid units format" << endl;
167 return kFALSE;
168 }
169
170 if(strncmp(momentumUnit, "GEV", 3) == 0)
171 {
172 fMomentumCoefficient = 1.0;
173 }
174 else if(strncmp(momentumUnit, "MEV", 3) == 0)
175 {
176 fMomentumCoefficient = 0.001;
177 }
178
179 if(strncmp(positionUnit, "MM", 3) == 0)
180 {
181 fPositionCoefficient = 1.0;
182 }
183 else if(strncmp(positionUnit, "CM", 3) == 0)
184 {
185 fPositionCoefficient = 10.0;
186 }
187 }
188 else if(key == 'F')
189 {
190 rc = bufferStream.ReadInt(fID1)
191 && bufferStream.ReadInt(fID2)
192 && bufferStream.ReadDbl(fX1)
193 && bufferStream.ReadDbl(fX2)
194 && bufferStream.ReadDbl(fScalePDF)
195 && bufferStream.ReadDbl(fPDF1)
196 && bufferStream.ReadDbl(fPDF2);
197
198 if(!rc)
199 {
200 cerr << "** ERROR: " << "invalid PDF format" << endl;
201 return kFALSE;
202 }
203 }
204 else if(key == 'V' && fVertexCounter > 0)
205 {
206 rc = bufferStream.ReadInt(fOutVertexCode)
207 && bufferStream.ReadInt(fVertexID)
208 && bufferStream.ReadDbl(fX)
209 && bufferStream.ReadDbl(fY)
210 && bufferStream.ReadDbl(fZ)
211 && bufferStream.ReadDbl(fT)
212 && bufferStream.ReadInt(fInCounter)
213 && bufferStream.ReadInt(fOutCounter);
214
215 if(!rc)
216 {
217 cerr << "** ERROR: " << "invalid vertex format" << endl;
218 return kFALSE;
219 }
220 --fVertexCounter;
221 }
222 else if(key == 'P' && fOutCounter > 0)
223 {
224 rc = bufferStream.ReadInt(fParticleCode)
225 && bufferStream.ReadInt(fPID)
226 && bufferStream.ReadDbl(fPx)
227 && bufferStream.ReadDbl(fPy)
228 && bufferStream.ReadDbl(fPz)
229 && bufferStream.ReadDbl(fE)
230 && bufferStream.ReadDbl(fMass)
231 && bufferStream.ReadInt(fStatus)
232 && bufferStream.ReadDbl(fTheta)
233 && bufferStream.ReadDbl(fPhi)
234 && bufferStream.ReadInt(fInVertexCode);
235
236 if(!rc)
237 {
238 cerr << "** ERROR: " << "invalid particle format" << endl;
239 return kFALSE;
240 }
241
242 if(fInVertexCode < 0)
243 {
244 itMotherMap = fMotherMap.find(fInVertexCode);
245 if(itMotherMap == fMotherMap.end())
246 {
247 fMotherMap[fInVertexCode] = make_pair(fParticleCounter, -1);
248 }
249 else
250 {
251 itMotherMap->second.second = fParticleCounter;
252 }
253 }
254
255 if(fInCounter <= 0)
256 {
257 itDaughterMap = fDaughterMap.find(fOutVertexCode);
258 if(itDaughterMap == fDaughterMap.end())
259 {
260 fDaughterMap[fOutVertexCode] = make_pair(fParticleCounter, fParticleCounter);
261 }
262 else
263 {
264 itDaughterMap->second.second = fParticleCounter;
265 }
266 }
267
268 AnalyzeParticle(factory, allParticleOutputArray,
269 stableParticleOutputArray, partonOutputArray);
270
271 if(fInCounter > 0)
272 {
273 --fInCounter;
274 }
275 else
276 {
277 --fOutCounter;
278 }
279
280 ++fParticleCounter;
281 }
282
283 if(EventReady())
284 {
285 FinalizeParticles(allParticleOutputArray);
286 }
287
288 return kTRUE;
289}
290
291//---------------------------------------------------------------------------
292
293void DelphesHepMCReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
294 TStopwatch *readStopWatch, TStopwatch *procStopWatch)
295{
296 HepMCEvent *element;
297
298 element = static_cast<HepMCEvent *>(branch->NewEntry());
299 element->Number = fEventNumber;
300
301 element->ProcessID = fProcessID;
302 element->MPI = fMPI;
303 element->Weight = fWeight.size() > 0 ? fWeight[0] : 1.0;
304 element->Scale = fScale;
305 element->AlphaQED = fAlphaQED;
306 element->AlphaQCD = fAlphaQCD;
307
308 element->ID1 = fID1;
309 element->ID2 = fID2;
310 element->X1 = fX1;
311 element->X2 = fX2;
312 element->ScalePDF = fScalePDF;
313 element->PDF1 = fPDF1;
314 element->PDF2 = fPDF2;
315
316 element->ReadTime = readStopWatch->RealTime();
317 element->ProcTime = procStopWatch->RealTime();
318}
319
320//---------------------------------------------------------------------------
321
322void DelphesHepMCReader::AnalyzeParticle(DelphesFactory *factory,
323 TObjArray *allParticleOutputArray,
324 TObjArray *stableParticleOutputArray,
325 TObjArray *partonOutputArray)
326{
327 Candidate *candidate;
328 TParticlePDG *pdgParticle;
329 int pdgCode;
330
331 candidate = factory->NewCandidate();
332
333 candidate->PID = fPID;
334 pdgCode = TMath::Abs(candidate->PID);
335
336 candidate->Status = fStatus;
337
338 pdgParticle = fPDG->GetParticle(fPID);
339 candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
340 candidate->Mass = fMass;
341
342 candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
343 if(fMomentumCoefficient != 1.0)
344 {
345 candidate->Momentum *= fMomentumCoefficient;
346 }
347
348 candidate->M2 = 1;
349 candidate->D2 = 1;
350 if(fInCounter > 0)
351 {
352 candidate->M1 = 1;
353 candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
354 }
355 else
356 {
357 candidate->M1 = fOutVertexCode;
358 candidate->Position.SetXYZT(fX, fY, fZ, fT);
359 if(fPositionCoefficient != 1.0)
360 {
361 candidate->Position *= fPositionCoefficient;
362 }
363 }
364 if(fInVertexCode < 0)
365 {
366 candidate->D1 = fInVertexCode;
367 }
368 else
369 {
370 candidate->D1 = 1;
371 }
372
373 allParticleOutputArray->Add(candidate);
374
375 if(!pdgParticle) return;
376
377 if(fStatus == 1 && pdgParticle->Stable())
378 {
379 stableParticleOutputArray->Add(candidate);
380 }
381 else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
382 {
383 partonOutputArray->Add(candidate);
384 }
385}
386
387//---------------------------------------------------------------------------
388
389void DelphesHepMCReader::FinalizeParticles(TObjArray *allParticleOutputArray)
390{
391 Candidate *candidate;
392 map< int, pair< int, int > >::iterator itMotherMap;
393 map< int, pair< int, int > >::iterator itDaughterMap;
394 int i;
395
396 for(i = 0; i < allParticleOutputArray->GetEntriesFast(); ++i)
397 {
398 candidate = static_cast<Candidate *>(allParticleOutputArray->At(i));
399
400 if(candidate->M1 > 0)
401 {
402 candidate->M1 = -1;
403 candidate->M2 = -1;
404 }
405 else
406 {
407 itMotherMap = fMotherMap.find(candidate->M1);
408 if(itMotherMap == fMotherMap.end())
409 {
410 candidate->M1 = -1;
411 candidate->M2 = -1;
412 }
413 else
414 {
415 candidate->M1 = itMotherMap->second.first;
416 candidate->M2 = itMotherMap->second.second;
417 }
418 }
419 if(candidate->D1 > 0)
420 {
421 candidate->D1 = -1;
422 candidate->D2 = -1;
423 }
424 else
425 {
426 itDaughterMap = fDaughterMap.find(candidate->D1);
427 if(itDaughterMap == fDaughterMap.end())
428 {
429 candidate->D1 = -1;
430 candidate->D2 = -1;
431 }
432 else
433 {
434 candidate->D1 = itDaughterMap->second.first;
435 candidate->D2 = itDaughterMap->second.second;
436 }
437 }
438 }
439}
440
441//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.