CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EvtPlaneProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: EvtPlaneProducer
4 // Class: EvtPlaneProducer
5 //
13 //
14 // Original Author: Sergey Petrushanko
15 // Created: Fri Jul 11 10:05:00 2008
16 // $Id: EvtPlaneProducer.cc,v 1.21 2012/02/15 11:04:09 eulisse Exp $
17 //
18 //
19 #define TRACKCOLLECTION 1
20 //#define RECOCHARGEDCANDIDATECOLLECTION 1
21 
22 
23 // system include files
24 #include <memory>
25 #include <iostream>
26 #include <time.h>
27 #include "TMath.h"
28 
29 // user include files
30 #include "FWCore/Framework/interface/Frameworkfwd.h"
31 #include "FWCore/Framework/interface/EDProducer.h"
32 #include "FWCore/Framework/interface/Event.h"
33 #include "FWCore/Framework/interface/EventSetup.h"
34 
35 #include "FWCore/Framework/interface/MakerMacros.h"
36 #include "FWCore/Framework/interface/ESHandle.h"
37 #include "FWCore/ParameterSet/interface/ParameterSet.h"
38 
39 #include "DataFormats/Candidate/interface/Candidate.h"
40 #include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
41 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
42 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
43 
44 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
45 #include "DataFormats/Common/interface/EDProduct.h"
46 #include "DataFormats/Common/interface/Ref.h"
47 
48 #include "DataFormats/Common/interface/Handle.h"
49 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
50 
51 #include "FWCore/ServiceRegistry/interface/Service.h"
52 #include "CommonTools/UtilAlgos/interface/TFileService.h"
53 #include "DataFormats/TrackReco/interface/TrackFwd.h"
54 #include "DataFormats/TrackReco/interface/Track.h"
55 #include "DataFormats/VertexReco/interface/Vertex.h"
56 #include "DataFormats/VertexReco/interface/VertexFwd.h"
57 
58 #include <cstdlib>
59 #include "RecoHI/HiEvtPlaneAlgos/interface/HiEvtPlaneList.h"
60 
61 using namespace std;
62 using namespace hi;
63 
64 //
65 // class decleration
66 //
67 
68 class EvtPlaneProducer : public edm::EDProducer {
69 public:
70  explicit EvtPlaneProducer(const edm::ParameterSet&);
71  ~EvtPlaneProducer();
72 
73 private:
74  //edm::Service<TFileService> fs;
75  class GenPlane {
76  public:
77  GenPlane(string name,double etaminval1,double etamaxval1,double etaminval2,double etamaxval2,int orderval){
78  epname=name;
79  etamin1=etaminval1;
80  etamax1=etamaxval1;
81  etamin2=etaminval2;
82  etamax2=etamaxval2;
83  sumsin=0;
84  sumcos=0;
85  order = (double) orderval;
86  }
87  ~GenPlane(){;}
88  void addParticle(double w, double s, double c, double eta) {
89  if(w < 0.001) return;
90  if((eta>=etamin1 && eta<etamax1) ||
91  (etamin2!= etamax2 && eta>=etamin2 && eta<etamax2 )) {
92  sumsin+=w*s;
93  sumcos+=w*c;
94  //For tracking, w=1 and mult is the track multiplicity.
95  //For calorimetors, w is an energy that needs to be subsequently
96  //converted if an multiplicity if needed
97  mult+=w;
98  }
99  }
100 
101  double getAngle(double &ang, double &sv, double &cv){
102  ang = -10;
103  sv = 0;
104  cv = 0;
105  sv = sumsin;
106  cv = sumcos;
107  double q = sv*sv+cv*cv;
108  if(q>0) ang = atan2(sv,cv)/order;
109  return ang;
110  }
111  void reset() {
112  sumsin=0;
113  sumcos=0;
114  mult = 0;
115  }
116  private:
117  string epname;
118  double etamin1;
119  double etamax1;
120 
121  double etamin2;
122  double etamax2;
123  double sumsin;
124  double sumcos;
125  int mult;
126  double order;
127  };
128 
129 
130  GenPlane *rp[NumEPNames];
131 
132  virtual void beginJob() ;
133  virtual void produce(edm::Event&, const edm::EventSetup&);
134  virtual void endJob() ;
135 
136  // ----------member data ---------------------------
137  edm::InputTag vtxCollection_;
138  edm::InputTag caloCollection_;
139  edm::InputTag trackCollection_;
140  bool useECAL_;
141  bool useHCAL_;
142  bool useTrack_;
143  bool useTrackPtWeight_;
144  double minet_;
145  double maxet_;
146  double minpt_;
147  double maxpt_;
148  double minvtx_;
149  double maxvtx_;
150  double dzerr_;
151  double chi2_;
152 
153  bool storeNames_;
154 };
155 
156 //
157 // constants, enums and typedefs
158 //
159 
160 
161 //
162 // static data member definitions
163 //
164 
165 //
166 // constructors and destructor
167 //
168 EvtPlaneProducer::EvtPlaneProducer(const edm::ParameterSet& iConfig)
169 {
170 
171 
172  //register your products
173  vtxCollection_ = iConfig.getParameter<edm::InputTag>("vtxCollection_");
174  caloCollection_ = iConfig.getParameter<edm::InputTag>("caloCollection_");
175  trackCollection_ = iConfig.getParameter<edm::InputTag>("trackCollection_");
176  useECAL_ = iConfig.getUntrackedParameter<bool>("useECAL_",true);
177  useHCAL_ = iConfig.getUntrackedParameter<bool>("useHCAL_",true);
178  useTrack_ = iConfig.getUntrackedParameter<bool>("useTrack",true);
179  useTrackPtWeight_ = iConfig.getUntrackedParameter<bool>("useTrackPtWeight_",true);
180  minet_ = iConfig.getUntrackedParameter<double>("minet_",0.2);
181  maxet_ = iConfig.getUntrackedParameter<double>("maxet_",500.);
182  minpt_ = iConfig.getUntrackedParameter<double>("minpt_",0.3);
183  maxpt_ = iConfig.getUntrackedParameter<double>("maxpt_",2.5);
184  minvtx_ = iConfig.getUntrackedParameter<double>("minvtx_",-50.);
185  maxvtx_ = iConfig.getUntrackedParameter<double>("maxvtx_",50.);
186  dzerr_ = iConfig.getUntrackedParameter<double>("dzerr_",10.);
187  chi2_ = iConfig.getUntrackedParameter<double>("chi2_",40.);
188 
189  storeNames_ = 1;
190  produces<reco::EvtPlaneCollection>("recoLevel");
191  for(int i = 0; i<NumEPNames; i++ ) {
192  rp[i] = new GenPlane(EPNames[i].data(),EPEtaMin1[i],EPEtaMax1[i],EPEtaMin2[i],EPEtaMax2[i],EPOrder[i]);
193  }
194 }
195 
196 
197 EvtPlaneProducer::~EvtPlaneProducer()
198 {
199 
200  // do anything here that needs to be done at desctruction time
201  // (e.g. close files, deallocate resources etc.)
202 
203 }
204 
205 
206 //
207 // member functions
208 //
209 
210 // ------------ method called to produce the data ------------
211 void
212 EvtPlaneProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
213 {
214  using namespace edm;
215  using namespace reco;
216 
217  int vs_sell;
218  float vzr_sell;
219  //
220  //Get Vertex
221  //
222  edm::Handle<reco::VertexCollection> vertexCollection3;
223  iEvent.getByLabel(vtxCollection_,vertexCollection3);
224  const reco::VertexCollection * vertices3 = vertexCollection3.product();
225  vs_sell = vertices3->size();
226  if(vs_sell>0) {
227  vzr_sell = vertices3->begin()->z();
228  } else
229  vzr_sell = -999.9;
230  //
231  for(int i = 0; i<NumEPNames; i++) rp[i]->reset();
232  if(vzr_sell>minvtx_ && vzr_sell<maxvtx_) {
233  //calorimetry part
234 
235  double tower_eta, tower_phi;
236  double tower_energyet, tower_energyet_e, tower_energyet_h;
237  double s1, s2, s13,s23,s14,s24,s15,s25,s16,s26;
238  Handle<CaloTowerCollection> calotower;
239  iEvent.getByLabel(caloCollection_,calotower);
240 
241  if(calotower.isValid()){
242 
243  for (CaloTowerCollection::const_iterator j = calotower->begin();j !=calotower->end(); j++) {
244  tower_eta = j->eta();
245  tower_phi = j->phi();
246  tower_energyet_e = j->emEt();
247  tower_energyet_h = j->hadEt();
248  tower_energyet = tower_energyet_e + tower_energyet_h;
249 
250  s1 = sin(2.*tower_phi);
251  s2 = cos(2.*tower_phi);
252  s13 = sin(3.*tower_phi);
253  s23 = cos(3.*tower_phi);
254  s14 = sin(4.*tower_phi);
255  s24 = cos(4.*tower_phi);
256  s15 = sin(5.*tower_phi);
257  s25 = cos(5.*tower_phi);
258  s16 = sin(6.*tower_phi);
259  s26 = cos(6.*tower_phi);
260 
261  if(tower_energyet<minet_) continue;
262  if(tower_energyet>maxet_) continue;;
263 
264  rp[etEcal ]->addParticle (tower_energyet_e, s1, s2, tower_eta);
265  rp[etEcalP ]->addParticle (tower_energyet_e, s1, s2, tower_eta);
266  rp[etEcalM ]->addParticle (tower_energyet_e, s1, s2, tower_eta);
267  rp[etHcal ]->addParticle (tower_energyet_h, s1, s2, tower_eta);
268  rp[etHcalP ]->addParticle (tower_energyet_h, s1, s2, tower_eta);
269  rp[etHcalM ]->addParticle (tower_energyet_h, s1, s2, tower_eta);
270  rp[etHF ]->addParticle (tower_energyet, s1, s2, tower_eta);
271  rp[etHFp ]->addParticle (tower_energyet, s1, s2, tower_eta);
272  rp[etHFm ]->addParticle (tower_energyet, s1, s2, tower_eta);
273  rp[etCaloHFP ]->addParticle (tower_energyet, s1, s2, tower_eta);
274  rp[etCaloHFM ]->addParticle (tower_energyet, s1, s2, tower_eta);
275 
276  rp[etHF3 ]->addParticle (tower_energyet, s13, s23, tower_eta);
277  rp[etHFp3 ]->addParticle (tower_energyet, s13, s23, tower_eta);
278  rp[etHFm3 ]->addParticle (tower_energyet, s13, s23, tower_eta);
279 
280  rp[etHF4 ]->addParticle (tower_energyet, s14, s24, tower_eta);
281  rp[etHFp4 ]->addParticle (tower_energyet, s14, s24, tower_eta);
282  rp[etHFm4 ]->addParticle (tower_energyet, s14, s24, tower_eta);
283 
284  rp[etHF5 ]->addParticle (tower_energyet, s15, s25, tower_eta);
285  rp[etHFp5 ]->addParticle (tower_energyet, s15, s25, tower_eta);
286  rp[etHFm5 ]->addParticle (tower_energyet, s15, s25, tower_eta);
287 
288  rp[etHF6 ]->addParticle (tower_energyet, s16, s26, tower_eta);
289  rp[etHFp6 ]->addParticle (tower_energyet, s16, s26, tower_eta);
290  rp[etHFm6 ]->addParticle (tower_energyet, s16, s26, tower_eta);
291  }
292  }
293  //Tracking part
294 
295  double track_eta;
296  double track_phi;
297  double track_pt;
298 #ifdef TRACKCOLLECTION
299  Handle<reco::TrackCollection> tracks;
300  iEvent.getByLabel(trackCollection_, tracks);
301 
302  if(tracks.isValid()){
303  for(reco::TrackCollection::const_iterator j = tracks->begin(); j != tracks->end(); j++){
304 
305  //Find possible collections with command line: edmDumpEventContent *.root
306 #endif
307 #ifdef RECOCHARGEDCANDIDATECOLLECTION
308  edm::Handle<reco::RecoChargedCandidateCollection> trackCollection;
309  iEvent.getByLabel("allMergedPtSplit12Tracks",trackCollection);
310  // iEvent.getByLabel("hiGoodMergedTracks",trackCollection);
311  if(trackCollection.isValid()){
312 
313  const reco::RecoChargedCandidateCollection * tracks = trackCollection.product();
314  for(reco::RecoChargedCandidateCollection::const_iterator j = tracks->begin(); j != tracks->end(); j++){
315 #endif
316  edm::Handle<reco::VertexCollection> vertex;
317  iEvent.getByLabel(vtxCollection_, vertex);
318 
319 // find the vertex point and error
320 
321  math::XYZPoint vtxPoint(0.0,0.0,0.0);
322  double vzErr =0.0, vxErr=0.0, vyErr=0.0;
323  if(vertex->size()>0) {
324  vtxPoint=vertex->begin()->position();
325  vzErr=vertex->begin()->zError();
326  vxErr=vertex->begin()->xError();
327  vyErr=vertex->begin()->yError();
328  }
329  bool accepted = true;
330  bool isPixel = false;
331  // determine if the track is a pixel track
332  if ( j->numberOfValidHits() < 7 ) isPixel = true;
333 
334  // determine the vertex significance
335  double d0=0.0, dz=0.0, d0sigma=0.0, dzsigma=0.0;
336  d0 = -1.*j->dxy(vtxPoint);
337  dz = j->dz(vtxPoint);
338  d0sigma = sqrt(j->d0Error()*j->d0Error()+vxErr*vyErr);
339  dzsigma = sqrt(j->dzError()*j->dzError()+vzErr*vzErr);
340 
341  // cuts for pixel tracks
342  if( isPixel )
343  {
344  // dz significance cut
345  if ( fabs(dz/dzsigma) > dzerr_ ) accepted = false;
346 
347  // chi2/ndof cut
348  if ( j->normalizedChi2() > chi2_ ) accepted = false;
349  }
350 
351  // cuts for full tracks
352  if ( ! isPixel)
353  {
354  // dz and d0 significance cuts
355  if ( fabs(dz/dzsigma) > 3 ) accepted = false;
356  if ( fabs(d0/d0sigma) > 3 ) accepted = false;
357 
358  // pt resolution cut
359  if ( j->ptError()/j->pt() > 0.05 ) accepted = false;
360 
361  // number of valid hits cut
362  if ( j->numberOfValidHits() < 12 ) accepted = false;
363  }
364  if( accepted ) {
365  track_eta = j->eta();
366  track_phi = j->phi();
367  track_pt = j->pt();
368  double s =sin(2*track_phi);
369  double c =cos(2*track_phi);
370  double s3 =sin(3*track_phi);
371  double c3 =cos(3*track_phi);
372  double s4 =sin(4*track_phi);
373  double c4 =cos(4*track_phi);
374  double s5 =sin(5*track_phi);
375  double c5 =cos(5*track_phi);
376  double s6 =sin(6*track_phi);
377  double c6 =cos(6*track_phi);
378  double w = 1;
379  if(useTrackPtWeight_) {
380  w = track_pt;
381  if(w>2.5) w=2.0; //v2 starts decreasing above ~2.5 GeV/c
382  }
383  if(track_pt<minpt_) continue;
384  if(track_pt>maxpt_) continue;
385  rp[EvtPlaneFromTracksMidEta]->addParticle(w,s,c,track_eta);
386  rp[EvtPTracksPosEtaGap]->addParticle(w,s,c,track_eta);
387  rp[EvtPTracksNegEtaGap]->addParticle(w,s,c,track_eta);
388  rp[EPTracksMid3]->addParticle(w,s3,c3,track_eta);
389  rp[EPTracksPos3]->addParticle(w,s3,c3,track_eta);
390  rp[EPTracksNeg3]->addParticle(w,s3,c3,track_eta);
391  rp[EPTracksMid4]->addParticle(w,s4,c4,track_eta);
392  rp[EPTracksPos4]->addParticle(w,s4,c4,track_eta);
393  rp[EPTracksNeg4]->addParticle(w,s4,c4,track_eta);
394  rp[EPTracksMid5]->addParticle(w,s5,c5,track_eta);
395  rp[EPTracksPos5]->addParticle(w,s5,c5,track_eta);
396  rp[EPTracksNeg5]->addParticle(w,s5,c5,track_eta);
397  rp[EPTracksMid6]->addParticle(w,s6,c6,track_eta);
398  rp[EPTracksPos6]->addParticle(w,s6,c6,track_eta);
399  rp[EPTracksNeg6]->addParticle(w,s6,c6,track_eta);
400 
401  }
402 
403  }
404  }
405  }
406  std::auto_ptr<EvtPlaneCollection> evtplaneOutput(new EvtPlaneCollection);
407 
408  EvtPlane *ep[NumEPNames];
409 
410  double ang=-10;
411  double sv = 0;
412  double cv = 0;
413 
414  for(int i = 0; i<NumEPNames; i++) {
415  rp[i]->getAngle(ang,sv,cv);
416  if(storeNames_) ep[i] = new EvtPlane(ang,sv,cv,EPNames[i]);
417  else ep[i] = new EvtPlane(ang,sv,cv,"");
418  }
419  if(useTrack_) {
420  for(int i = 0; i<9; i++) {
421  evtplaneOutput->push_back(*ep[i]);
422  }
423  }
424  for(int i = 9; i<NumEPNames; i++) {
425  if(useECAL_ && !useHCAL_) {
426  if(EPNames[i].rfind("Ecal")!=string::npos) {
427  evtplaneOutput->push_back(*ep[i]);
428  }
429  } else if (useHCAL_ && !useECAL_) {
430  if(EPNames[i].rfind("Hcal")!=string::npos) {
431  evtplaneOutput->push_back(*ep[i]);
432  }
433  }else if (useECAL_ && useHCAL_) {
434  evtplaneOutput->push_back(*ep[i]);
435  }
436  }
437 
438  iEvent.put(evtplaneOutput, "recoLevel");
439  // storeNames_ = 0;
440 }
441 
442  // ------------ method called once each job just before starting event loop ------------
443 void
444 EvtPlaneProducer::beginJob()
445 {
446 }
447 
448 // ------------ method called once each job just after ending the event loop ------------
449 void
450 EvtPlaneProducer::endJob() {
451 }
452 
453 //define this as a plug-in
454 DEFINE_FWK_MODULE(EvtPlaneProducer);
EvtPlaneProducer::GenPlane::getAngle
double getAngle(double &ang, double &sv, double &cv)
Definition: EvtPlaneProducer.cc:101
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
hi::EPOrder
const int EPOrder[]
Definition: HiEvtPlaneList.h:120
i
int i
Definition: DBlmapReader.cc:9
hi::EPEtaMin2
const double EPEtaMin2[]
Definition: HiEvtPlaneList.h:167
EvtPlaneProducer::GenPlane::GenPlane
GenPlane(string name, double etaminval1, double etamaxval1, double etaminval2, double etamaxval2, int orderval)
Definition: EvtPlaneProducer.cc:77
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
edm::SortedCollection< CaloTower >::const_iterator
std::vector< CaloTower >::const_iterator const_iterator
Definition: SortedCollection.h:80
reco::VertexCollection
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
eta
T eta() const
Definition: Basic3DVectorLD.h:179
hi::EPNames
const std::string EPNames[]
Definition: HiEvtPlaneList.h:103
bk::beginJob
void beginJob()
Definition: Breakpoints.cc:15
indexGen.s2
tuple s2
Definition: indexGen.py:106
hi::EPEtaMin1
const double EPEtaMin1[]
Definition: HiEvtPlaneList.h:136
EvtPlaneProducer::endJob
virtual void endJob()
Definition: EvtPlaneProducer.cc:450
iEvent
int iEvent
Definition: GenABIO.cc:243
EvtPlaneProducer::beginJob
virtual void beginJob()
Definition: EvtPlaneProducer.cc:444
edm::Event::put
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:94
L1Trigger_dataformats.reco
dictionary reco
Definition: L1Trigger_dataformats.py:9
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
reco::EvtPlaneCollection
std::vector< EvtPlane > EvtPlaneCollection
Definition: EvtPlane.h:34
EvtPlaneProducer::produce
virtual void produce(edm::Event &, const edm::EventSetup &)
Definition: EvtPlaneProducer.cc:212
j
int j
Definition: DBlmapReader.cc:9
EvtPlaneProducer::EvtPlaneProducer
EvtPlaneProducer(const edm::ParameterSet &)
Definition: EvtPlaneProducer.cc:168
edm::HandleBase::isValid
bool isValid() const
Definition: HandleBase.h:76
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
cuy.cv
dictionary cv
Definition: cuy.py:362
reco::RecoChargedCandidateCollection
std::vector< RecoChargedCandidate > RecoChargedCandidateCollection
collectin of RecoChargedCandidate objects
Definition: RecoChargedCandidateFwd.h:9
EvtPlaneProducer::GenPlane::addParticle
void addParticle(double w, double s, double c, double eta)
Definition: EvtPlaneProducer.cc:88
testEve_cfg.tracks
tuple tracks
Definition: testEve_cfg.py:39
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
edm::Handle::product
T const * product() const
Definition: Handle.h:74
data
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
EvtPlaneProducer::caloCollection_
edm::InputTag caloCollection_
Definition: EvtPlaneProducer.cc:138
EvtPlaneProducer::vtxCollection_
edm::InputTag vtxCollection_
Definition: EvtPlaneProducer.cc:137
hi::EPEtaMax2
const double EPEtaMax2[]
Definition: HiEvtPlaneList.h:182
hi::NumEPNames
static const int NumEPNames
Definition: HiEvtPlaneList.h:198
w
T w() const
Definition: extBasic3DVector.h:245
EvtPlaneProducer::trackCollection_
edm::InputTag trackCollection_
Definition: EvtPlaneProducer.cc:139
reset
void reset(double vett[256])
Definition: TPedValues.cc:11
hi::EPEtaMax1
const double EPEtaMax1[]
Definition: HiEvtPlaneList.h:152
Generated for CMSSW Reference Manual by doxygen 1.8.5