ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/CmsHi/JetAnalysis/src/HeavyIonJetAnalyzer.cc
Revision: 1.3
Committed: Wed May 6 19:09:22 2009 UTC (16 years ago) by yilmaz
Content type: text/plain
Branch: MAIN
CVS Tags: cmshi_22X
Branch point for: BRANCH22X
Changes since 1.2: +13 -9 lines
Log Message:
Some fixes

File Contents

# User Rev Content
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 yilmaz 1.3 ncoll = hi->Ncoll();
310     nhard = hi->Ncoll_hard();
311 yilmaz 1.1
312     if(printLists_){
313     out_b<<b<<endl;
314     out_n<<npart<<endl;
315     }
316    
317     }
318    
319     // edm::Handle<reco::GenJetCollection> genjets;
320     edm::Handle<reco::JetView> genjets;
321     // edm::Handle<vector<reco::Jet> > genjets;
322     iEvent.getByLabel(jetSrc_[0],genjets);
323     for(int ijet = 0; ijet < genjets->size(); ++ijet){
324    
325     const reco::Jet* jet = &((*genjets)[ijet]);
326     // const reco::GenJet* jet = dynamic_cast<const reco::GenJet*>(&((*genjets)[ijet]));
327    
328     cout<<"Jet Quantities : "<<jet->pt()<<" "<<jet->eta()<<" "<<jet->phi()<<" "<<jet->jetArea()<<endl;
329    
330     hev_.et[hev_.njet] = jet->pt();
331     hev_.eta[hev_.njet] = jet->eta();
332     hev_.phi[hev_.njet] = jet->phi();
333     // hev_.area[hev_.njet] = jet->jetArea();
334    
335     ++(hev_.njet);
336 yilmaz 1.2
337 yilmaz 1.1 }
338    
339 yilmaz 1.3
340     edm::Handle<edm::SimVertexContainer> simVertices;
341 yilmaz 1.1 iEvent.getByType<edm::SimVertexContainer>(simVertices);
342 yilmaz 1.3
343 yilmaz 1.1 if (! simVertices.isValid() ) throw cms::Exception("FatalError") << "No vertices found\n";
344     int inum = 0;
345 yilmaz 1.3
346 yilmaz 1.1 edm::SimVertexContainer::const_iterator it=simVertices->begin();
347     SimVertex vertex = (*it);
348     cout<<" Vertex position "<< inum <<" " << vertex.position().rho()<<" "<<vertex.position().z()<<endl;
349     vx = vertex.position().x();
350     vy = vertex.position().y();
351     vz = vertex.position().z();
352     vr = vertex.position().rho();
353 yilmaz 1.3
354 yilmaz 1.1
355 yilmaz 1.3 /*
356     for(int i = 0; i<3; ++i){
357     hev_.ptav[i] = hev_.ptav[i]/hev_.n[i];
358     }
359     */
360    
361 yilmaz 1.1 hev_.b = b;
362     hev_.npart = npart;
363     hev_.ncoll = ncoll;
364     hev_.nhard = nhard;
365     hev_.vx = vx;
366     hev_.vy = vy;
367     hev_.vz = vz;
368     hev_.vr = vr;
369    
370     hydjetTree_->Fill();
371    
372     }
373    
374    
375     // ------------ method called once each job just before starting event loop ------------
376     void
377     HeavyIonJetAnalyzer::beginJob(const edm::EventSetup& iSetup)
378     {
379     iSetup.getData(pdt);
380    
381     if(printLists_){
382     out_b.open(fBFileName.c_str());
383     if(out_b.good() == false)
384     throw cms::Exception("BadFile") << "Can\'t open file " << fBFileName;
385     out_n.open(fNFileName.c_str());
386     if(out_n.good() == false)
387     throw cms::Exception("BadFile") << "Can\'t open file " << fNFileName;
388     out_m.open(fMFileName.c_str());
389     if(out_m.good() == false)
390     throw cms::Exception("BadFile") << "Can\'t open file " << fMFileName;
391     }
392    
393     if(doAnalysis_){
394     hydjetTree_ = f->make<TTree>("hi","Tree of Hydjet Events");
395    
396     hydjetTree_->Branch("event",&hev_.event,"event/I");
397     hydjetTree_->Branch("b",&hev_.b,"b/F");
398     hydjetTree_->Branch("npart",&hev_.npart,"npart/F");
399     hydjetTree_->Branch("ncoll",&hev_.ncoll,"ncoll/F");
400     hydjetTree_->Branch("nhard",&hev_.nhard,"nhard/F");
401    
402     if(doParticles_){
403     hydjetTree_->Branch("n",hev_.n,"n[3]/I");
404     hydjetTree_->Branch("ptav",hev_.ptav,"ptav[3]/F");
405     hydjetTree_->Branch("np",&hev_.np,"np/I");
406     hydjetTree_->Branch("par_pt",hev_.par_pt,"par_pt[np]/F");
407     hydjetTree_->Branch("par_eta",hev_.par_eta,"par_eta[np]/F");
408     hydjetTree_->Branch("par_phi",hev_.par_phi,"par_phi[np]/F");
409     hydjetTree_->Branch("pdg",hev_.pdg,"pdg[np]/I");
410     hydjetTree_->Branch("chg",hev_.chg,"chg[np]/I");
411     }
412    
413     hydjetTree_->Branch("vx",&hev_.vx,"vx/F");
414     hydjetTree_->Branch("vy",&hev_.vy,"vy/F");
415     hydjetTree_->Branch("vz",&hev_.vz,"vz/F");
416     hydjetTree_->Branch("vr",&hev_.vr,"vr/F");
417    
418     hydjetTree_->Branch("njet",&hev_.njet,"njet/I");
419     hydjetTree_->Branch("et",hev_.et,"et[njet]/F");
420     hydjetTree_->Branch("eta",hev_.eta,"eta[njet]/F");
421     hydjetTree_->Branch("phi",hev_.phi,"phi[njet]/F");
422     hydjetTree_->Branch("area",hev_.area,"area[njet]/F");
423    
424     }
425    
426     }
427    
428     // ------------ method called once each job just after ending the event loop ------------
429     void
430     HeavyIonJetAnalyzer::endJob() {
431     }
432    
433     //define this as a plug-in
434     DEFINE_FWK_MODULE(HeavyIonJetAnalyzer);