Fork me on GitHub

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

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

first commit

File size: 4.6 KB
Line 
1#include <iostream>
2#include <utility>
3#include <deque>
4#include <sstream>
5#include <fstream>
6
7#include "TROOT.h"
8#include "TApplication.h"
9#include "TFile.h"
10#include "TTree.h"
11#include "TBranch.h"
12#include "TLeaf.h"
13#include "TLorentzVector.h"
14#include "TChain.h"
15
16#include "Utilities/ExRootAnalysis/interface/BlockClasses.h"
17#include "Utilities/ExRootAnalysis/interface/ExRootTreeReader.h"
18#include "Utilities/ExRootAnalysis/interface/ExRootTreeWriter.h"
19#include "Utilities/ExRootAnalysis/interface/ExRootTreeBranch.h"
20#include "Utilities/ExRootAnalysis/interface/LHEF.h"
21#include "Utilities/stdhep/stdhep_mcfio.h"
22#include "Utilities/stdhep/stdhep_declarations.h"
23
24#include "interface/STDHEPConverter.h"
25
26
27
28using namespace std;
29
30//---------------------------------------------------------------------------
31
32void STDHEPConverter::AnalyseEvent(ExRootTreeBranch *branch, const Long64_t eventNumber)
33{
34 TRootGenEvent *element;
35
36 element = static_cast<TRootGenEvent*>(branch->NewEntry());
37
38 element->Number = eventNumber;
39}
40
41//---------------------------------------------------------------------------
42
43void STDHEPConverter::AnalyseParticles(ExRootTreeBranch *branch)
44{
45 TRootGenParticle *element;
46
47 Double_t signPz;
48 TLorentzVector momentum;
49 Int_t number;
50
51 for(number = 0; number < myhepevt.nhep; ++number)
52 {
53
54 element = static_cast<TRootGenParticle*>(branch->NewEntry());
55
56 element->PID = myhepevt.idhep[number];
57 element->Status = myhepevt.isthep[number];
58 element->M1 = myhepevt.jmohep[number][0] - 1;
59 element->M2 = myhepevt.jmohep[number][1] - 1;
60 element->D1 = myhepevt.jdahep[number][0] - 1;
61 element->D2 = myhepevt.jdahep[number][1] - 1;
62// element->Charge = myhepevt.hepchg(element->PID);
63
64 element->E = myhepevt.phep[number][3];
65 element->Px = myhepevt.phep[number][0];
66 element->Py = myhepevt.phep[number][1];
67 element->Pz = myhepevt.phep[number][2];
68
69 momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
70 element->PT = momentum.Perp();
71 signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
72 element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta();
73 element->Phi = momentum.Phi();
74
75 element->T = myhepevt.vhep[number][3];
76 element->X = myhepevt.vhep[number][0];
77 element->Y = myhepevt.vhep[number][1];
78 element->Z = myhepevt.vhep[number][2];
79 }
80}
81
82
83//------------------------------------------------------------------------------
84
85STDHEPConverter::~STDHEPConverter()
86{
87}
88
89//------------------------------------------------------------------------------
90
91STDHEPConverter::STDHEPConverter(const string& inputFileList, const string& outputFileName)
92{
93
94 int ierr, entryType;
95 int istr = 0;
96 int nevt = 0;
97
98 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
99
100 // information about generated event
101 ExRootTreeBranch *branchGenEvent = treeWriter->NewBranch("Event", TRootGenEvent::Class());
102 // generated particles from HEPEVT
103 ExRootTreeBranch *branchGenParticle = treeWriter->NewBranch("Particle", TRootGenParticle::Class());
104
105 // Open a stream connected to an event file:
106 ifstream infile(inputFileList.c_str());
107 string filename;
108 if(!infile.is_open())
109 {
110 cerr << "** ERROR: Can't open '" << inputFileList << "' for input" << endl;
111 exit(1);
112 }
113 while(1)
114 {
115 infile >> filename;
116 if(!infile.good()) break;
117
118 ierr = StdHepXdrReadInit(const_cast<char*>(filename.c_str()), &nevt, istr);
119
120 if(ierr != 0)
121 {
122 cerr << "** ERROR: Can't open file for input" << endl;
123 break;
124 }
125
126 Long64_t allEntries = nevt;
127 cout << "** Input file contains " << allEntries << " entries" << endl;
128
129 if(allEntries > 0)
130 {
131 // Loop over all objects
132 Long64_t entry = 0;
133 Long64_t recordNumber = 1;
134 for(entry = 0; entry < allEntries; ++entry)
135 {
136 ierr = StdHepXdrRead(&entryType, istr);
137
138 if(ierr != 0)
139 {
140 cerr << "** ERROR: Unexpected end of file after " << entry << " entries" << endl;
141 break;
142 }
143
144 // analyse only entries with standard HEPEVT common block
145 if(entryType == 1)
146 {
147 // add empty events for missing event numbers
148 while(recordNumber < myhepevt.nevhep)
149 {
150 treeWriter->Clear();
151 AnalyseEvent(branchGenEvent, recordNumber);
152 treeWriter->Fill();
153 ++recordNumber;
154 }
155
156 treeWriter->Clear();
157
158 AnalyseEvent(branchGenEvent, myhepevt.nevhep);
159 AnalyseParticles(branchGenParticle);
160
161 treeWriter->Fill();
162
163 ++recordNumber;
164
165 }
166 }
167 }
168 }
169 treeWriter->Write();
170 cout << "** Exiting..." << endl;
171
172 delete treeWriter;
173 StdHepXdrEnd(istr);
174
175}
176
177
178
179
Note: See TracBrowser for help on using the repository browser.