Fork me on GitHub

source: git/classes/DelphesLHEFReader.cc @ 01f9722

ImprovedOutputFileTimingdual_readoutllp
Last change on this file since 01f9722 was 01f9722, checked in by Michele Selvaggi <michele.selvaggi@…>, 2 years ago

read xsec from LHE files produced with Wizard

  • Property mode set to 100644
File size: 7.4 KB
Line 
1/*
2 *  Delphes: a framework for fast simulation of a generic collider experiment
3 *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
4 *
5 *  This program is free software: you can redistribute it and/or modify
6 *  it under the terms of the GNU General Public License as published by
7 *  the Free Software Foundation, either version 3 of the License, or
8 *  (at your option) any later version.
9 *
10 *  This program is distributed in the hope that it will be useful,
11 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
12 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13 *  GNU General Public License for more details.
14 *
15 *  You should have received a copy of the GNU General Public License
16 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
17 */
18
19
20/** \class DelphesLHEFReader
21 *
22 *  Reads LHEF file
23 *
24 *  \author P. Demin - UCL, Louvain-la-Neuve
25 *
26 */
27
28#include "classes/DelphesLHEFReader.h"
29
30#include <stdexcept>
31#include <iostream>
32#include <sstream>
33#include <string>
34
35#include <stdio.h>
36
37#include "TObjArray.h"
38#include "TStopwatch.h"
39#include "TDatabasePDG.h"
40#include "TParticlePDG.h"
41#include "TLorentzVector.h"
42
43#include "classes/DelphesClasses.h"
44#include "classes/DelphesFactory.h"
45#include "classes/DelphesStream.h"
46
47#include "ExRootAnalysis/ExRootTreeBranch.h"
48
49using namespace std;
50
51static const int kBufferSize = 16384;
52
53//---------------------------------------------------------------------------
54
55DelphesLHEFReader::DelphesLHEFReader() :
56  fInputFile(0), fBuffer(0), fPDG(0),
57  fEventReady(kFALSE), fEventCounter(-1), fParticleCounter(-1),   fCrossSection(1)
58
59{
60  fBuffer = new char[kBufferSize];
61
62  fPDG = TDatabasePDG::Instance();
63}
64
65//---------------------------------------------------------------------------
66
67DelphesLHEFReader::~DelphesLHEFReader()
68{
69  if(fBuffer) delete[] fBuffer;
70}
71
72//---------------------------------------------------------------------------
73
74void DelphesLHEFReader::SetInputFile(FILE *inputFile)
75{
76  fInputFile = inputFile;
77}
78
79//---------------------------------------------------------------------------
80
81void DelphesLHEFReader::Clear()
82{
83  fEventReady = kFALSE;
84  fEventCounter = -1;
85  fParticleCounter = -1;
86  fWeightList.clear();
87}
88
89//---------------------------------------------------------------------------
90
91bool DelphesLHEFReader::EventReady()
92{
93  return fEventReady;
94}
95
96//---------------------------------------------------------------------------
97
98bool DelphesLHEFReader::ReadBlock(DelphesFactory *factory,
99  TObjArray *allParticleOutputArray,
100  TObjArray *stableParticleOutputArray,
101  TObjArray *partonOutputArray)
102{
103  int rc, id;
104  char *pch;
105  double weight;
106
107  if(!fgets(fBuffer, kBufferSize, fInputFile)) return kFALSE;
108
109  // read cross section string from LHE file header (meant to work on Wizard output only)
110  if(strstr(fBuffer, "<xsecinfo ") && fEventCounter < 0)
111  {
112    DelphesStream bufferStream(fBuffer);
113
114    string buf = fBuffer;   
115    string xsecstr="totxsec";
116   
117    if (buf.find(xsecstr) != string::npos)
118    {
119        buf = buf.substr(buf.find(xsecstr) + xsecstr.length()+1);
120        unsigned first = buf.find_first_of("\"")+1;
121        unsigned last = buf.find_last_of("\"");
122        xsecstr = buf.substr(first,last-first);
123        fCrossSection = stod(xsecstr);
124    } 
125  }
126 
127  if(strstr(fBuffer, "<event>"))
128  {
129    Clear();
130    fEventCounter = 1;
131  }
132 
133  else if(fEventCounter > 0)
134  {
135    DelphesStream bufferStream(fBuffer);
136
137    rc = bufferStream.ReadInt(fParticleCounter)
138      && bufferStream.ReadInt(fProcessID)
139      && bufferStream.ReadDbl(fWeight)
140      && bufferStream.ReadDbl(fScalePDF)
141      && bufferStream.ReadDbl(fAlphaQED)
142      && bufferStream.ReadDbl(fAlphaQCD);
143
144    if(!rc)
145    {
146      cerr << "** ERROR: " << "invalid event format" << endl;
147      return kFALSE;
148    }
149
150    --fEventCounter;
151  }
152  else if(fParticleCounter > 0)
153  {
154    DelphesStream bufferStream(fBuffer);
155
156    rc = bufferStream.ReadInt(fPID)
157      && bufferStream.ReadInt(fStatus)
158      && bufferStream.ReadInt(fM1)
159      && bufferStream.ReadInt(fM2)
160      && bufferStream.ReadInt(fC1)
161      && bufferStream.ReadInt(fC2)
162      && bufferStream.ReadDbl(fPx)
163      && bufferStream.ReadDbl(fPy)
164      && bufferStream.ReadDbl(fPz)
165      && bufferStream.ReadDbl(fE)
166      && bufferStream.ReadDbl(fMass);
167
168    if(!rc)
169    {
170      cerr << "** ERROR: " << "invalid particle format" << endl;
171      return kFALSE;
172    }
173
174    AnalyzeParticle(factory, allParticleOutputArray,
175      stableParticleOutputArray, partonOutputArray);
176
177    --fParticleCounter;
178  }
179  else if(strstr(fBuffer, "<wgt"))
180  {
181    pch = strpbrk(fBuffer, "\"'");
182    if(!pch)
183    {
184      cerr << "** ERROR: " << "invalid weight format" << endl;
185      return kFALSE;
186    }
187
188    DelphesStream idStream(pch + 1);
189    rc = idStream.ReadInt(id);
190
191    pch = strchr(fBuffer, '>');
192    if(!pch)
193    {
194      cerr << "** ERROR: " << "invalid weight format" << endl;
195      return kFALSE;
196    }
197
198    DelphesStream weightStream(pch + 1);
199    rc = weightStream.ReadDbl(weight);
200
201    if(!rc)
202    {
203      cerr << "** ERROR: " << "invalid weight format" << endl;
204      return kFALSE;
205    }
206
207    fWeightList.push_back(make_pair(id, weight));
208  }
209  else if(strstr(fBuffer, "</event>"))
210  {
211    fEventReady = kTRUE;
212  }
213
214  return kTRUE;
215}
216
217//---------------------------------------------------------------------------
218
219void DelphesLHEFReader::AnalyzeEvent(ExRootTreeBranch *branch, long long eventNumber,
220  TStopwatch *readStopWatch, TStopwatch *procStopWatch)
221{
222  LHEFEvent *element;
223
224  element = static_cast<LHEFEvent *>(branch->NewEntry());
225  element->Number = eventNumber;
226
227  element->ProcessID = fProcessID;
228  element->Weight = fWeight;
229  element->CrossSection = fCrossSection;
230
231  element->ScalePDF = fScalePDF;
232  element->AlphaQED = fAlphaQED;
233  element->AlphaQCD = fAlphaQCD;
234
235  element->ReadTime = readStopWatch->RealTime();
236  element->ProcTime = procStopWatch->RealTime();
237}
238
239//---------------------------------------------------------------------------
240
241void DelphesLHEFReader::AnalyzeWeight(ExRootTreeBranch *branch)
242{
243  LHEFWeight *element;
244  vector< pair< int, double > >::const_iterator itWeightList;
245
246  for(itWeightList = fWeightList.begin(); itWeightList != fWeightList.end(); ++itWeightList)
247  {
248    element = static_cast<LHEFWeight *>(branch->NewEntry());
249
250    element->ID = itWeightList->first;
251    element->Weight = itWeightList->second;
252  }
253}
254
255//---------------------------------------------------------------------------
256
257void DelphesLHEFReader::AnalyzeParticle(DelphesFactory *factory,
258  TObjArray *allParticleOutputArray,
259  TObjArray *stableParticleOutputArray,
260  TObjArray *partonOutputArray)
261{
262  Candidate *candidate;
263  TParticlePDG *pdgParticle;
264  int pdgCode;
265
266  candidate = factory->NewCandidate();
267
268  candidate->PID = fPID;
269  pdgCode = TMath::Abs(candidate->PID);
270
271  candidate->Status = fStatus;
272
273  pdgParticle = fPDG->GetParticle(fPID);
274  candidate->Charge = pdgParticle ? int(pdgParticle->Charge()/3.0) : -999;
275  candidate->Mass = fMass;
276
277  candidate->Momentum.SetPxPyPzE(fPx, fPy, fPz, fE);
278  candidate->Position.SetXYZT(0.0, 0.0, 0.0, 0.0);
279
280  candidate->M1 = fM1 - 1;
281  candidate->M2 = fM2 - 1;
282
283  candidate->D1 = -1;
284  candidate->D2 = -1;
285
286  allParticleOutputArray->Add(candidate);
287
288  if(!pdgParticle) return;
289
290  if(fStatus == 1)
291  {
292    stableParticleOutputArray->Add(candidate);
293  }
294  else if(pdgCode <= 5 || pdgCode == 21 || pdgCode == 15)
295  {
296    partonOutputArray->Add(candidate);
297  }
298}
299
300//---------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.