source: trunk/test/ExRootSTDHEPConverter.cpp

Last change on this file was 2, checked in by Pavel Demin, 16 years ago

first commit

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