Fork me on GitHub

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

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

checks that the input file is present

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