1 |
yilmaz |
1.1 |
// -*- C++ -*-
|
2 |
|
|
//
|
3 |
|
|
// Package: HeavyIonJetAnalyzer
|
4 |
|
|
// Class: HeavyIonJetAnalyzer
|
5 |
|
|
//
|
6 |
|
|
/**\class HeavyIonJetAnalyzer HeavyIonJetAnalyzer.cc yetkin/HeavyIonJetAnalyzer/src/HeavyIonJetAnalyzer.cc
|
7 |
|
|
|
8 |
|
|
Description: <one line class summary>
|
9 |
|
|
|
10 |
|
|
Implementation:
|
11 |
|
|
<Notes on implementation>
|
12 |
|
|
*/
|
13 |
|
|
//
|
14 |
|
|
// Original Author: Yetkin Yilmaz
|
15 |
|
|
// Created: Tue Dec 18 09:44:41 EST 2007
|
16 |
yilmaz |
1.2 |
// $Id: HeavyIonJetAnalyzer.cc,v 1.2 2008/12/19 19:04:12 yilmaz Exp $
|
17 |
yilmaz |
1.1 |
//
|
18 |
|
|
//
|
19 |
|
|
|
20 |
|
|
|
21 |
|
|
// system include files
|
22 |
|
|
#include <memory>
|
23 |
|
|
#include <vector>
|
24 |
|
|
#include <iostream>
|
25 |
|
|
#include <string>
|
26 |
|
|
#include <fstream>
|
27 |
|
|
|
28 |
|
|
// user include files
|
29 |
|
|
#include "FWCore/Framework/interface/Frameworkfwd.h"
|
30 |
|
|
#include "FWCore/Framework/interface/EDAnalyzer.h"
|
31 |
|
|
|
32 |
|
|
#include "FWCore/Framework/interface/Event.h"
|
33 |
|
|
#include "FWCore/Framework/interface/MakerMacros.h"
|
34 |
|
|
#include "FWCore/Framework/interface/EventSetup.h"
|
35 |
|
|
|
36 |
|
|
#include "FWCore/ParameterSet/interface/ParameterSet.h"
|
37 |
|
|
#include "FWCore/ParameterSet/interface/InputTag.h"
|
38 |
|
|
#include "FWCore/Framework/interface/ESHandle.h"
|
39 |
|
|
|
40 |
|
|
#include "FWCore/ServiceRegistry/interface/Service.h"
|
41 |
|
|
|
42 |
|
|
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
|
43 |
|
|
|
44 |
|
|
#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"
|
45 |
|
|
#include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
|
46 |
|
|
#include "SimDataFormats/Vertex/interface/SimVertex.h"
|
47 |
|
|
#include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
|
48 |
|
|
|
49 |
|
|
#include "DataFormats/JetReco/interface/JetCollection.h"
|
50 |
|
|
#include "DataFormats/JetReco/interface/Jet.h"
|
51 |
|
|
#include "DataFormats/JetReco/interface/GenJetCollection.h"
|
52 |
|
|
#include "DataFormats/JetReco/interface/GenJet.h"
|
53 |
|
|
|
54 |
|
|
|
55 |
|
|
#include "HepMC/GenEvent.h"
|
56 |
|
|
#include "HepMC/HeavyIon.h"
|
57 |
|
|
|
58 |
|
|
#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
|
59 |
|
|
|
60 |
|
|
// root include file
|
61 |
|
|
#include "TFile.h"
|
62 |
|
|
#include "TNtuple.h"
|
63 |
|
|
|
64 |
|
|
using namespace std;
|
65 |
|
|
|
66 |
|
|
|
67 |
|
|
#define PI 3.14159265358979
|
68 |
|
|
|
69 |
|
|
#define MAXPARTICLES 5000000
|
70 |
|
|
#define MAXJETS 500000
|
71 |
|
|
|
72 |
|
|
#define MAXHITS 50000
|
73 |
|
|
#define MAXVTX 1000
|
74 |
|
|
#define ETABINS 3 // Fix also in branch string
|
75 |
|
|
|
76 |
|
|
//
|
77 |
|
|
// class decleration
|
78 |
|
|
//
|
79 |
|
|
|
80 |
|
|
struct HydjetEvent{
|
81 |
|
|
|
82 |
|
|
int event;
|
83 |
|
|
float b;
|
84 |
|
|
float npart;
|
85 |
|
|
float ncoll;
|
86 |
|
|
float nhard;
|
87 |
|
|
|
88 |
|
|
int n[ETABINS];
|
89 |
|
|
float ptav[ETABINS];
|
90 |
|
|
|
91 |
|
|
int np;
|
92 |
|
|
float par_pt[MAXPARTICLES];
|
93 |
|
|
float par_eta[MAXPARTICLES];
|
94 |
|
|
float par_phi[MAXPARTICLES];
|
95 |
|
|
int pdg[MAXPARTICLES];
|
96 |
|
|
int chg[MAXPARTICLES];
|
97 |
|
|
|
98 |
|
|
int algos;
|
99 |
|
|
int njet;
|
100 |
|
|
|
101 |
|
|
float et[MAXJETS];
|
102 |
|
|
float eta[MAXJETS];
|
103 |
|
|
float phi[MAXJETS];
|
104 |
|
|
float area[MAXJETS];
|
105 |
|
|
|
106 |
|
|
float vx;
|
107 |
|
|
float vy;
|
108 |
|
|
float vz;
|
109 |
|
|
float vr;
|
110 |
|
|
|
111 |
|
|
};
|
112 |
|
|
|
113 |
|
|
class HeavyIonJetAnalyzer : public edm::EDAnalyzer {
|
114 |
|
|
public:
|
115 |
|
|
explicit HeavyIonJetAnalyzer(const edm::ParameterSet&);
|
116 |
|
|
~HeavyIonJetAnalyzer();
|
117 |
|
|
|
118 |
|
|
|
119 |
|
|
private:
|
120 |
|
|
virtual void beginJob(const edm::EventSetup&) ;
|
121 |
|
|
virtual void analyze(const edm::Event&, const edm::EventSetup&);
|
122 |
|
|
virtual void endJob() ;
|
123 |
|
|
|
124 |
|
|
// ----------member data ---------------------------
|
125 |
|
|
|
126 |
|
|
std::ofstream out_b;
|
127 |
|
|
std::string fBFileName;
|
128 |
|
|
|
129 |
|
|
std::ofstream out_n;
|
130 |
|
|
std::string fNFileName;
|
131 |
|
|
|
132 |
|
|
std::ofstream out_m;
|
133 |
|
|
std::string fMFileName;
|
134 |
|
|
|
135 |
|
|
|
136 |
|
|
TTree* hydjetTree_;
|
137 |
|
|
HydjetEvent hev_;
|
138 |
|
|
|
139 |
|
|
TNtuple *ntpart;
|
140 |
|
|
|
141 |
|
|
std::string output; // Output filename
|
142 |
|
|
|
143 |
|
|
bool doAnalysis_;
|
144 |
|
|
bool doParticles_;
|
145 |
|
|
bool printLists_;
|
146 |
|
|
bool doCF_;
|
147 |
|
|
double etaMax_;
|
148 |
|
|
double ptMin_;
|
149 |
|
|
|
150 |
|
|
vector<string> jetSrc_;
|
151 |
|
|
|
152 |
|
|
edm::ESHandle < ParticleDataTable > pdt;
|
153 |
|
|
edm::Service<TFileService> f;
|
154 |
|
|
|
155 |
|
|
|
156 |
|
|
};
|
157 |
|
|
|
158 |
|
|
//
|
159 |
|
|
// constants, enums and typedefs
|
160 |
|
|
//
|
161 |
|
|
|
162 |
|
|
//
|
163 |
|
|
// static data member definitions
|
164 |
|
|
//
|
165 |
|
|
|
166 |
|
|
//
|
167 |
|
|
// constructors and destructor
|
168 |
|
|
//
|
169 |
|
|
HeavyIonJetAnalyzer::HeavyIonJetAnalyzer(const edm::ParameterSet& iConfig)
|
170 |
|
|
{
|
171 |
|
|
//now do what ever initialization is needed
|
172 |
|
|
fBFileName = iConfig.getUntrackedParameter<std::string>("output_b", "b_values.txt");
|
173 |
|
|
fNFileName = iConfig.getUntrackedParameter<std::string>("output_n", "n_values.txt");
|
174 |
|
|
fMFileName = iConfig.getUntrackedParameter<std::string>("output_m", "m_values.txt");
|
175 |
|
|
doAnalysis_ = iConfig.getUntrackedParameter<bool>("doAnalysis", true);
|
176 |
|
|
doParticles_ = iConfig.getUntrackedParameter<bool>("doParticles", true);
|
177 |
|
|
printLists_ = iConfig.getUntrackedParameter<bool>("printLists", false);
|
178 |
|
|
doCF_ = iConfig.getUntrackedParameter<bool>("doMixed", false);
|
179 |
|
|
jetSrc_ = iConfig.getParameter<vector<string> >("jetSrc");
|
180 |
|
|
|
181 |
|
|
etaMax_ = iConfig.getUntrackedParameter<double>("etaMax", 2);
|
182 |
|
|
ptMin_ = iConfig.getUntrackedParameter<double>("ptMin", 0);
|
183 |
|
|
|
184 |
|
|
// Output
|
185 |
|
|
|
186 |
|
|
}
|
187 |
|
|
|
188 |
|
|
|
189 |
|
|
HeavyIonJetAnalyzer::~HeavyIonJetAnalyzer()
|
190 |
|
|
{
|
191 |
|
|
// do anything here that needs to be done at desctruction time
|
192 |
|
|
// (e.g. close files, deallocate resources etc.)
|
193 |
|
|
|
194 |
|
|
}
|
195 |
|
|
|
196 |
|
|
|
197 |
|
|
//
|
198 |
|
|
// member functions
|
199 |
|
|
//
|
200 |
|
|
|
201 |
|
|
// ------------ method called to for each event ------------
|
202 |
|
|
void
|
203 |
|
|
HeavyIonJetAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
|
204 |
|
|
{
|
205 |
|
|
using namespace edm;
|
206 |
|
|
using namespace HepMC;
|
207 |
|
|
|
208 |
|
|
hev_.event = iEvent.id().event();
|
209 |
|
|
for(int ieta = 0; ieta < ETABINS; ++ieta){
|
210 |
|
|
hev_.n[ieta] = 0;
|
211 |
|
|
hev_.ptav[ieta] = 0;
|
212 |
|
|
}
|
213 |
|
|
|
214 |
|
|
hev_.njet = 0;
|
215 |
|
|
hev_.np = 0;
|
216 |
|
|
|
217 |
|
|
double b = -1;
|
218 |
|
|
int npart = -1;
|
219 |
|
|
int ncoll = -1;
|
220 |
|
|
int nhard = -1;
|
221 |
|
|
double vx = -99;
|
222 |
|
|
double vy = -99;
|
223 |
|
|
double vz = -99;
|
224 |
|
|
double vr = -99;
|
225 |
|
|
const GenEvent* evt;
|
226 |
yilmaz |
1.2 |
|
227 |
yilmaz |
1.1 |
if(doCF_){
|
228 |
|
|
|
229 |
|
|
Handle<CrossingFrame<HepMCProduct> > cf;
|
230 |
|
|
iEvent.getByLabel(InputTag("mix","source"),cf);
|
231 |
|
|
|
232 |
|
|
/*
|
233 |
|
|
|
234 |
|
|
MixCollection<HepMCProduct> mix(cf.product());
|
235 |
|
|
|
236 |
|
|
int mixsize = mix.size();
|
237 |
|
|
|
238 |
|
|
cout<<"Mix Collection Size: "<<mixsize<<endl;
|
239 |
|
|
evt = mix.getObject(mixsize-1).GetEvent();
|
240 |
|
|
|
241 |
|
|
MixCollection<HepMCProduct>::iterator begin = mix.begin();
|
242 |
|
|
MixCollection<HepMCProduct>::iterator end = mix.end();
|
243 |
|
|
|
244 |
|
|
for(MixCollection<HepMCProduct>::iterator mixit = begin; mixit != end; ++mixit){
|
245 |
|
|
|
246 |
|
|
const GenEvent* subevt = (*mixit).GetEvent();
|
247 |
|
|
int all = subevt->particles_size();
|
248 |
|
|
HepMC::GenEvent::particle_const_iterator begin = subevt->particles_begin();
|
249 |
|
|
HepMC::GenEvent::particle_const_iterator end = subevt->particles_end();
|
250 |
|
|
for(HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it){
|
251 |
|
|
if((*it)->status() == 1){
|
252 |
|
|
float pdg_id = (*it)->pdg_id();
|
253 |
|
|
float eta = (*it)->momentum().eta();
|
254 |
|
|
float pt = (*it)->momentum().perp();
|
255 |
|
|
const ParticleData * part = pdt->particle(pdg_id );
|
256 |
|
|
float charge = part->charge();
|
257 |
|
|
}
|
258 |
|
|
}
|
259 |
|
|
}
|
260 |
|
|
|
261 |
|
|
}
|
262 |
|
|
|
263 |
|
|
*/
|
264 |
|
|
|
265 |
|
|
}else{
|
266 |
|
|
|
267 |
|
|
Handle<HepMCProduct> mc;
|
268 |
|
|
iEvent.getByLabel("source",mc);
|
269 |
|
|
evt = mc->GetEvent();
|
270 |
|
|
|
271 |
|
|
if(doParticles_){
|
272 |
yilmaz |
1.2 |
int all = evt->particles_size();
|
273 |
|
|
HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
|
274 |
|
|
HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
|
275 |
|
|
for(HepMC::GenEvent::particle_const_iterator it = begin; it != end; ++it){
|
276 |
|
|
if((*it)->status() == 1){
|
277 |
|
|
int pdg_id = (*it)->pdg_id();
|
278 |
|
|
float eta = (*it)->momentum().eta();
|
279 |
|
|
float phi = (*it)->momentum().phi();
|
280 |
|
|
float pt = (*it)->momentum().perp();
|
281 |
|
|
const ParticleData * part = pdt->particle(pdg_id );
|
282 |
|
|
int charge = part->charge();
|
283 |
|
|
|
284 |
|
|
hev_.par_pt[hev_.np] = pt;
|
285 |
|
|
hev_.par_eta[hev_.np] = eta;
|
286 |
|
|
hev_.par_phi[hev_.np] = phi;
|
287 |
|
|
hev_.pdg[hev_.np] = pdg_id;
|
288 |
|
|
hev_.chg[hev_.np] = charge;
|
289 |
|
|
|
290 |
|
|
eta = fabs(eta);
|
291 |
|
|
int etabin = 0;
|
292 |
|
|
if(eta > 0.5) etabin = 1;
|
293 |
|
|
if(eta > 1.) etabin = 2;
|
294 |
|
|
if(eta < 2.){
|
295 |
|
|
hev_.ptav[etabin] += pt;
|
296 |
|
|
++(hev_.n[etabin]);
|
297 |
|
|
}
|
298 |
|
|
++(hev_.np);
|
299 |
|
|
}
|
300 |
|
|
}
|
301 |
|
|
}
|
302 |
yilmaz |
1.1 |
|
303 |
|
|
}
|
304 |
|
|
|
305 |
|
|
const HeavyIon* hi = evt->heavy_ion();
|
306 |
|
|
if(hi){
|
307 |
|
|
b = hi->impact_parameter();
|
308 |
|
|
npart = hi->Npart_proj()+hi->Npart_targ();
|
309 |
|
|
|
310 |
|
|
if(printLists_){
|
311 |
|
|
out_b<<b<<endl;
|
312 |
|
|
out_n<<npart<<endl;
|
313 |
|
|
}
|
314 |
|
|
|
315 |
|
|
}
|
316 |
|
|
|
317 |
|
|
// edm::Handle<reco::GenJetCollection> genjets;
|
318 |
|
|
edm::Handle<reco::JetView> genjets;
|
319 |
|
|
// edm::Handle<vector<reco::Jet> > genjets;
|
320 |
|
|
iEvent.getByLabel(jetSrc_[0],genjets);
|
321 |
|
|
for(int ijet = 0; ijet < genjets->size(); ++ijet){
|
322 |
|
|
|
323 |
|
|
const reco::Jet* jet = &((*genjets)[ijet]);
|
324 |
|
|
// const reco::GenJet* jet = dynamic_cast<const reco::GenJet*>(&((*genjets)[ijet]));
|
325 |
|
|
|
326 |
|
|
cout<<"Jet Quantities : "<<jet->pt()<<" "<<jet->eta()<<" "<<jet->phi()<<" "<<jet->jetArea()<<endl;
|
327 |
|
|
|
328 |
|
|
hev_.et[hev_.njet] = jet->pt();
|
329 |
|
|
hev_.eta[hev_.njet] = jet->eta();
|
330 |
|
|
hev_.phi[hev_.njet] = jet->phi();
|
331 |
|
|
// hev_.area[hev_.njet] = jet->jetArea();
|
332 |
|
|
|
333 |
|
|
++(hev_.njet);
|
334 |
yilmaz |
1.2 |
|
335 |
yilmaz |
1.1 |
}
|
336 |
|
|
|
337 |
|
|
/*
|
338 |
|
|
edm::Handle<edm::SimVertexContainer> simVertices;
|
339 |
|
|
iEvent.getByType<edm::SimVertexContainer>(simVertices);
|
340 |
|
|
|
341 |
|
|
if (! simVertices.isValid() ) throw cms::Exception("FatalError") << "No vertices found\n";
|
342 |
|
|
int inum = 0;
|
343 |
|
|
|
344 |
|
|
edm::SimVertexContainer::const_iterator it=simVertices->begin();
|
345 |
|
|
SimVertex vertex = (*it);
|
346 |
|
|
cout<<" Vertex position "<< inum <<" " << vertex.position().rho()<<" "<<vertex.position().z()<<endl;
|
347 |
|
|
vx = vertex.position().x();
|
348 |
|
|
vy = vertex.position().y();
|
349 |
|
|
vz = vertex.position().z();
|
350 |
|
|
vr = vertex.position().rho();
|
351 |
|
|
|
352 |
|
|
for(int i = 0; i<3; ++i){
|
353 |
|
|
hev_.ptav[i] = hev_.ptav[i]/hev_.n[i];
|
354 |
|
|
}
|
355 |
|
|
*/
|
356 |
|
|
|
357 |
|
|
hev_.b = b;
|
358 |
|
|
hev_.npart = npart;
|
359 |
|
|
hev_.ncoll = ncoll;
|
360 |
|
|
hev_.nhard = nhard;
|
361 |
|
|
hev_.vx = vx;
|
362 |
|
|
hev_.vy = vy;
|
363 |
|
|
hev_.vz = vz;
|
364 |
|
|
hev_.vr = vr;
|
365 |
|
|
|
366 |
|
|
hydjetTree_->Fill();
|
367 |
|
|
|
368 |
|
|
}
|
369 |
|
|
|
370 |
|
|
|
371 |
|
|
// ------------ method called once each job just before starting event loop ------------
|
372 |
|
|
void
|
373 |
|
|
HeavyIonJetAnalyzer::beginJob(const edm::EventSetup& iSetup)
|
374 |
|
|
{
|
375 |
|
|
iSetup.getData(pdt);
|
376 |
|
|
|
377 |
|
|
if(printLists_){
|
378 |
|
|
out_b.open(fBFileName.c_str());
|
379 |
|
|
if(out_b.good() == false)
|
380 |
|
|
throw cms::Exception("BadFile") << "Can\'t open file " << fBFileName;
|
381 |
|
|
out_n.open(fNFileName.c_str());
|
382 |
|
|
if(out_n.good() == false)
|
383 |
|
|
throw cms::Exception("BadFile") << "Can\'t open file " << fNFileName;
|
384 |
|
|
out_m.open(fMFileName.c_str());
|
385 |
|
|
if(out_m.good() == false)
|
386 |
|
|
throw cms::Exception("BadFile") << "Can\'t open file " << fMFileName;
|
387 |
|
|
}
|
388 |
|
|
|
389 |
|
|
if(doAnalysis_){
|
390 |
|
|
hydjetTree_ = f->make<TTree>("hi","Tree of Hydjet Events");
|
391 |
|
|
|
392 |
|
|
hydjetTree_->Branch("event",&hev_.event,"event/I");
|
393 |
|
|
hydjetTree_->Branch("b",&hev_.b,"b/F");
|
394 |
|
|
hydjetTree_->Branch("npart",&hev_.npart,"npart/F");
|
395 |
|
|
hydjetTree_->Branch("ncoll",&hev_.ncoll,"ncoll/F");
|
396 |
|
|
hydjetTree_->Branch("nhard",&hev_.nhard,"nhard/F");
|
397 |
|
|
|
398 |
|
|
if(doParticles_){
|
399 |
|
|
hydjetTree_->Branch("n",hev_.n,"n[3]/I");
|
400 |
|
|
hydjetTree_->Branch("ptav",hev_.ptav,"ptav[3]/F");
|
401 |
|
|
hydjetTree_->Branch("np",&hev_.np,"np/I");
|
402 |
|
|
hydjetTree_->Branch("par_pt",hev_.par_pt,"par_pt[np]/F");
|
403 |
|
|
hydjetTree_->Branch("par_eta",hev_.par_eta,"par_eta[np]/F");
|
404 |
|
|
hydjetTree_->Branch("par_phi",hev_.par_phi,"par_phi[np]/F");
|
405 |
|
|
hydjetTree_->Branch("pdg",hev_.pdg,"pdg[np]/I");
|
406 |
|
|
hydjetTree_->Branch("chg",hev_.chg,"chg[np]/I");
|
407 |
|
|
}
|
408 |
|
|
|
409 |
|
|
hydjetTree_->Branch("vx",&hev_.vx,"vx/F");
|
410 |
|
|
hydjetTree_->Branch("vy",&hev_.vy,"vy/F");
|
411 |
|
|
hydjetTree_->Branch("vz",&hev_.vz,"vz/F");
|
412 |
|
|
hydjetTree_->Branch("vr",&hev_.vr,"vr/F");
|
413 |
|
|
|
414 |
|
|
hydjetTree_->Branch("njet",&hev_.njet,"njet/I");
|
415 |
|
|
hydjetTree_->Branch("et",hev_.et,"et[njet]/F");
|
416 |
|
|
hydjetTree_->Branch("eta",hev_.eta,"eta[njet]/F");
|
417 |
|
|
hydjetTree_->Branch("phi",hev_.phi,"phi[njet]/F");
|
418 |
|
|
hydjetTree_->Branch("area",hev_.area,"area[njet]/F");
|
419 |
|
|
|
420 |
|
|
}
|
421 |
|
|
|
422 |
|
|
}
|
423 |
|
|
|
424 |
|
|
// ------------ method called once each job just after ending the event loop ------------
|
425 |
|
|
void
|
426 |
|
|
HeavyIonJetAnalyzer::endJob() {
|
427 |
|
|
}
|
428 |
|
|
|
429 |
|
|
//define this as a plug-in
|
430 |
|
|
DEFINE_FWK_MODULE(HeavyIonJetAnalyzer);
|