Fork me on GitHub

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

Last change on this file since 199 was 56, checked in by severine ovyn, 16 years ago

add M

File size: 4.8 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 element->M = myhepevt.phep[number][4];
69
70 momentum.SetPxPyPzE(element->Px, element->Py, element->Pz, element->E);
71 element->PT = momentum.Perp();
72 signPz = (element->Pz >= 0.0) ? 1.0 : -1.0;
73 element->Eta = element->PT == 0.0 ? signPz*999.9 : momentum.Eta();
74 element->Phi = momentum.Phi();
75
76 element->T = myhepevt.vhep[number][3];
77 element->X = myhepevt.vhep[number][0];
78 element->Y = myhepevt.vhep[number][1];
79 element->Z = myhepevt.vhep[number][2];
80 }
81}
82
83
84//------------------------------------------------------------------------------
85
86STDHEPConverter::~STDHEPConverter()
87{
88}
89
90//------------------------------------------------------------------------------
91
92STDHEPConverter::STDHEPConverter(const string& inputFileList, const string& outputFileName)
93{
94
95 int ierr, entryType;
96 int istr = 0;
97 int nevt = 0;
98
99 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
100
101 // information about generated event
102 ExRootTreeBranch *branchGenEvent = treeWriter->NewBranch("Event", TRootGenEvent::Class());
103 // generated particles from HEPEVT
104 ExRootTreeBranch *branchGenParticle = treeWriter->NewBranch("Particle", TRootGenParticle::Class());
105
106 // Open a stream connected to an event file:
107 ifstream infile(inputFileList.c_str());
108 string filename;
109 if(!infile.is_open())
110 {
111 cerr << "** ERROR: Can't open '" << inputFileList << "' for input" << endl;
112 exit(1);
113 }
114 while(1)
115 {
116 infile >> filename;
117 if(!infile.good()) break;
118 ifstream checking_the_file(filename.c_str());
119 if(!checking_the_file.good()) { cout << filename << ": file not found\n"; continue;}
120 else checking_the_file.close();
121
122 ierr = StdHepXdrReadInit(const_cast<char*>(filename.c_str()), &nevt, istr);
123
124 if(ierr != 0)
125 {
126 cerr << "** ERROR: Can't open file for input" << endl;
127 break;
128 }
129
130 Long64_t allEntries = nevt;
131 cout << "** Input file contains " << allEntries << " entries" << endl;
132
133 if(allEntries > 0)
134 {
135 // Loop over all objects
136 Long64_t entry = 0;
137 Long64_t recordNumber = 1;
138 for(entry = 0; entry < allEntries; ++entry)
139 {
140 ierr = StdHepXdrRead(&entryType, istr);
141
142 if(ierr != 0)
143 {
144 cerr << "** ERROR: Unexpected end of file after " << entry << " entries" << endl;
145 break;
146 }
147
148 // analyse only entries with standard HEPEVT common block
149 if(entryType == 1)
150 {
151 // add empty events for missing event numbers
152 while(recordNumber < myhepevt.nevhep)
153 {
154 treeWriter->Clear();
155 AnalyseEvent(branchGenEvent, recordNumber);
156 treeWriter->Fill();
157 ++recordNumber;
158 }
159
160 treeWriter->Clear();
161
162 AnalyseEvent(branchGenEvent, myhepevt.nevhep);
163 AnalyseParticles(branchGenParticle);
164
165 treeWriter->Fill();
166
167 ++recordNumber;
168
169 }
170 }
171 }
172 }
173 treeWriter->Write();
174 cout << "** Exiting..." << endl;
175
176 delete treeWriter;
177 StdHepXdrEnd(istr);
178
179}
180
181
182
183
Note: See TracBrowser for help on using the repository browser.