Fork me on GitHub

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

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

minor change: cleaning of the 'cout' and 'cerr'

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