Fork me on GitHub

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

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

bug fixed in the number of desired events to process

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 cerr << left << setw(30) <<"** ERROR: Can't open "<<""
132 << left << setw(20) << inputFileList <<""
133 << right << setw(19) <<"for input **"<<""<<endl;
134 exit(1);
135 }
136
137 // Parsing the input list
138 while(1)
139 {
140 infile >> filename;
141 if(!infile.good()) break; // end of the list
142 if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
143
144 ifstream checking_the_file(filename.c_str());
145 if(!checking_the_file.good()) {
146 cerr << left << setw(30) <<"** ERROR: Can't find file "<<""
147 << left << setw(20) << filename <<""
148 << right << setw(19) <<"for input **"<<""<<endl;
149 continue;
150 } else checking_the_file.close();
151
152 ierr = StdHepXdrReadInit(const_cast<char*>(filename.c_str()), &nevt, istr);
153
154 if(ierr != 0) {
155 cerr <<"** ERROR: Can't open file for input **"<<endl;
156 break;
157 }
158
159
160 Long64_t allEntries = (Nevt<0)?nevt:min(nevt,Nevt+1); // do not miss the "+1"
161 if(allEntries > 0)
162 {
163 // Loop over all objects
164 Long64_t entry = 0;
165 Long64_t recordNumber = 1;
166 for(entry = 0; entry < allEntries; ++entry)
167 {
168 ierr = StdHepXdrRead(&entryType, istr);
169 if (Nevt>0 && nevt_already_processed >=Nevt) break; // enough events already processed
170
171 if(ierr != 0) {
172 cerr << left << setw(49) <<"** ERROR: Unexpected end of file after"<<""
173 << left << setw(10) << entry <<""
174 << right << setw(10) <<"entries **"<<""<<endl;
175 break;
176 } // if error
177
178 // analyse only entries with standard HEPEVT common block
179 if(entryType == 1) {
180 // add empty events for missing event numbers
181 while(recordNumber < myhepevt.nevhep) {
182 treeWriter->Clear();
183 AnalyseEvent(branchGenEvent, recordNumber);
184 treeWriter->Fill();
185 ++recordNumber;
186 ++nevt_already_processed;
187 } // while
188
189 treeWriter->Clear();
190 AnalyseEvent(branchGenEvent, myhepevt.nevhep);
191 AnalyseParticles(branchGenParticle);
192 treeWriter->Fill();
193 ++recordNumber;
194 ++nevt_already_processed;
195 } // if
196 } // for
197 } // if allEntries >0
198 StdHepXdrEnd(istr);
199 } // while inputlist not empty
200 treeWriter->Write();
201
202 delete treeWriter;
203
204}
205
206
207
208
Note: See TracBrowser for help on using the repository browser.