CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SelectReplacementCandidates.cc
Go to the documentation of this file.
3 
5 {
6  using namespace edm;
7  using namespace std;
8  produces< vector<uint32_t> >();
9  produces< vector<uint32_t> >("assocHitsWithHCAL");
10  produces< std::vector<reco::Muon> >();
11 
12  targetParticleMass_=iConfig.getUntrackedParameter<double>("targetParticlesMass",1.77690);
13  targetParticlePdgID_=iConfig.getUntrackedParameter<int>("targetParticlesPdgID",15);
14 
15  muonInputTag_ = iConfig.getParameter<edm::InputTag>("muonInputTag");
16  edm::ParameterSet parameters = iConfig.getParameter<edm::ParameterSet>("TrackAssociatorParameters");
19 
21 }
22 
24 {
25 
26 }
27 
29 {
30  using namespace std;
31  using namespace edm;
32  reco::Muon muon1, muon2;
33 
34  // determine the muons to be used
35  // these can be reconstructed muons or muons from a HepMC::Event
36  if (determineMuonsToUse(iEvent, iSetup, &muon1, &muon2)!=0)
37  return;
38 
39  vector<uint32_t> hits ;
40  getRawIDsAdvanced(iEvent, iSetup, &hits, &muon1, false);
41  getRawIDsAdvanced(iEvent, iSetup, &hits, &muon2, false);
42  std::auto_ptr< vector<uint32_t> > selectedHitsAutoPtr(new vector<uint32_t>(hits) );
43  iEvent.put( selectedHitsAutoPtr );
44 
45  vector<uint32_t> assoc_hits_withHCAL;
46  getRawIDsAdvanced(iEvent, iSetup, &assoc_hits_withHCAL, &muon1, true);
47  getRawIDsAdvanced(iEvent, iSetup, &assoc_hits_withHCAL, &muon2, true);
48  std::cout << "found in total " << assoc_hits_withHCAL.size() << " cells with associated hits with hcal\n";
49  std::auto_ptr< vector<uint32_t> > selectedAssocHitsWithHCALAutoPtr(new vector<uint32_t>(assoc_hits_withHCAL) );
50  iEvent.put( selectedAssocHitsWithHCALAutoPtr, "assocHitsWithHCAL");
51 
52  std::vector<reco::Muon> muons;
53  transformMuMu2TauTau(&muon1, &muon2);
54  muons.push_back(muon1);
55  muons.push_back(muon2);
56  std::auto_ptr< std::vector<reco::Muon> > selectedMuonsPtr(new std::vector<reco::Muon>(muons) );
57  iEvent.put( selectedMuonsPtr );
58 
59 }
60 
62 {
63 
64 }
65 
67 {
68 
69 }
70 
71 
73 {
74  using namespace edm;
75  using namespace reco;
76  using namespace std;
77 
79  if (!iEvent.getByLabel("dimuonsGlobal", zCandidate_handle))
80  {
81  std::cout << "Could not find product: " << "dimuonsGlobal" << "\n";
82  std::vector< edm::Handle< edm::RefToBaseVector<reco::Candidate> > > allHandles;
83  iEvent.getManyByType(allHandles);
84  std::vector< edm::Handle< edm::RefToBaseVector<reco::Candidate> > >::iterator it;
85  for (it = allHandles.begin(); it != allHandles.end(); it++)
86  {
87  std::cout << "available product: " << (*it).provenance()->moduleLabel() << ", " << (*it).provenance()->productInstanceName() << ", " << (*it).provenance()->processName();
88  }
89 
90  std::cout << "Objekt nicht gefunden: dimuonsGloal\n";
91  return -1;
92  }
93 // std::cout << zCandidate_handle->size() << " Kandidaten gefunden!\n";
94 
95  unsigned int nMuons = zCandidate_handle->size();
96  if (nMuons==0)
97  return -1;
98 
99  for (edm::RefToBaseVector<reco::Candidate>::const_iterator z = zCandidate_handle->begin(); z!=zCandidate_handle->end(); ++z)
100  {
101  reco::Particle::LorentzVector muon1p4 = z->get()->daughter(0)->p4();
102  reco::Particle::LorentzVector muon2p4 = z->get()->daughter(1)->p4();
103 
104  edm::Handle<edm::View<reco::Muon> > trackCollection;
105  iEvent.getByLabel(muonInputTag_, trackCollection);
106  const edm::View<reco::Muon>& muons = * trackCollection;
107 
108  bool found1=false, found2=false;
109  for (unsigned int i=0;i<muons.size() && !(found1 && found2);i++)
110  {
111  if (deltaR(muon1p4,muons[i].p4())<0.1)
112  {
113  *muon1 = muons[i];
114  found1=true;
115  }
116  if (deltaR(muon2p4,muons[i].p4())<0.1)
117  {
118  *muon2 = muons[i];
119  found2=true;
120  }
121  }
122 
123  break;
124  }
125  return 0;
126 }
127 
128 
130 {
131  using namespace edm;
132  using namespace reco;
133  using namespace std;
134 
135  Handle<MuonCollection> muonHandle;
136  if (!iEvent.getByLabel(muonInputTag_,muonHandle))
137  ProductNotFound<reco::MuonCollection>(iEvent, muonInputTag_);
138 
139  edm::Handle<edm::View<reco::Muon> > trackCollection;
140  iEvent.getByLabel(muonInputTag_, trackCollection);
141  const edm::View<reco::Muon>& muons = * trackCollection;
142 
143  unsigned int nMuons = muons.size();
144  if (nMuons<2)
145  return -1;
146 
147  *muon1 = muons[0];
148  *muon2 = muons[1];
149  std::cout << muons[0].p4() << "\n";
150  std::cout << muons[1].p4() << "\n";
151 
152  return 0;
153 }
154 
155 void SelectReplacementCandidates::getRawIDsAdvanced(const edm::Event& iEvent, const edm::EventSetup& iConfig, std::vector<uint32_t> * L, reco::Muon * muon, bool includeHCAL)
156 {
157  using namespace edm;
158  using namespace reco;
159  using namespace std;
160 
161  if (muon->bestTrackRef().isNonnull())
162  {
163  TrackDetMatchInfo info = trackAssociator_.associate(iEvent, iConfig, *(muon->bestTrackRef().get()), parameters_);
164  if (includeHCAL)
165  {
166  for(std::vector<const EcalRecHit*>::const_iterator hit = info.ecalRecHits.begin(); hit != info.ecalRecHits.end(); ++hit)
167  L->push_back((*hit)->detid().rawId());
168 
169  for(std::vector<const HBHERecHit*>::const_iterator hit = info.hcalRecHits.begin(); hit != info.hcalRecHits.end(); ++hit)
170  L->push_back((*hit)->detid().rawId());
171  }
172  int recs = (muon->bestTrackRef().get())->recHitsSize();
173  for (int i=0; i<recs; i++)
174  L->push_back(((muon->bestTrackRef().get())->recHit(i).get())->geographicalId().rawId());
175  }
176  else
177  std::cout << "ERROR: Muon has no bestTrackRef!!\n";
178 
179  std::cout << " with " << L->size() << " muon hits found\n";
180 // return L;
181 }
182 
183 
184 template < typename T > void SelectReplacementCandidates::ProductNotFound(const edm::Event& iEvent, edm::InputTag inputTag)
185 {
186  std::cout << "--- TauHitSelector ------------------------------------\n";
187  std::cout << "Could not find product with:\n"<< inputTag << "\n";
188  std::vector< edm::Handle< T > > allHandles;
189  iEvent.getManyByType(allHandles);
190  typename std::vector< edm::Handle< T > >::iterator it;
191  for (it = allHandles.begin(); it != allHandles.end(); it++)
192  {
193  std::cout << "module label: " << (*it).provenance()->moduleLabel() << "\n";
194  std::cout << "productInstanceName: " << (*it).provenance()->productInstanceName() << "\n";
195  std::cout << "processName: " << (*it).provenance()->processName() << "\n\n";
196  }
197  std::cout << "-------------------------------------------------------\n";
198  return;
199 }
200 
201 
204 {
205  using namespace edm;
206  using namespace reco;
207  using namespace std;
208 
209  reco::Particle::LorentzVector muon1_momentum = muon1->p4();
210  reco::Particle::LorentzVector muon2_momentum = muon2->p4();
211  reco::Particle::LorentzVector z_momentum = muon1_momentum + muon2_momentum;
212 
213  ROOT::Math::Boost booster(z_momentum.BoostToCM());
214  ROOT::Math::Boost invbooster(booster.Inverse());
215 
216  reco::Particle::LorentzVector Zb = booster(z_momentum);
217 
218  reco::Particle::LorentzVector muon1b = booster(muon1_momentum);
219  reco::Particle::LorentzVector muon2b = booster(muon2_momentum);
220 
221  double tau_mass2 = targetParticleMass_*targetParticleMass_;
222 
223  double muonxb_mom2 = muon1b.x()*muon1b.x() + muon1b.y()*muon1b.y() + muon1b.z() * muon1b.z();
224  double tauxb_mom2 = 0.25 * Zb.t() * Zb.t() - tau_mass2;
225 
226  float scaling1 = sqrt(tauxb_mom2 / muonxb_mom2);
227  float scaling2 = scaling1;
228 
229  float tauEnergy= Zb.t() / 2.;
230 
231  if (tauEnergy*tauEnergy<tau_mass2)
232  return;
233 
234  reco::Particle::LorentzVector tau1b_mom = reco::Particle::LorentzVector(scaling1*muon1b.x(),scaling1*muon1b.y(),scaling1*muon1b.z(),tauEnergy);
235  reco::Particle::LorentzVector tau2b_mom = reco::Particle::LorentzVector(scaling2*muon2b.x(),scaling2*muon2b.y(),scaling2*muon2b.z(),tauEnergy);
236 
237  // some checks
238  // the following test guarantees a deviation
239  // of less than 0.1% for phi and theta for the
240  // original muons and the placed taus
241  // (in the centre-of-mass system of the z boson)
242  assert((muon1b.phi()-tau1b_mom.phi())/muon1b.phi()<0.001);
243  assert((muon2b.phi()-tau2b_mom.phi())/muon2b.phi()<0.001);
244  assert((muon1b.theta()-tau1b_mom.theta())/muon1b.theta()<0.001);
245  assert((muon2b.theta()-tau2b_mom.theta())/muon2b.theta()<0.001);
246 
247  reco::Particle::LorentzVector tau1_mom = (invbooster(tau1b_mom));
248  reco::Particle::LorentzVector tau2_mom = (invbooster(tau2b_mom));
249 
250  // some additional checks
251  // the following tests guarantee a deviation of less
252  // than 0.1% for the following values of the original
253  // muons and the placed taus
254  // invariant mass
255  // transverse momentum
256  assert(((muon1_momentum+muon1_momentum).mass()-(tau1_mom+tau2_mom).mass())/(muon1_momentum+muon1_momentum).mass()<0.001);
257  assert(((muon1_momentum+muon2_momentum).pt()-(tau1_mom+tau2_mom).pt())/(muon1_momentum+muon1_momentum).pt()<0.001);
258 
259  muon1->setP4(tau1_mom);
260  muon2->setP4(tau2_mom);
261 
262  muon1->setPdgId(targetParticlePdgID_*muon1->pdgId()/abs(muon1->pdgId()));
263  muon2->setPdgId(targetParticlePdgID_*muon2->pdgId()/abs(muon2->pdgId()));
264 
265  muon1->setStatus(1);
266  muon2->setStatus(1);
267 
268  return;
269 }
270 
272 
274 
void getManyByType(std::vector< Handle< PROD > > &results) const
Definition: Event.h:436
value_type const * get() const
Definition: RefToBase.h:213
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
virtual int pdgId() const
PDG identifier.
dictionary parameters
Definition: Parameters.py:2
static const TGPicture * info(bool iBackgroundIsBlack)
TrackDetectorAssociator trackAssociator_
virtual void setStatus(int status)
set status word
void ProductNotFound(const edm::Event &iEvent, edm::InputTag inputTag)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
std::vector< const EcalRecHit * > ecalRecHits
hits in the cone
void useDefaultPropagator()
use the default propagator
void loadParameters(const edm::ParameterSet &, edm::ConsumesCollector &)
virtual void setP4(const LorentzVector &p4)
set 4-momentum
size_type size() const
bool isNonnull() const
Checks for non-null.
Definition: RefToBase.h:280
float float float z
SelectReplacementCandidates(const edm::ParameterSet &iSetup)
int iEvent
Definition: GenABIO.cc:230
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:113
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
T sqrt(T t)
Definition: SSEVec.h:48
double p4[4]
Definition: TauolaWrapper.h:92
std::vector< const HBHERecHit * > hcalRecHits
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:402
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
int determineMuonsToUse(const edm::Event &iEvent, const edm::EventSetup &iConfig, reco::Muon *muon1, reco::Muon *muon2)
void getRawIDsAdvanced(const edm::Event &iEvent, const edm::EventSetup &iConfig, std::vector< uint32_t > *L, reco::Muon *muon, bool includeHCAL)
void transformMuMu2TauTau(reco::Muon *muon1, reco::Muon *muon2)
transform muon into tau
tuple muons
Definition: patZpeak.py:38
virtual void setPdgId(int pdgId)
TrackDetMatchInfo associate(const edm::Event &, const edm::EventSetup &, const FreeTrajectoryState &, const AssociatorParameters &)
tuple cout
Definition: gather_cfg.py:121
TrackAssociatorParameters parameters_
virtual void produce(edm::Event &iEvent, const edm::EventSetup &iConfig)
virtual const LorentzVector & p4() const
four-momentum Lorentz vector
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:27
int determineMuonsToUse_old(const edm::Event &iEvent, const edm::EventSetup &iConfig, reco::Muon *muon1, reco::Muon *muon2)
virtual TrackBaseRef bestTrackRef() const
best track RefToBase
Definition: Muon.h:62