Fork me on GitHub

source: svn/trunk/src/STDHEPConverter.cc@ 600

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

adding #include<cstdlib> for compatibility with newer versions of g++

File size: 7.9 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#include <cmath>
37#include <cstdlib>
38#include "TLorentzVector.h"
39#include "BlockClasses.h"
40#include "SmearUtil.h"
41#include "ExRootTreeWriter.h"
42#include "ExRootTreeBranch.h"
43#include "stdhep_mcfio.h"
44#include "stdhep_declarations.h"
45#include "STDHEPConverter.h"
46#include "PdgParticle.h"
47
48using namespace std;
49
50//---------------------------------------------------------------------------
51
52void STDHEPConverter::AnalyseEvent(ExRootTreeBranch *branch, const Long64_t eventNumber)
53{
54 TRootGenEvent *element;
55
56 element = static_cast<TRootGenEvent*>(branch->NewEntry());
57
58 element->Number = eventNumber;
59}
60
61//---------------------------------------------------------------------------
62
63void STDHEPConverter::AnalyseParticles(ExRootTreeBranch *branch)
64{
65 TRootC::GenParticle *element;
66
67 Double_t signPz;
68 TLorentzVector momentum;
69 Int_t number;
70
71 for(number = 0; number < myhepevt.nhep; ++number)
72 {
73 element = static_cast<TRootC::GenParticle*>(branch->NewEntry());
74
75 element->PID = myhepevt.idhep[number];
76 element->Status = myhepevt.isthep[number];
77 element->M1 = myhepevt.jmohep[number][0] - 1;
78 element->M2 = myhepevt.jmohep[number][1] - 1;
79 element->D1 = myhepevt.jdahep[number][0] - 1;
80 element->D2 = myhepevt.jdahep[number][1] - 1;
81 PdgParticle pdg_part = pdg_table[element->PID];
82 element->Charge = pdg_part.charge();
83 element->M = pdg_part.mass();
84
85 element->E = myhepevt.phep[number][3];
86 element->Px = myhepevt.phep[number][0];
87 element->Py = myhepevt.phep[number][1];
88 element->Pz = myhepevt.phep[number][2];
89 //element->M = myhepevt.phep[number][4];
90
91 momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
92 element->PT = momentum.Perp();
93 signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
94 //element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta(); to avoid a warning from ROOT, replace the "==0" by "< 1e-6"
95 element->Eta = element->PT <1e-6 ? signPz*999.9 : momentum.Eta();
96 element->Phi = momentum.Phi();
97
98 element->T = myhepevt.vhep[number][3];
99 element->X = myhepevt.vhep[number][0];
100 element->Y = myhepevt.vhep[number][1];
101 element->Z = myhepevt.vhep[number][2];
102 }
103}
104
105
106//------------------------------------------------------------------------------
107
108STDHEPConverter::~STDHEPConverter()
109{
110}
111
112//------------------------------------------------------------------------------
113
114STDHEPConverter::STDHEPConverter(const string& inputFileList, const string& outputFileName, const PdgTable& pdg, const int& Nevents) : DataConverter(pdg, Nevents)
115{
116 int ierr, entryType;
117 int istr = 0;
118 int nevt = 0, nevt_already_processed=0;
119
120 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
121
122 // information about generated event
123 ExRootTreeBranch *branchGenEvent = treeWriter->NewBranch("Event", TRootGenEvent::Class());
124 // generated particles from HEPEVT
125 ExRootTreeBranch *branchGenParticle = treeWriter->NewBranch("Particle", TRootC::GenParticle::Class());
126
127 // Open a stream connected to an event file:
128 ifstream infile(inputFileList.c_str());
129 string filename;
130 if(!infile.is_open())
131 {
132 string text = inputFileList + " for input ";
133 cerr << left << setw(26) <<"** ERROR: Can't open "<<""
134 << left << setw(40) << text <<""
135 << right << setw(3) <<" **"<<""<<endl;
136 exit(1);
137 }
138
139 // Parsing the input list
140 while(1)
141 {
142 infile >> filename;
143 if(!infile.good()) break; // end of the list
144 if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
145
146 ifstream checking_the_file(filename.c_str());
147 if(!checking_the_file.good()) {
148 cerr << left << setw(28) <<"** ERROR: Can't find file "<<""
149 << left << setw(38) << filename <<""
150 << right << setw(3) <<" **"<<""<<endl;
151 continue;
152 } else checking_the_file.close();
153
154 ierr = StdHepXdrReadInit(const_cast<char*>(filename.c_str()), &nevt, istr);
155
156 if(ierr != 0) {
157 cerr <<"** ERROR: Can't open file for input **"<<endl;
158 break;
159 }
160
161
162 Long64_t allEntries = (Nevt<0)?nevt:min(nevt,Nevt+1); // do not miss the "+1"
163 if(allEntries > 0)
164 {
165 // Loop over all objects
166 Long64_t entry = 0;
167 Long64_t recordNumber = 1;
168 for(entry = 0; entry < allEntries; ++entry)
169 {
170 ierr = StdHepXdrRead(&entryType, istr);
171 if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
172
173 if(ierr != 0) {
174 cerr << left << setw(49) <<"** ERROR: Unexpected end of file after"<<""
175 << left << setw(10) << entry <<""
176 << right << setw(10) <<"entries **"<<""<<endl;
177 break;
178 } // if error
179
180 // analyse only entries with standard HEPEVT common block
181 if(entryType == 1) {
182 // add empty events for missing event numbers
183 while(recordNumber < myhepevt.nevhep) {
184 treeWriter->Clear();
185 AnalyseEvent(branchGenEvent, recordNumber);
186 treeWriter->Fill();
187 ++recordNumber;
188 ++nevt_already_processed;
189 } // while
190
191 treeWriter->Clear();
192 AnalyseEvent(branchGenEvent, myhepevt.nevhep);
193 AnalyseParticles(branchGenParticle);
194 treeWriter->Fill();
195 ++recordNumber;
196 ++nevt_already_processed;
197 } // if
198 } // for
199 } // if allEntries >0
200 StdHepXdrEnd(istr);
201 } // while inputlist not empty
202 treeWriter->Write();
203
204 delete treeWriter;
205
206}
207
208
209
210
Note: See TracBrowser for help on using the repository browser.