Fork me on GitHub

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

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

vers un Nevnt implemente correctement. Pas fini

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