Fork me on GitHub

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

Last change on this file since 376 was 369, checked in by Xavier Rouby, 16 years ago

avoid the ROOT warning when PT is very small
File size: 6.4 KB
Line 
1/***********************************************************************
2** **
3** /----------------------------------------------\ **
4** | Delphes, a framework for the fast simulation | **
5** | of a generic collider experiment | **
6** \----------------------------------------------/ **
7** **
8** **
9** This package uses: **
10** ------------------ **
11** FastJet algorithm: Phys. Lett. B641 (2006) [hep-ph/0512210] **
12** Hector: JINST 2:P09005 (2007) [physics.acc-ph:0707.1198v2] **
13** FROG: [hep-ex/0901.2718v1] **
14** **
15** ------------------------------------------------------------------ **
16** **
17** Main authors: **
18** ------------- **
19** **
20** Severine Ovyn Xavier Rouby **
21** severine.ovyn@uclouvain.be xavier.rouby@cern **
22** **
23** Center for Particle Physics and Phenomenology (CP3) **
24** Universite catholique de Louvain (UCL) **
25** Louvain-la-Neuve, Belgium **
26** **
27** Copyright (C) 2008-2009, **
28** All rights reserved. **
29** **
30***********************************************************************/
31
32
33#include <iostream>
34#include <fstream>
35
36#include "TLorentzVector.h"
37
38#include "ExRootTreeWriter.h"
39#include "BlockClasses.h"
40#include "LHEFConverter.h"
41
42using namespace std;
43
44
45//------------------------------------------------------------------------------
46
47void LHEFConverter::AnalyseEvent(LHEF::Reader *reader, ExRootTreeBranch *branch, const Long64_t eventNumber)
48{
49 const LHEF::HEPEUP &hepeup = reader->hepeup;
50
51 TRootLHEFEvent *element;
52
53 element = (TRootLHEFEvent*) branch->NewEntry();
54
55 element->Number = eventNumber;
56 element->Nparticles = hepeup.NUP;
57 element->ProcessID = hepeup.IDPRUP;
58 element->Weight = hepeup.XWGTUP;
59 element->ScalePDF = hepeup.SCALUP;
60 element->CouplingQED = hepeup.AQEDUP;
61 element->CouplingQCD = hepeup.AQCDUP;
62}
63
64//---------------------------------------------------------------------------
65
66void LHEFConverter::AnalyseParticles(LHEF::Reader *reader, ExRootTreeBranch *branch)
67{
68 const LHEF::HEPEUP &hepeup = reader->hepeup;
69 Double_t signPz;
70 TLorentzVector momentum;
71 TRootLHEFParticle *element;
72
73 for(Int_t particle = 0; particle < hepeup.NUP; ++particle)
74 {
75 element = (TRootLHEFParticle*) branch->NewEntry();
76
77 element->PID = hepeup.IDUP[particle];
78 element->Status = hepeup.ISTUP[particle];
79 element->Mother1 = hepeup.MOTHUP[particle].first;
80 element->Mother2 = hepeup.MOTHUP[particle].second;
81 element->ColorLine1 = hepeup.ICOLUP[particle].first;
82 element->ColorLine2 = hepeup.ICOLUP[particle].second;
83 element->Px = hepeup.PUP[particle][0];
84 element->Py = hepeup.PUP[particle][1];
85 element->Pz = hepeup.PUP[particle][2];
86 element->E = hepeup.PUP[particle][3];
87 element->M = hepeup.PUP[particle][4];
88
89 momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
90 element->PT = momentum.Perp();
91 signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
92 //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
93 element->Eta = element->PT < 1e-6 ? signPz*999.9 : momentum.Eta();
94 element->Phi = momentum.Phi();
95 element->Rapidity = element->PT == 0.0 ? signPz*999.9 : momentum.Rapidity();
96 element->LifeTime = hepeup.VTIMUP[particle];
97 element->Spin = hepeup.SPINUP[particle];
98 }
99}
100
101LHEFConverter::~LHEFConverter()
102{
103}
104
105//------------------------------------------------------------------------------
106//Nevents not yet implemented! 08.03.2009
107LHEFConverter::LHEFConverter(const string& inputFileList, const string& outputFileName, const int& Nevents) :
108 DataConverter(Nevents) {
109
110 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
111
112 // generated event from LHEF
113 ExRootTreeBranch *branchEvent = treeWriter->NewBranch("Event", TRootLHEFEvent::Class());
114
115 // generated partons from LHEF
116 ExRootTreeBranch *branchParticle = treeWriter->NewBranch("Particle", TRootLHEFParticle::Class());
117
118 // Open a stream connected to an event file:
119 ifstream infile(inputFileList.c_str());
120 string filename;
121 if(!infile.is_open()) {
122 cerr << left << setw(30) <<"** ERROR: Can't open "<<""
123 << left << setw(20) << inputFileList <<""
124 << right << setw(19) <<"for input **"<<""<<endl;
125 exit(1);
126 }
127
128 while(1) { // parsing the list of files
129
130 infile >> filename;
131 if(!infile.good()) break;
132 ifstream checking_the_file(filename.c_str());
133 if(!checking_the_file.good())
134 {
135 cerr << left << setw(30) <<"** ERROR: Can't find file "<<""
136 << left << setw(20) << filename <<""
137 << right << setw(19) <<"for input **"<<""<<endl;
138 continue;
139 }
140 else checking_the_file.close();
141
142 // Create the Reader object:
143 LHEF::Reader *inputReader = new LHEF::Reader(filename);
144
145 Long64_t allEntries = inputReader->getNumberOfEvents();
146
147 if(allEntries > 0) {
148 // Loop over all events
149 Long64_t entry = 0;
150 while(inputReader->readEvent()) {
151 treeWriter->Clear();
152 AnalyseEvent(inputReader, branchEvent, entry + 1);
153 AnalyseParticles(inputReader, branchParticle);
154 treeWriter->Fill();
155 ++entry;
156 }
157 }
158 }
159 treeWriter->Write();
160
161
162 delete treeWriter;
163 //delete inputReader;
164}
165
Note: See TracBrowser for help on using the repository browser.