Fork me on GitHub

source: svn/trunk/src/HEPEVTConverter.cc@ 256

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

OK for error message

File size: 6.2 KB
Line 
1#include <iostream>
2#include <utility>
3#include <deque>
4#include <sstream>
5#include <fstream>
6
7#include "TROOT.h"
8#include "TTree.h"
9#include "TBranch.h"
10#include "TLeaf.h"
11#include "TString.h"
12#include "TLorentzVector.h"
13#include "TChain.h"
14
15#include "BlockClasses.h"
16#include "ExRootTreeWriter.h"
17#include "ExRootTreeBranch.h"
18
19#include "HEPEVTConverter.h"
20
21using namespace std;
22
23//------------------------------------------------------------------------------
24
25//------------------------------------------------------------------------------
26
27HEPTreeReader::HEPTreeReader(TTree *tree, HEPEvent *event) : fChain(tree), fCurrentTree(-1), fEvent(event)
28{
29 if(!fChain) return;
30
31 TBranch *branch;
32 TLeaf *leaf;
33 TString name;
34 Int_t i;
35
36 name = "Nhep";
37 branch = fChain->GetBranch(name);
38 branch->SetAddress(&(event->Nhep));
39 fBranches.push_back(make_pair(name, branch));
40
41 TString intNames[3] = {"Idhep", "Jsmhep", "Jsdhep"};
42 Int_t **intData[3] = {&event->Idhep, &event->Jsmhep, &event->Jsdhep};
43
44 for(i = 0; i < 3; ++i)
45 {
46 name = intNames[i];
47 branch = fChain->GetBranch(name);
48 leaf = branch->GetLeaf(name);
49 *intData[i] = new Int_t[leaf->GetNdata()];
50 branch->SetAddress(*intData[i]);
51 fBranches.push_back(make_pair(name, branch));
52 }
53
54 TString floatNames[2] = {"Phep", "Vhep"};
55 Float_t **floatData[2] = {&event->Phep, &event->Vhep};
56
57 for(i = 0; i < 2; ++i)
58 {
59 name = floatNames[i];
60 branch = fChain->GetBranch(name);
61 leaf = branch->GetLeaf(name);
62 *floatData[i] = new Float_t[leaf->GetNdata()];
63 branch->SetAddress(*floatData[i]);
64 fBranches.push_back(make_pair(name, branch));
65 }
66}
67
68//------------------------------------------------------------------------------
69
70HEPTreeReader::~HEPTreeReader()
71{
72 if(!fChain) return;
73
74 delete[] fEvent->Idhep;
75 delete[] fEvent->Jsmhep;
76 delete[] fEvent->Jsdhep;
77 delete[] fEvent->Phep;
78 delete[] fEvent->Vhep;
79}
80
81//------------------------------------------------------------------------------
82
83Bool_t HEPTreeReader::ReadEntry(Long64_t element)
84{
85 if(!fChain) return kFALSE;
86
87 Int_t treeEntry = fChain->LoadTree(element);
88 if(treeEntry < 0) return kFALSE;
89
90 if(fChain->IsA() == TChain::Class())
91 {
92 TChain *chain = static_cast<TChain*>(fChain);
93 if(chain->GetTreeNumber() != fCurrentTree)
94 {
95 fCurrentTree = chain->GetTreeNumber();
96 Notify();
97 }
98 }
99
100 deque< pair<TString, TBranch*> >::iterator it_deque;
101 TBranch *branch;
102
103 for(it_deque = fBranches.begin(); it_deque != fBranches.end(); ++it_deque)
104 {
105 branch = it_deque->second;
106 if(branch)
107 {
108 branch->GetEntry(treeEntry);
109 }
110 }
111
112 return kTRUE;
113}
114
115//------------------------------------------------------------------------------
116
117void HEPTreeReader::Notify()
118{
119 // Called when loading a new file.
120 // Get branch pointers.
121 if(!fChain) return;
122
123 deque< pair<TString, TBranch*> >::iterator it_deque;
124
125 for(it_deque = fBranches.begin(); it_deque != fBranches.end(); ++it_deque)
126 {
127 it_deque->second = fChain->GetBranch(it_deque->first);
128 if(!it_deque->second)
129 {
130 cerr << left << setw(40) <<"** WARNING: cannot get branch: "<<""
131 << left << setw(20) << it_deque->first <<""
132 << right << setw(9) <<" **"<<""<<endl;
133 }
134 }
135}
136
137
138
139HEPEVTConverter::HEPEVTConverter(const string& inputFileList, const string& outputFileName)
140{
141 string buffer;
142
143 gROOT->SetBatch();
144
145 ExRootTreeWriter *treeWriter = new ExRootTreeWriter(outputFileName, "GEN");
146 ExRootTreeBranch *branchGen = treeWriter->NewBranch("Particle", TRootGenParticle::Class());
147
148 ifstream infile(inputFileList.c_str());
149 if(!infile.is_open())
150 {
151 cerr << left << setw(30) <<"** ERROR: Can't open "<<""
152 << left << setw(20) << inputFileList <<""
153 << right << setw(19) <<"for input **"<<""<<endl;
154 exit(1);
155 }
156 while(1)
157 {
158
159 TChain *chain = new TChain("h101");
160 infile >> buffer;
161 ifstream checking_the_file(buffer.c_str());
162 if(!checking_the_file.good())
163 {
164 cerr << left << setw(30) <<"** ERROR: Can't find file "<<""
165 << left << setw(20) << buffer <<""
166 << right << setw(19) <<"for input **"<<""<<endl;
167 continue;
168 }
169 else checking_the_file.close();
170
171 if(!infile.good()) break;
172 chain->Add(buffer.c_str());
173
174 HEPEvent event;
175 HEPTreeReader *treeReader = new HEPTreeReader(chain, &event);
176
177 Long64_t entry, allEntries = treeReader->GetEntries();
178 Int_t address, particle;
179 Float_t signEta;
180
181 TRootGenParticle *element;
182
183 // Loop over all events
184 for(entry = 0; entry < allEntries; ++entry)
185 {
186 // Load selected branches with data from specified event
187 treeReader->ReadEntry(entry);
188 treeWriter->Clear();
189
190 for(particle = 0; particle < event.Nhep; ++particle)
191 {
192 element = (TRootGenParticle*) branchGen->NewEntry();
193
194 element->PID = event.Idhep[particle];
195 element->Status = event.Jsmhep[particle]/16000000
196 + event.Jsdhep[particle]/16000000*100;
197 element->M1 = (event.Jsmhep[particle]%16000000)/4000 -1;
198 element->M2 = event.Jsmhep[particle]%4000 -1;
199 element->D1 = (event.Jsdhep[particle]%16000000)/4000 -1;
200 element->D2 = event.Jsdhep[particle]%4000 -1;
201 address = particle*5;
202 element->E = event.Phep[address + 3];
203 element->Px = event.Phep[address + 0];
204 element->Py = event.Phep[address + 1];
205 element->Pz = event.Phep[address + 2];
206 element->M = event.Phep[address + 4];
207
208 TLorentzVector vector(element->Px, element->Py, element->Pz, element->E);
209 element->PT = vector.Perp();
210 signEta = (element->Pz >= 0.0) ? 1.0 : -1.0;
211 element->Eta = vector.CosTheta()*vector.CosTheta() == 1.0 ? signEta*999.9 : vector.Eta();
212 element->Phi = vector.Phi();
213
214 address = particle*4;
215 element->T = event.Vhep[address + 3];
216 element->X = event.Vhep[address + 0];
217 element->Y = event.Vhep[address + 1];
218 element->Z = event.Vhep[address + 2];
219 }
220 treeWriter->Fill();
221 }
222 delete chain;
223 delete treeReader;
224 }
225 treeWriter->Write();
226
227 delete treeWriter;
228
229}
230
231
232
233
Note: See TracBrowser for help on using the repository browser.