CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ParticleReplacerParticleGun.cc
Go to the documentation of this file.
3 
5 
6 #include "HepMC/PythiaWrapper.h"
7 #include "HepMC/IO_HEPEVT.h"
8 
10  ParticleReplacerBase(iConfig),
11  tauola_(gen::TauolaInterface::getInstance()),
12  pythia_(iConfig),
13  particleOrigin_(iConfig.getParameter<std::string>("particleOrigin")),
14  forceTauPolarization_(iConfig.getParameter<std::string>("forceTauPolarization")),
15  forceTauDecay_(iConfig.getParameter<std::string>("forceTauDecay")),
16  generatorMode_(iConfig.getParameter<std::string>("generatorMode")),
17  gunParticle_(iConfig.getParameter<int>("gunParticle")),
18  forceTauPlusHelicity_(iConfig.getParameter<int>("forceTauPlusHelicity")),
19  forceTauMinusHelicity_(iConfig.getParameter<int>("forceTauMinusHelicity")),
20  printout_(verbose) {
21  tauola_->setPSet(iConfig.getParameter<edm::ParameterSet>("ExternalDecays").getParameter<edm::ParameterSet>("Tauola"));
22  srand(time(NULL)); // Should we use RandomNumberGenerator service?
23 
24  if(forceTauPlusHelicity_ != 0)
25  edm::LogInfo("MuonReplacement") << "[ParticleReplacer::ParticleReplacer] "
26  << "Forcing tau+ to helicity " << forceTauPlusHelicity_ << std::endl;
27  if(forceTauMinusHelicity_ != 0)
28  edm::LogInfo("MuonReplacement") << "[ParticleReplacer::ParticleReplacer] "
29  << "Forcing tau- to helicity " << forceTauMinusHelicity_ << std::endl;
31  edm::LogInfo("MuonReplacement") << "[ParticleReplacer::ParticleReplacer] "
32  << " Forcing tau helicity as decayed from a " << forceTauPolarization_ << std::endl;
33  if(forceTauDecay_ != "" && forceTauDecay_ != "none")
34  edm::LogInfo("MuonReplacement") << "[ParticleReplacer::ParticleReplacer] "
35  << "Forcing tau decaying into " << forceTauDecay_ << std::endl;
36 
37  std::memset(pol1_, 0, 4*sizeof(float));
38  std::memset(pol2_, 0, 4*sizeof(float));
39 
40  if(generatorMode_ != "Tauola")
41  throw cms::Exception("Configuration") << "Generator mode other than Tauola is not supported" << std::endl;
42 
43  throw cms::Exception("UnimplementedFeature") << "ParticleReplacerParticleGun is not usable yet." << std::endl;
44 }
45 
47 
50 
52 
53  if(abs(gunParticle_) == 15) {
54  /* FIXME
55  call_tauola(-1,1);
56  */
57  }
58 }
59 
61  if(abs(gunParticle_) == 15) {
62  /* FIXME
63  call_tauola(1,1);
64  */
65  }
66 }
67 
68 std::auto_ptr<HepMC::GenEvent> ParticleReplacerParticleGun::produce(const reco::MuonCollection& muons, const reco::Vertex *pvtx, const HepMC::GenEvent *genEvt) {
69  if(genEvt != 0)
70  throw cms::Exception("UnimplementedFeature") << "ParticleReplacerParticleGun does NOT support merging at HepMC level" << std::endl;
71 
72  std::auto_ptr<HepMC::GenEvent> evt(0);
73  std::vector<HepMC::FourVector> muons_corrected;
74  muons_corrected.reserve(muons.size());
75  correctTauMass(muons, muons_corrected);
76 
78 
79  for(unsigned int i=0; i<muons_corrected.size(); ++i) {
80  HepMC::FourVector& muon = muons_corrected[i];
81  call_py1ent(i+1, gunParticle_*muons[i].charge(), muon.e(), muon.theta(), muon.phi());
82  }
83 
84  // Let's not do a call_pyexec here because it's unnecessary
85 
86  if(printout_) {
87  std::cout << " ///////////////////// After py1ent, before pyhepc /////////////////////" << std::endl;
88  call_pylist(3);
89  }
90 
91  // Vertex shift
92  call_pyhepc(1); // pythia -> hepevt
93 
94  if(printout_) {
95  std::cout << " ///////////////////// After pyhepc, before vertex shift /////////////////////" << std::endl;
96  HepMC::HEPEVT_Wrapper::print_hepevt();
97  }
98  // Search for HepMC/HEPEVT_Wrapper.h for the wrapper interface
99  int nparticles = HepMC::HEPEVT_Wrapper::number_entries();
100  HepMC::ThreeVector shift(0,0,0);
101  if(particleOrigin_ == "primaryVertex") {
102  if(!pvtx)
103  throw cms::Exception("LogicError") << "Particle origin set to primaryVertex, but pvtx is null!" << std::endl;
104 
105  shift.setX(pvtx->x()*10); // cm -> mm
106  shift.setY(pvtx->y()*10); // cm -> mm
107  shift.setZ(pvtx->z()*10); // cm -> mm
108  }
109  for(int i=1; i <= nparticles; ++i) {
110  if(abs(HepMC::HEPEVT_Wrapper::id(i)) != abs(gunParticle_)) {
111  throw cms::Exception("LogicError") << "Particle in HEPEVT is " << HepMC::HEPEVT_Wrapper::id(i)
112  << " is not the same as gunParticle " << gunParticle_
113  << " for index " << i << std::endl;
114  }
115 
116  if(particleOrigin_ == "muonReferencePoint") {
117  const reco::Muon& muon = muons[i-1];
118  shift.setX(muon.vx()*10); // cm -> mm
119  shift.setY(muon.vy()*10); // cm -> mm
120  shift.setZ(muon.vz()*10); // cm -> mm
121  }
122 
123  HepMC::HEPEVT_Wrapper::set_position(i,
124  HepMC::HEPEVT_Wrapper::x(i) + shift.x(),
125  HepMC::HEPEVT_Wrapper::y(i) + shift.y(),
126  HepMC::HEPEVT_Wrapper::z(i) + shift.z(),
128  }
129 
130  if(printout_) {
131  std::cout << " ///////////////////// After vertex shift, before pyhepc/tauola /////////////////////" << std::endl;
132  HepMC::HEPEVT_Wrapper::print_hepevt();
133  }
134 
135  if(abs(gunParticle_) == 15){
136  // Code example from TauolaInterface::processEvent()
137  /* FIXME
138  int dummy = -1;
139  int numGenParticles_beforeTAUOLA = call_ihepdim(dummy);
140  */
141 
143 
144  for(unsigned int i=0; i<muons_corrected.size(); ++i) {
145  tauola_forParticleGun(i+1, gunParticle_*muons[i].charge(), muons_corrected[i]);
146  }
147 
148  if(printout_) {
149  std::cout << " ///////////////////// After tauola, before pyhepc /////////////////////" << std::endl;
150  HepMC::HEPEVT_Wrapper::print_hepevt();
151  }
152  call_pyhepc(2); // hepevt->pythia
153  if(printout_) {
154  std::cout << " ///////////////////// After pyhepc, before vertex fix /////////////////////" << std::endl;
155  call_pylist(3);
156  }
157 
158  // Fix tau decay vertex position
159  /* FIXME
160  int numGenParticles_afterTAUOLA = call_ihepdim(dummy);
161  tauola_.setDecayVertex(numGenParticles_beforeTAUOLA, numGenParticles_afterTAUOLA);
162  */
163 
164  if(printout_) {
165  /* FIXME
166  std::cout << " numGenParticles_beforeTAUOLA " << numGenParticles_beforeTAUOLA << std::endl
167  << " numGenParticles_afterTAUOLA " << numGenParticles_afterTAUOLA << std::endl;
168  */
169  std::cout << " ///////////////////// After vertex fix, before pyexec /////////////////////" << std::endl;
170  call_pylist(3);
171  }
172  }
173  else {
174  call_pyhepc(2); // hepevt->pythia
175  if(printout_) {
176  std::cout << " ///////////////////// After pyhepc, before pyexec /////////////////////" << std::endl;
177  call_pylist(3);
178  }
179  }
180 
181  call_pyexec(); // taus: decay pi0's etc; others: decay whatever
182 
183  if(printout_) {
184  std::cout << " ///////////////////// After pyexec, before pyhepc /////////////////////" << std::endl;
185  call_pylist(3);
186  }
187 
188  call_pyhepc(1); // pythia -> hepevt
189 
190  HepMC::IO_HEPEVT conv;
191  evt = std::auto_ptr<HepMC::GenEvent>(new HepMC::GenEvent(*conv.read_next_event()));
192 
193  if(printout_) {
194  evt->print();
195 
196  std::cout << std::endl << "Vertices: " << std::endl;
197  for(HepMC::GenEvent::vertex_const_iterator iter = evt->vertices_begin(); iter != evt->vertices_end(); ++iter) {
198  std::cout << "Vertex " << (*iter)->id() << ", barcode " << (*iter)->barcode() << std::endl;
199 
200  HepMC::ThreeVector point = (*iter)->point3d();
201  std::cout << " Point (" << point.x() << ", " << point.y() << ", " << point.z() << ")" << std::endl;
202 
203  std::cout << " Incoming particles: ";
204  for(HepMC::GenVertex::particles_in_const_iterator pi = (*iter)->particles_in_const_begin(); pi != (*iter)->particles_in_const_end(); ++pi) {
205  std::cout << (*pi)->pdg_id() << " ";
206  }
207  std::cout << std::endl;
208  std::cout << " Outgoing particles: ";
209  for(HepMC::GenVertex::particles_out_const_iterator pi = (*iter)->particles_out_const_begin(); pi != (*iter)->particles_out_const_end(); ++pi) {
210  std::cout << (*pi)->pdg_id() << " ";
211  }
212  std::cout << std::endl << std::endl;
213  }
214  }
215 
216  return evt;
217 }
218 
219 void ParticleReplacerParticleGun::correctTauMass(const reco::MuonCollection& muons, std::vector<HepMC::FourVector>& corrected) {
220  if(abs(gunParticle_) == 15) {
221  // Correct energy for tau
222  for(reco::MuonCollection::const_iterator iMuon = muons.begin(); iMuon != muons.end(); ++iMuon) {
223  const reco::Muon::LorentzVector& vec = iMuon->p4();
224 
225  double E = sqrt(tauMass*tauMass + vec.x()*vec.x() + vec.y()*vec.y() + vec.z()*vec.z());
226 
227  corrected.push_back(HepMC::FourVector(vec.x(), vec.y(), vec.z(), E));
228  }
229  }
230  else {
231  // Just copy the LorentzVector
232  for(reco::MuonCollection::const_iterator iMuon = muons.begin(); iMuon != muons.end(); ++iMuon) {
233  corrected.push_back(iMuon->p4());
234  }
235  }
236 }
237 
239  if(forceTauDecay_ == "" || forceTauDecay_ == "none") return;
240 
241  // for documentation look at Tauola file tauola_photos_ini.f
242  // COMMON block COMMON / TAUBRA / GAMPRT(30),JLIST(30),NCHAN
243 
244  if(forceTauDecay_ == "hadrons"){
245  /* FIXME
246  taubra.gamprt[0] = 0; // disable branchings to electrons in Tauola
247  taubra.gamprt[1] = 0; // disable branchings to muons in Tauola
248  */
249  }
250  else if(forceTauDecay_ == "1prong"){
251  /* FIXME
252  taubra.gamprt[0] = 0; // disable branchings to electrons in Tauola
253  taubra.gamprt[1] = 0; // disable branchings to muons in Tauola
254  taubra.gamprt[7] = 0;
255  taubra.gamprt[9] = 0;
256  taubra.gamprt[10] = 0;
257  taubra.gamprt[11] = 0;
258  taubra.gamprt[12] = 0;
259  taubra.gamprt[13] = 0;
260  taubra.gamprt[17] = 0;
261  */
262  }
263  else if(forceTauDecay_ == "3prong"){
264  /* FIXME
265  taubra.gamprt[0] = 0; // disable branchings to electrons in Tauola
266  taubra.gamprt[1] = 0; // disable branchings to muons in Tauola
267  taubra.gamprt[2] = 0;
268  taubra.gamprt[3] = 0;
269  taubra.gamprt[4] = 0;
270  taubra.gamprt[5] = 0;
271  taubra.gamprt[6] = 0;
272  taubra.gamprt[8] = 0;
273  taubra.gamprt[14] = 0;
274  taubra.gamprt[15] = 0;
275  taubra.gamprt[16] = 0;
276  taubra.gamprt[18] = 0;
277  taubra.gamprt[19] = 0;
278  taubra.gamprt[20] = 0;
279  taubra.gamprt[21] = 0;
280  */
281  }
282  else
283  edm::LogError("MuonReplacement") << "[ParticleReplacerAlgoParticleGun::forceTauoladecays] "
284  << "Unknown value for forcing tau decays: " << forceTauDecay_ << std::endl;
285 }
286 
287 void ParticleReplacerParticleGun::tauola_forParticleGun(int tau_idx, int pdg_id, const HepMC::FourVector& particle_momentum) {
288  if(abs(pdg_id) != 15) {
289  edm::LogError("MuonReplacement") << "[ParticleReplacerAlgoParticleGun::tauola_forParticleGuns] "
290  << "Trying to decay something else than tau: pdg_id = " << pdg_id << std::endl;
291  return;
292  }
293 
294  // By default the index of tau+ is 3 and tau- is 4. The TAUOLA
295  // routine takes internally care of finding the correct
296  // position of the tau which is decayed. However, since we are
297  // calling DEXAY routine directly by ourselves, we must set
298  // the index manually by ourselves. Fortunately, this seems to
299  // be simple.
300  /* FIXME
301  if(printout_)
302  std::cout << " Tauola taupos common block: np1 " << taupos.np1 << " np2 " << taupos.np2 << std::endl;
303  taupos.np1 = tau_idx;
304  taupos.np2 = tau_idx;
305  if(printout_)
306  std::cout << " Resetting taupos common block to: np1 " << taupos.np1 << " np2 " << taupos.np2 << std::endl;
307  */
308 
309  if(pdg_id == -15){ // tau+
310 
311  pol1_[2] = tauHelicity(pdg_id);
312 
313  /* FIXME
314  momdec.p1[0] = particle_momentum.x();
315  momdec.p1[1] = particle_momentum.y();
316  momdec.p1[2] = particle_momentum.z();
317  momdec.p1[3] = particle_momentum.e();
318 
319  momdec.p2[0] = -momdec.p1[0];
320  momdec.p2[1] = -momdec.p1[1];
321  momdec.p2[2] = -momdec.p1[2];
322  momdec.p2[3] = momdec.p1[3];
323 
324  // "mother" p4
325  momdec.q1[0] = 0;
326  momdec.q1[1] = 0;
327  momdec.q1[2] = 0;
328  momdec.q1[3] = momdec.p1[3] + momdec.p2[3];
329 
330  call_dexay(1,pol1);
331  */
332  }
333  else if (pdg_id == 15){ // tau-
334 
335  pol2_[2] = tauHelicity(pdg_id);
336 
337  /* FIXME
338  momdec.p2[0] = particle_momentum.x();
339  momdec.p2[1] = particle_momentum.y();
340  momdec.p2[2] = particle_momentum.z();
341  momdec.p2[3] = particle_momentum.e();
342 
343  momdec.p1[0] = -momdec.p2[0];
344  momdec.p1[1] = -momdec.p2[1];
345  momdec.p1[2] = -momdec.p2[2];
346  momdec.p1[3] = momdec.p2[3];
347 
348  // "mother" p4
349  momdec.q1[0] = 0;
350  momdec.q1[1] = 0;
351  momdec.q1[2] = 0;
352  momdec.q1[3] = momdec.p1[3] + momdec.p2[3];
353 
354  call_dexay(2,pol2);
355  */
356  }
357 }
358 
360  /* tau polarization summary from Tauola source code:
361  in all cases W : pol1 = -1, pol2 = -1 (tau or neutrino)
362  H+: pol1 = +1, pol2 = +1 (tau or neutrino)
363 
364  if neutral higgs : pol1 = +1, pol2 = -1 OR pol1 = -1, pol2 = +1 (tau tau)
365  if Z or undetermined: pol1 = +1, pol2 = +1 OR pol1 = -1, pol2 = -1 (tau tau)
366  */
367 
368  if(pdg_id < 0) { // tau+
369  if(forceTauPlusHelicity_ != 0) {
370  return forceTauPlusHelicity_;
371  }
372  if(forceTauPolarization_ == "W") return -1;
373  if(forceTauPolarization_ == "H+") return 1;
374  if(pol2_[2] != 0) {
375  if(forceTauPolarization_ == "h" ||
376  forceTauPolarization_ == "H" ||
377  forceTauPolarization_ == "A") return -pol2_[2];
378  else return pol2_[2];
379  }
380  else {
381  return randomPolarization(); //all other cases random, when first tau
382  }
383  }
384  else { // tau-
385  if(forceTauMinusHelicity_ != 0) {
386  return forceTauMinusHelicity_;
387  }
388  if(forceTauPolarization_ == "W") return -1;
389  if(forceTauPolarization_ == "H+") return 1;
390  if(pol1_[2] != 0){
391  if(forceTauPolarization_ == "h" ||
392  forceTauPolarization_ == "H" ||
393  forceTauPolarization_ == "A") return -pol1_[2];
394  else return pol1_[2];
395  }
396  else {
397  return randomPolarization(); //all other cases random, when first tau
398  }
399  }
400 
401  edm::LogError("MuonReplacement") << "[ParticleReplacerAlgoParticleGun::tauHelicity] "
402  << "tauHelicity undetermined, returning 0" << std::endl;
403  return 0;
404 }
405 
407  uint32_t randomNumber = rand();
408  if(randomNumber%2 > 0.5) return 1;
409  return -1;
410 }
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
static HepMC::IO_HEPEVT conv
double y() const
y coordinate
Definition: Vertex.h:97
void call_pyexec()
Definition: extraPythia.cc:3
void setPSet(const edm::ParameterSet &)
void correctTauMass(const reco::MuonCollection &muons, std::vector< HepMC::FourVector > &corrected)
#define abs(x)
Definition: mlp_lapack.h:159
#define NULL
Definition: scimark2.h:8
virtual double vy() const
y coordinate of vertex position
double charge(const std::vector< uint8_t > &Ampls)
float float float z
void tauola_forParticleGun(int tau_idx, int pdg_id, const HepMC::FourVector &particle_momentum)
std::vector< Muon > MuonCollection
collection of Muon objects
Definition: MuonFwd.h:9
void call_pylist(int mode)
T sqrt(T t)
Definition: SSEVec.h:48
ParticleReplacerParticleGun(const edm::ParameterSet &, bool)
double z() const
y coordinate
Definition: Vertex.h:99
virtual double vz() const
z coordinate of vertex position
double x() const
x coordinate
Definition: Vertex.h:95
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Candidate.h:41
tuple muons
Definition: patZpeak.py:38
void call_py1ent(int, int, double, double, double)
Definition: extraPythia.cc:4
void tauola_(int *, int *)
Signal rand(Signal arg)
Definition: vlib.cc:442
std::auto_ptr< HepMC::GenEvent > produce(const reco::MuonCollection &muons, const reco::Vertex *pvtx=0, const HepMC::GenEvent *genEvt=0)
static unsigned int const shift
virtual double vx() const
x coordinate of vertex position
tuple cout
Definition: gather_cfg.py:121
double pi
Definition: DDAxes.h:10
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5