Fork me on GitHub

source: svn/trunk/src/LHEFConverter.cc@ 507

Last change on this file since 507 was 480, checked in by Xavier Rouby, 15 years ago

minor change: cleaning of the 'cout' and 'cerr'
File size: 9.0 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \------------- arXiv:0903.2225v1 ------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** ROOT: Nucl. Inst. & Meth. in Phys. Res. A389 (1997) 81-86 **
12** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
13** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
14** FROG: [hep-ex/0901.2718v1] **
15** HepMC: Comput. Phys. Commun.134 (2001) 41 **
16** **
17** ------------------------------------------------------------------ **
18** **
19** Main authors: **
20** ------------- **
21** **
22** Severine Ovyn Xavier Rouby **
23** severine.ovyn@uclouvain.be xavier.rouby@cern **
24** **
25** Center for Particle Physics and Phenomenology (CP3) **
26** Universite catholique de Louvain (UCL) **
27** Louvain-la-Neuve, Belgium **
28** **
29** Copyright (C) 2008-2009, **
30** All rights reserved. **
31** **
32***********************************************************************/
33
34#include <iostream>
35#include <fstream>
36
37#include "TLorentzVector.h"
38
39#include "ExRootTreeWriter.h"
40#include "BlockClasses.h"
41#include "LHEFConverter.h"
42#include "SmearUtil.h"
43#include "PdgParticle.h"
44
45using namespace std;
46
47
48//------------------------------------------------------------------------------
49
50void LHEFConverter::AnalyseEvent(LHEF::Reader *reader, ExRootTreeBranch *branch, const Long64_t eventNumber)
51{
52 const LHEF::HEPEUP &hepeup = reader->hepeup;
53
54 TRootLHEFEvent *element;
55
56 element = (TRootLHEFEvent*) branch->NewEntry();
57
58 element->Number = eventNumber;
59 element->Nparticles = hepeup.NUP;
60 element->ProcessID = hepeup.IDPRUP;
61 element->Weight = hepeup.XWGTUP;
62 element->ScalePDF = hepeup.SCALUP;
63 element->CouplingQED = hepeup.AQEDUP;
64 element->CouplingQCD = hepeup.AQCDUP;
65}
66
67//---------------------------------------------------------------------------
68
69void LHEFConverter::AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
70{
71 const LHEF::HEPEUP &hepeup = reader->hepeup;
72 Double_t signPz;
73 TLorentzVector momentum;
74 TRootLHEFParticle *element;
75
76 for(Int_t particle = 0; particle < hepeup.NUP; ++particle)
77 {
78 element = (TRootLHEFParticle*) branch->NewEntry();
79
80 element->PID = hepeup.IDUP[particle];
81 element->Status = hepeup.ISTUP[particle];
82 element->M1 = hepeup.MOTHUP[particle].first;
83 element->M2 = hepeup.MOTHUP[particle].second;
84 element->D1 = 0;
85 element->D2 = 0;
86 element->ColorLine1 = hepeup.ICOLUP[particle].first;
87 element->ColorLine2 = hepeup.ICOLUP[particle].second;
88 element->Px = hepeup.PUP[particle][0];
89 element->Py = hepeup.PUP[particle][1];
90 element->Pz = hepeup.PUP[particle][2];
91 element->E = hepeup.PUP[particle][3];
92 //element->M = hepeup.PUP[particle][4];
93 //element->Charge = ChargeVal(element->PID);
94
95 PdgParticle pdg_part = pdg_table[element->PID];
96 element->Charge = pdg_part.charge();
97 element->M = pdg_part.mass();
98
99 momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
100 element->PT = momentum.Perp();
101 signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
102 //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
103 element->Eta = element->PT < 1E-6 ? signPz*999.9 : momentum.Eta();
104 element->Phi = momentum.Phi();
105 element->Rapidity = element->PT < 1E-6 ? signPz*999.9 : momentum.Rapidity();
106 element->LifeTime = hepeup.VTIMUP[particle];
107 element->Spin = hepeup.SPINUP[particle];
108 }
109}
110
111
112//---------------------------------------------------------------------------
113/*
114void LHEFConverter::AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
115{
116 const LHEF::HEPEUP &hepeup = reader->hepeup;
117 Double_t signPz;
118 TLorentzVector momentum;
119 TRootC::GenParticle *element;
120
121 for(Int_t particle = 0; particle < hepeup.NUP; ++particle)
122 {
123 element = (TRootC::GenParticle*) branch->NewEntry();
124
125 element->PID = hepeup.IDUP[particle];
126 element->Status = hepeup.ISTUP[particle];
127 element->M1 = hepeup.MOTHUP[particle].first;
128 element->M2 = hepeup.MOTHUP[particle].second;
129 //element->ColorLine1 = hepeup.ICOLUP[particle].first;
130 //element->ColorLine2 = hepeup.ICOLUP[particle].second;
131 element->D1 = 0;
132 element->D2 = 0;
133
134 PdgParticle pdg_part = pdg_table[element->PID];
135 element->Charge = pdg_part.charge();
136 element->M = pdg_part.mass();
137 //element->Charge = ChargeVal(element->PID);
138 //element->M = hepeup.PUP[particle][4];
139
140 element->E = hepeup.PUP[particle][3];
141 element->Px = hepeup.PUP[particle][0];
142 element->Py = hepeup.PUP[particle][1];
143 element->Pz = hepeup.PUP[particle][2];
144
145 momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
146 element->PT = momentum.Perp();
147 signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
148 //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
149 element->Eta = element->PT < 1E-6 ? signPz*999.9 : momentum.Eta();
150 element->Phi = momentum.Phi();
151// element->Rapidity = element->PT < 1E-6 ? signPz*999.9 : momentum.Rapidity();
152// element->LifeTime = hepeup.VTIMUP[particle];
153// element->Spin = hepeup.SPINUP[particle];
154 }
155}
156*/
157
158LHEFConverter::~LHEFConverter()
159{
160}
161
162//------------------------------------------------------------------------------
163LHEFConverter::LHEFConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg,Nevents) {
164
165 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
166
167 // generated event from LHEF
168 ExRootTreeBranch *branchEvent = treeWriter->NewBranch("Event", TRootLHEFEvent::Class());
169
170 // generated partons from LHEF
171 ExRootTreeBranch *branchParticle = treeWriter->NewBranch("Particle", TRootLHEFParticle::Class());
172 //ExRootTreeBranch *branchParticle = treeWriter->NewBranch("Particle", TRootC::GenParticle::Class());
173
174 int nevt_already_processed=0;
175
176 // Open a stream connected to an event file:
177 ifstream infile(inputFileList.c_str());
178 string filename;
179 if(!infile.is_open()) {
180 cerr << left << setw(30) <<"** ERROR: Can't open "<<""
181 << left << setw(20) << inputFileList <<""
182 << right << setw(19) <<"for input **"<<""<<endl;
183 exit(1);
184 }
185
186 while(1) { // parsing the list of files
187
188 infile >> filename;
189 if(!infile.good()) break;
190 if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
191
192 ifstream checking_the_file(filename.c_str());
193 if(!checking_the_file.good())
194 {
195 cerr << left << setw(30) <<"** ERROR: Can't find file "<<""
196 << left << setw(36) << filename <<""
197 << right << setw(3) <<" **"<<""<<endl;
198 continue;
199 }
200 else checking_the_file.close();
201
202 // Create the Reader object:
203 LHEF::Reader *inputReader = new LHEF::Reader(filename);
204
205 Long64_t allEntries = inputReader->getNumberOfEvents();
206 allEntries = (Nevt<0)?allEntries:min((int)allEntries,Nevt); // do not miss the "+1"
207
208 if(allEntries > 0) {
209 // Loop over all events
210 Long64_t entry = 0;
211 while(inputReader->readEvent()) {
212 if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
213 treeWriter->Clear();
214 AnalyseEvent(inputReader, branchEvent, entry + 1);
215 AnalyseParticles(inputReader, branchParticle);
216 treeWriter->Fill();
217 ++entry;
218 ++nevt_already_processed;
219 if(allEntries<entry+1) break;
220 }
221 }
222 }
223 treeWriter->Write();
224 delete treeWriter;
225 //delete inputReader;
226}
227
Note: See TracBrowser for help on using the repository browser.