Fork me on GitHub

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

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

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