CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Hector.cc
Go to the documentation of this file.
1 #include "FWCore/MessageLogger/interface/MessageLogger.h"
2 #include "FWCore/ParameterSet/interface/FileInPath.h"
3 #include "FWCore/ServiceRegistry/interface/Service.h"
4 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
5 #include "IOMC/RandomEngine/src/TRandomAdaptor.h"
6 
7 #include "SimTransport/HectorProducer/interface/Hector.h"
8 
9 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
10 #include "CLHEP/Units/GlobalSystemOfUnits.h"
11 #include "CLHEP/Units/GlobalPhysicalConstants.h"
12 #include "HepMC/SimpleVector.h"
13 
14 #include "H_Parameters.h"
15 
16 #include <math.h>
17 
18 Hector::Hector(const edm::ParameterSet & param, bool verbosity, bool FP420Transport,bool ZDCTransport) :
19  m_verbosity(verbosity),
20  m_FP420Transport(FP420Transport),
21  m_ZDCTransport(ZDCTransport),
22  rootEngine_(0)
23 {
24 
25  // Create LHC beam line
26  edm::ParameterSet hector_par = param.getParameter<edm::ParameterSet>("Hector");
27 
28  // User definitons
29  lengthfp420 = hector_par.getParameter<double>("BeamLineLengthFP420" );
30  m_rpp420_f = (float) hector_par.getParameter<double>("RP420f");
31  m_rpp420_b = (float) hector_par.getParameter<double>("RP420b");
32  lengthzdc = hector_par.getParameter<double>("BeamLineLengthZDC" );
33  lengthd1 = hector_par.getParameter<double>("BeamLineLengthD1" );
34  beam1filename = hector_par.getParameter<string>("Beam1");
35  beam2filename = hector_par.getParameter<string>("Beam2");
36  m_rppzdc = (float) lengthzdc ;
37  m_rppd1 = (float) lengthd1 ;
38  m_smearAng = hector_par.getParameter<bool>("smearAng");
39  m_sigmaSTX = hector_par.getParameter<double>("sigmaSTX" );
40  m_sigmaSTY = hector_par.getParameter<double>("sigmaSTY" );
41  m_smearE = hector_par.getParameter<bool>("smearEnergy");
42  m_sig_e = hector_par.getParameter<double>("sigmaEnergy");
43  etacut = hector_par.getParameter<double>("EtaCutForHector" );
44 
45  edm::Service<edm::RandomNumberGenerator> rng;
46  if ( ! rng.isAvailable() ) {
47  throw cms::Exception("Configuration")
48  << "LHCTransport (Hector) requires the RandomNumberGeneratorService\n"
49  "which is not present in the configuration file. You must add the service\n"
50  "in the configuration file or remove the modules that require it.";
51  }
52  if ( (rng->getEngine()).name() == "TRandom3" ) {
53  rootEngine_ = ( (edm::TRandomAdaptor*) &(rng->getEngine()) )->getRootEngine();
54  if ( m_verbosity ) LogDebug("HectorSetup") << "LHCTransport seed = " << rootEngine_->GetSeed();
55  }
56  else {
57  edm::LogError("WrongRandomNumberGenerator") << "The TRandom3 engine must be used, Random Number Generator Service not correctly initialized!";
58  rootEngine_ = new TRandom3();
59  }
60 
61  theCorrespondenceMap.clear();
62 
63  if(m_verbosity) {
64  edm::LogInfo("HectorSetup") << "===================================================================\n"
65  << " * * * * * * * * * * * * * * * * * * * * * * * * * * * * \n"
66  << " * * \n"
67  << " * --<--<-- A fast simulator --<--<-- * \n"
68  << " * | --<--<-- of particle --<--<-- * \n"
69  << " * ----HECTOR----< * \n"
70  << " * | -->-->-- transport through-->-->-- * \n"
71  << " * -->-->-- generic beamlines -->-->-- * \n"
72  << " * * \n"
73  << " * JINST 2:P09005 (2007) * \n"
74  << " * X Rouby, J de Favereau, K Piotrzkowski (CP3) * \n"
75  << " * http://www.fynu.ucl.ac.be/hector.html * \n"
76  << " * * \n"
77  << " * Center for Cosmology, Particle Physics and Phenomenology * \n"
78  << " * Universite catholique de Louvain * \n"
79  << " * Louvain-la-Neuve, Belgium * \n"
80  << " * * \n"
81  << " * * * * * * * * * * * * * * * * * * * * * * * * * * * * \n"
82  << " Hector configuration: \n"
83  << " m_FP420Transport = " << m_FP420Transport << "\n"
84  << " m_ZDCTransport = " << m_ZDCTransport << "\n"
85  << " lengthfp420 = " << lengthfp420 << "\n"
86  << " m_rpp420_f = " << m_rpp420_f << "\n"
87  << " m_rpp420_b = " << m_rpp420_b << "\n"
88  << " lengthzdc = " << lengthzdc << "\n"
89  << " lengthd1 = " << lengthd1 << "\n"
90  << "===================================================================\n";
91  }
92  edm::FileInPath b1(beam1filename.c_str());
93  edm::FileInPath b2(beam2filename.c_str());
94 
95  // construct beam line for FP420: .
96  if(m_FP420Transport && lengthfp420>0. ) {
97  m_beamlineFP4201 = new H_BeamLine( 1, lengthfp420 + 0.1 ); // (direction, length)
98  m_beamlineFP4202 = new H_BeamLine( -1, lengthfp420 + 0.1 ); //
99  m_beamlineFP4201->fill( b1.fullPath(), 1, "IP5" );
100  m_beamlineFP4202->fill( b2.fullPath(), -1, "IP5" );
101  m_beamlineFP4201->offsetElements( 120, -0.097 );
102  m_beamlineFP4202->offsetElements( 120, +0.097 );
103  m_beamlineFP4201->calcMatrix();
104  m_beamlineFP4202->calcMatrix();
105  }
106  else{
107  if ( m_verbosity ) LogDebug("HectorSetup") << "Hector: WARNING: lengthfp420= " << lengthfp420;
108  }
109 
110 
111  if (m_ZDCTransport && lengthzdc>0. && lengthd1>0.) {
112  // construct beam line for ZDC: .
113  m_beamlineZDC1 = new H_BeamLine( 1, lengthzdc + 0.1 ); // (direction, length)
114  m_beamlineZDC2 = new H_BeamLine( -1, lengthzdc + 0.1 ); //
115  m_beamlineZDC1->fill( b1.fullPath(), 1, "IP5" );
116  m_beamlineZDC2->fill( b2.fullPath(), -1, "IP5" );
117  m_beamlineZDC1->offsetElements( 120, -0.097 );
118  m_beamlineZDC2->offsetElements( 120, +0.097 );
119  m_beamlineZDC1->calcMatrix();
120  m_beamlineZDC2->calcMatrix();
121 
122 
123  // construct beam line for D1: .
124  m_beamlineD11 = new H_BeamLine( 1, lengthd1 + 0.1 ); // (direction, length)
125  m_beamlineD12 = new H_BeamLine( -1, lengthd1 + 0.1 ); //
126  m_beamlineD11->fill( b1.fullPath(), 1, "IP5" );
127  m_beamlineD12->fill( b2.fullPath(), -1, "IP5" );
128  m_beamlineD11->offsetElements( 120, -0.097 );
129  m_beamlineD12->offsetElements( 120, +0.097 );
130  m_beamlineD11->calcMatrix();
131  m_beamlineD12->calcMatrix();
132  }
133  else{
134  if ( m_verbosity ) LogDebug("HectorSetup") << "Hector: WARNING: lengthzdc= " << lengthzdc << "lengthd1= " << lengthd1;
135  }
136 
137 }
138 
139 Hector::~Hector(){
140 
141  for (std::map<unsigned int,H_BeamParticle*>::iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it ) {
142  delete (*it).second;
143  }
144 
145  delete m_beamlineFP4201;
146  delete m_beamlineFP4202;
147  delete m_beamlineZDC1;
148  delete m_beamlineZDC2;
149  delete m_beamlineD11;
150  delete m_beamlineD12;
151 
152 }
153 
154 void Hector::clearApertureFlags(){
155  m_isStoppedfp420.clear();
156  m_isStoppedzdc.clear();
157  m_isStoppedd1.clear();
158 }
159 
160 void Hector::clear(){
161  for ( std::map<unsigned int,H_BeamParticle*>::iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it ) {
162  delete (*it).second;
163  };
164  m_beamPart.clear();
165  m_direct.clear();
166  m_eta.clear();
167  m_pdg.clear();
168  m_pz.clear();
169  m_isCharged.clear();
170 }
171 
172 /*
173  bool Hector::isCharged(const HepMC::GenParticle * p){
174  const ParticleData * part = pdt->particle( p->pdg_id() );
175  if (part){
176  return part->charge()!=0;
177  }else{
178  // the new/improved particle table doesn't know anti-particles
179  return pdt->particle( -p->pdg_id() )!=0;
180  }
181  }
182 */
183 void Hector::add( const HepMC::GenEvent * evt ,const edm::EventSetup & iSetup) {
184 
185  H_BeamParticle * h_p;
186  unsigned int line;
187 
188  for (HepMC::GenEvent::particle_const_iterator eventParticle =evt->particles_begin();
189  eventParticle != evt->particles_end();
190  ++eventParticle ) {
191  if ( (*eventParticle)->status() == 1 ) {
192  if ( abs( (*eventParticle)->momentum().eta())>etacut){
193  line = (*eventParticle)->barcode();
194  if ( m_beamPart.find(line) == m_beamPart.end() ) {
195  double charge=1.;
196  m_isCharged[line] = false;// neutrals
197  HepMC::GenParticle * g = (*eventParticle);
198  iSetup.getData( pdt );
199  const ParticleData * part = pdt->particle( g->pdg_id() );
200  if (part){
201  charge = part->charge();
202  }
203  if(charge !=0) m_isCharged[line] = true;//charged
204  double mass = (*eventParticle)->generatedMass();
205 
206  h_p = new H_BeamParticle(mass,charge);
207 
208  double px,py,pz;
209  px = (*eventParticle)->momentum().px();
210  py = (*eventParticle)->momentum().py();
211  pz = (*eventParticle)->momentum().pz();
212 
213  h_p->set4Momentum( px, py, pz, (*eventParticle)->momentum().e() );
214 
215  // from mm to um
216  double XforPosition = (*eventParticle)->production_vertex()->position().x()/micrometer;//um
217  double YforPosition = (*eventParticle)->production_vertex()->position().y()/micrometer;//um
218  double ZforPosition = (*eventParticle)->production_vertex()->position().z()/meter;//m
219  // crossing angle (beam tilt) is not known a priory; keep now 0.0 but, in principle, can be entered as parameters
220  double TXforPosition=0.0, TYforPosition=0.0;//urad
221 
222  // Clears H_BeamParticle::positions and sets the initial one
223  h_p->setPosition(XforPosition, YforPosition, TXforPosition, TYforPosition, ZforPosition );
224 
225  m_beamPart[line] = h_p;
226  m_direct[line] = 0;
227  m_direct[line] = ( pz > 0 ) ? 1 : -1;
228 
229  m_eta[line] = (*eventParticle)->momentum().eta();
230  m_pdg[line] = (*eventParticle)->pdg_id();
231  m_pz[line] = (*eventParticle)->momentum().pz();
232 
233  if(m_verbosity) {
234  LogDebug("HectorEventProcessing") << "Hector:add: barcode = " << line
235  << " status = " << g->status()
236  << " PDG Id = " << g->pdg_id()
237  << " mass = " << mass
238  << " pz = " << pz
239  << " charge = " << charge
240  << " m_isCharged[line] = " << m_isCharged[line];
241  }
242  }// if find line
243  }// if eta > 8.2
244  }// if status
245  }// for loop
246 
247 }
248 
249 void Hector::filterFP420(){
250  unsigned int line;
251  H_BeamParticle * part;
252  std::map< unsigned int, H_BeamParticle* >::iterator it;
253 
254  bool is_stop;
255  int direction;
256 
257  float x1_420;
258  float y1_420;
259 
260  if ( m_beamPart.size() && lengthfp420>0. ) {
261 
262  for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it ) {
263  line = (*it).first;
264  part = (*it).second;
265 
266  if(m_verbosity) {
267  LogDebug("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line;
268  }
269  if ( (*m_isCharged.find( line )).second ) {
270  direction = (*m_direct.find( line )).second;
271  if (m_smearAng) {
272  // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
273  if ( m_sigmaSTX>0. && m_sigmaSTY>0.) {
274  part->smearAng(m_sigmaSTX,m_sigmaSTY,rootEngine_);
275  }
276  else {
277  // for smearAng() in urad, default are (STX=30.23, STY=30.23)
278  part->smearAng(STX,STY,rootEngine_);
279  }
280  }
281  if (m_smearE) {
282  if ( m_sig_e ) {
283  part->smearE(m_sig_e,rootEngine_);
284  }
285  else {
286  part->smearE(SBE,rootEngine_); // in GeV, default is SBE=0.79
287  }
288  }
289  if ( direction == 1 && m_beamlineFP4201 != 0 ) {
290  part->computePath( m_beamlineFP4201 );
291  is_stop = part->stopped( m_beamlineFP4201 );
292  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line << " positive is_stop= "<< is_stop;
293  }
294  else if ( direction == -1 && m_beamlineFP4202 != 0 ){
295  part->computePath( m_beamlineFP4202 );
296  is_stop = part->stopped( m_beamlineFP4202 );
297  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line << " negative is_stop= "<< is_stop;
298  }
299  else {
300  is_stop = true;
301  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line << " 0 is_stop= "<< is_stop;
302  }
303 
304  //propagating
305  m_isStoppedfp420[line] = is_stop;
306  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line << " isStopped=" << (*m_isStoppedfp420.find(line)).second;
307 
308  if (!is_stop) {
309  if ( direction == 1 ) part->propagate( m_rpp420_f );
310  if ( direction == -1 ) part->propagate( m_rpp420_b );
311  x1_420 = part->getX();
312  y1_420 = part->getY();
313  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line << " x= "<< x1_420 <<" y= " << y1_420;
314 
315  m_xAtTrPoint[line] = x1_420;
316  m_yAtTrPoint[line] = y1_420;
317  m_TxAtTrPoint[line] = part->getTX();
318  m_TyAtTrPoint[line] = part->getTY();
319  m_eAtTrPoint[line] = part->getE();
320 
321  }
322  }// if isCharged
323  else {
324  m_isStoppedfp420[line] = true;// imply that neutral particles stopped to reach 420m
325  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterFP420: barcode = " << line << " isStopped=" << (*m_isStoppedfp420.find(line)).second;
326  }
327 
328  } // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
329  } // if ( m_beamPart.size() )
330 
331 }
332 
333 void Hector::filterZDC(){
334  unsigned int line;
335  H_BeamParticle * part;
336  std::map< unsigned int, H_BeamParticle* >::iterator it;
337 
338  bool is_stop_zdc = false;
339  int direction;
340 
341  if ( m_beamPart.size() && lengthzdc>0. ) {
342 
343  for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it ) {
344  line = (*it).first;
345  part = (*it).second;
346  if(m_verbosity) {
347  LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode = " << line << " charge = " << (*m_isCharged.find(line)).second;
348  if (m_FP420Transport) LogDebug("HectorEventProcessing") << " isStoppedFP420 =" << (*m_isStoppedfp420.find(line)).second;
349  }
350  // if ( ((*m_isStoppedfp420.find(line)).second) && ((*m_isCharged.find(line)).second) ) {
351  if ( ((*m_isCharged.find(line)).second) ) {
352 
353  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode = " << line << " propagated ";
354 
355  direction = (*m_direct.find( line )).second;
356  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode = " << line << " direction = " << direction;
357  if (m_smearAng) {
358  if ( m_sigmaSTX>0. && m_sigmaSTY>0.) {
359  // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
360  part->smearAng(m_sigmaSTX,m_sigmaSTY,rootEngine_);
361  } else {
362  // for smearAng() in urad, default are (STX=30.23, STY=30.23)
363  part->smearAng(STX,STY,rootEngine_);
364  }
365  }
366  if (m_smearE) {
367  if ( m_sig_e ) {
368  part->smearE(m_sig_e,rootEngine_);
369  } else {
370  part->smearE(SBE,rootEngine_); // in GeV, default is SBE=0.79
371  }
372  }
373  if ( direction == 1 && m_beamlineZDC1 != 0 ){
374  part->computePath( m_beamlineZDC1 );
375  is_stop_zdc = part->stopped( m_beamlineZDC1 );
376  m_isStoppedzdc[line] = is_stop_zdc;
377  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode " << line << " positive is_stop_zdc= "<< is_stop_zdc;
378  } else if ( direction == -1 && m_beamlineZDC2 != 0 ){
379  part->computePath( m_beamlineZDC2 );
380  is_stop_zdc = part->stopped( m_beamlineZDC2 );
381  m_isStoppedzdc[line] = is_stop_zdc;
382  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode " << line << " negative is_stop_zdc= "<< is_stop_zdc;
383  } else {
384  m_isStoppedzdc[line] = true;
385  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode " << line << " 0 is_stop_zdc= "<< is_stop_zdc;
386  }
387  }
388  // if stopfp420 charged particles
389  /*
390  else if ( ((*m_isCharged.find(line)).second) ){
391  m_isStoppedzdc[line] = false;// not stopped in propagating to FP420 and therefore in propagation to ZDC too.
392  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode = " << line << " isStopped=" << (*m_isStoppedzdc.find(line)).second;
393  } */
394  else {
395  m_isStoppedzdc[line] = true;// neutrals particles considered as stopped in propagating to ZDC
396  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterZDC: barcode = " << line << " isStopped=" << (*m_isStoppedzdc.find(line)).second;
397  }
398 
399  } // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
400  } // if ( m_beamPart.size() )
401 
402 }
403 
404 void Hector::filterD1(){
405  unsigned int line;
406  H_BeamParticle * part;
407  std::map< unsigned int, H_BeamParticle* >::iterator it;
408 
409  bool is_stop_d1;
410  int direction;
411 
412  float x1_d1;
413  float y1_d1;
414 
415  if ( m_beamPart.size() && lengthd1>0.) {
416 
417  for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it ) {
418  line = (*it).first;
419  part = (*it).second;
420  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterD1: barcode = " << line << " isStoppedZDC =" << (*m_isStoppedzdc.find(line)).second;
421  if ( ((*m_isStoppedzdc.find(line)).second) || !((*m_isCharged.find( line )).second) ) {
422 
423  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterD1: barcode = " << line << " propagated ";
424 
425  direction = (*m_direct.find( line )).second;
426  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterD1:direction=" << direction;
427  if (m_smearAng) {
428  if ( m_sigmaSTX>0. && m_sigmaSTY>0.) {
429  // the beam transverse direction is centered on (TXforPosition, TYforPosition) at IP
430  part->smearAng(m_sigmaSTX,m_sigmaSTY,rootEngine_);
431  } else {
432  // for smearAng() in urad, default are (STX=30.23, STY=30.23)
433  part->smearAng(STX,STY,rootEngine_);
434  }
435  }
436  if (m_smearE) {
437  if ( m_sig_e ) {
438  part->smearE(m_sig_e,rootEngine_);
439  } else {
440  part->smearE(SBE,rootEngine_); // in GeV, default is SBE=0.79
441  }
442  }
443  if ( direction == 1 && m_beamlineD11 != 0 ) {
444  part->computePath( m_beamlineD11 );
445  is_stop_d1 = part->stopped( m_beamlineD11 );
446  m_isStoppedd1[line] = is_stop_d1;
447  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterD1 barcode " << line << " positive is_stop_d1 = "<< is_stop_d1;
448  } else if ( direction == -1 && m_beamlineD12 != 0 ){
449  part->computePath( m_beamlineD12 );
450  is_stop_d1 = part->stopped( m_beamlineD12 );
451  m_isStoppedd1[line] = is_stop_d1;
452  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterD1 barcode " << line << " negative is_stop_d1 = "<< is_stop_d1;
453  } else {
454  is_stop_d1 = true;
455  m_isStoppedd1[line] = is_stop_d1;
456  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterD1 barcode " << line << " 0 is_stop_d1 = "<< is_stop_d1;
457  }
458  //propagating
459  if (!is_stop_d1 ) {
460  part->propagate( lengthd1 );
461  x1_d1 = part->getX();
462  y1_d1 = part->getY();
463  m_xAtTrPoint[line] = x1_d1;
464  m_yAtTrPoint[line] = y1_d1;
465  m_TxAtTrPoint[line] = part->getTX();
466  m_TyAtTrPoint[line] = part->getTY();
467  m_eAtTrPoint[line] = part->getE();
468  }
469  }// if stopzdc
470  else {
471  m_isStoppedd1[line] = false;// not stopped in propagating to ZDC and therefore in propagation to D1 too.
472  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector:filterD1: barcode = " << line << " isStopped=" << (*m_isStoppedd1.find(line)).second;
473  }
474 
475  } // for (it = m_beamPart.begin(); it != m_beamPart.end(); it++ )
476  } // if ( m_beamPart.size() )
477 
478 }
479 
480 int Hector::getDirect( unsigned int part_n ) const {
481  std::map<unsigned int, int>::const_iterator it = m_direct.find( part_n );
482  if ( it != m_direct.end() ){
483  return (*it).second;
484  }
485  return 0;
486 }
487 
488 void Hector::print() const {
489  for (std::map<unsigned int,H_BeamParticle*>::const_iterator it = m_beamPart.begin(); it != m_beamPart.end(); ++it ) {
490  (*it).second->printProperties();
491  };
492 }
493 
494 
495 HepMC::GenEvent * Hector::addPartToHepMC( HepMC::GenEvent * evt ){
496 
497  theCorrespondenceMap.clear();
498 
499  unsigned int line;
500 
501  HepMC::GenParticle * gpart;
502  long double tx,ty,theta,fi,energy,time = 0;
503  std::map< unsigned int, H_BeamParticle* >::iterator it;
504 
505 
506  for (it = m_beamPart.begin(); it != m_beamPart.end(); ++it ) {
507  line = (*it).first;
508  if(!m_FP420Transport) m_isStoppedfp420[line] = true;
509  if(!m_ZDCTransport) {m_isStoppedzdc[line] = false;m_isStoppedd1[line] = true;}
510  if(m_verbosity) {
511  LogDebug("HectorEventProcessing") << "Hector:addPartToHepMC: barcode = " << line << "\n"
512  << "Hector:addPartToHepMC: isStoppedfp420=" << (*m_isStoppedfp420.find(line)).second << "\n"
513  << "Hector:addPartToHepMC: isStoppedzdc=" << (*m_isStoppedzdc.find(line)).second << "\n"
514  << "Hector:addPartToHepMC: isStoppedd1=" << (*m_isStoppedd1.find(line)).second;
515  }
516  if (!((*m_isStoppedfp420.find(line)).second) || (!((*m_isStoppedd1.find(line)).second) && ((*m_isStoppedzdc.find(line)).second))){
517 
518  gpart = evt->barcode_to_particle( line );
519  if ( gpart ) {
520  tx = (*m_TxAtTrPoint.find(line)).second / 1000000.;
521  ty = (*m_TyAtTrPoint.find(line)).second / 1000000.;
522  theta = sqrt((tx*tx) + (ty*ty));
523  double ddd = 0.;
524  if ( !((*m_isStoppedfp420.find(line)).second)) {
525  if( (*m_direct.find( line )).second >0 ) {
526  ddd = m_rpp420_f;
527  }
528  else if((*m_direct.find( line )).second <0 ) {
529  ddd = m_rpp420_b;
530  theta= pi-theta;
531  }
532  }
533  else {
534  ddd = lengthd1;
535  if((*m_direct.find( line )).second <0 ) theta= pi-theta;
536  }
537 
538  fi = std::atan2(tx,ty); // tx, ty never == 0?
539  energy = (*m_eAtTrPoint.find(line)).second;
540 
541  time = ( ddd*meter - gpart->production_vertex()->position().z()*mm ); // mm
542 
543  if(ddd != 0.) {
544  if(m_verbosity) {
545  LogDebug("HectorEventProcessing") <<"Hector:: x= "<< (*(m_xAtTrPoint.find(line))).second*0.001<< "\n"
546  <<"Hector:: y= "<< (*(m_yAtTrPoint.find(line))).second*0.001<< "\n"
547  <<"Hector:: z= "<< ddd * (*(m_direct.find( line ))).second*1000. << "\n"
548  <<"Hector:: t= "<< time;
549  }
550 
551  HepMC::GenVertex * vert = new HepMC::GenVertex( HepMC::FourVector( (*(m_xAtTrPoint.find(line))).second*0.001,
552  (*(m_yAtTrPoint.find(line))).second*0.001,
553  ddd * (*(m_direct.find( line ))).second*1000.,
554  time + .001*time ) );
555 
556  gpart->set_status( 2 );
557  vert->add_particle_in( gpart );
558  vert->add_particle_out( new HepMC::GenParticle( HepMC::FourVector(energy*std::sin(theta)*std::sin(fi),
559  energy*std::sin(theta)*std::cos(fi),
560  energy*std::cos(theta),
561  energy ),
562  gpart->pdg_id(), 1, gpart->flow() ) );
563  evt->add_vertex( vert );
564 
565  int ingoing = (*vert->particles_in_const_begin())->barcode();
566  int outgoing = (*vert->particles_out_const_begin())->barcode();
567  LHCTransportLink theLink(ingoing,outgoing);
568  if (m_verbosity) LogDebug("HectorEventProcessing") << "Hector:addPartToHepMC: LHCTransportLink " << theLink;
569  theCorrespondenceMap.push_back(theLink);
570 
571  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector::TRANSPORTED pz= " << gpart->momentum().pz()
572  << " eta= "<< gpart->momentum().eta()
573  << " status= "<< gpart->status();
574 
575 
576  }// ddd
577  }// if gpart
578  }// if !isStopped
579 
580  else {
581  gpart = evt->barcode_to_particle( line );
582  if ( gpart ) {
583  // HepMC::GenVertex * vert= new HepMC::GenVertex();
584  gpart->set_status( 2 );
585  // vert->add_particle_in( gpart );
586  // vert->add_particle_out( gpart );
587  // evt->add_vertex( vert );
588  if(m_verbosity) LogDebug("HectorEventProcessing") << "Hector::NON-transp. pz= " << gpart->momentum().pz()
589  << " eta= "<< gpart->momentum().eta()
590  << " status= "<< gpart->status();
591  }
592  }
593 
594  }//for
595 
596  return evt;
597 }
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Hector::m_direct
std::map< unsigned int, int > m_direct
Definition: Hector.h:113
Hector::m_isCharged
std::map< unsigned int, bool > m_isCharged
Definition: Hector.h:126
Hector::getDirect
int getDirect(unsigned int part_n) const
Definition: Hector.cc:480
Hector::m_rpp420_f
float m_rpp420_f
Definition: Hector.h:93
funct::sin
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Hector::m_ZDCTransport
bool m_ZDCTransport
Definition: Hector.h:133
theta
Geom::Theta< T > theta() const
Definition: Basic3DVectorLD.h:172
Hector::m_eta
std::map< unsigned int, double > m_eta
Definition: Hector.h:123
Hector::m_rpp420_b
float m_rpp420_b
Definition: Hector.h:94
Hector::pdt
edm::ESHandle< ParticleDataTable > pdt
Definition: Hector.h:98
Hector::m_isStoppedzdc
std::map< unsigned int, bool > m_isStoppedzdc
Definition: Hector.h:115
abs
#define abs(x)
Definition: mlp_lapack.h:159
Hector::filterFP420
void filterFP420()
Definition: Hector.cc:249
Hector::m_verbosity
bool m_verbosity
Definition: Hector.h:131
Hector::m_sig_e
double m_sig_e
Definition: Hector.h:88
Hector::m_beamlineZDC1
H_BeamLine * m_beamlineZDC1
Definition: Hector.h:103
Hector::m_beamlineD11
H_BeamLine * m_beamlineD11
Definition: Hector.h:105
Hector::m_isStoppedd1
std::map< unsigned int, bool > m_isStoppedd1
Definition: Hector.h:116
Hector::m_isStoppedfp420
std::map< unsigned int, bool > m_isStoppedfp420
Definition: Hector.h:114
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
Hector::clear
void clear()
Definition: Hector.cc:160
Hector::m_sigmaSTX
double m_sigmaSTX
Definition: Hector.h:90
g
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
edm::EventSetup::getData
void getData(T &iHolder) const
Definition: EventSetup.h:67
Hector::Hector
Hector(const edm::ParameterSet &ps, bool verbosity, bool FP420Transport, bool ZDCTransport)
Definition: Hector.cc:18
Hector::m_eAtTrPoint
std::map< unsigned int, double > m_eAtTrPoint
Definition: Hector.h:121
edm::second
U second(std::pair< T, U > const &p)
Definition: ParameterSet.cc:234
Hector::m_rppd1
float m_rppd1
Definition: Hector.h:96
PFRecoTauDiscriminationAgainstMuon2_cfi.verbosity
tuple verbosity
Definition: PFRecoTauDiscriminationAgainstMuon2_cfi.py:26
Hector::m_beamlineZDC2
H_BeamLine * m_beamlineZDC2
Definition: Hector.h:104
Hector::m_pz
std::map< unsigned int, double > m_pz
Definition: Hector.h:125
Hector::lengthfp420
double lengthfp420
Definition: Hector.h:82
Hector::beam2filename
string beam2filename
Definition: Hector.h:129
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
funct::cos
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Hector::m_TyAtTrPoint
std::map< unsigned int, double > m_TyAtTrPoint
Definition: Hector.h:120
Hector::theCorrespondenceMap
std::vector< LHCTransportLink > theCorrespondenceMap
Definition: Hector.h:135
Hector::m_sigmaSTY
double m_sigmaSTY
Definition: Hector.h:91
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:10
Hector::m_beamPart
std::map< unsigned int, H_BeamParticle * > m_beamPart
Definition: Hector.h:112
Hector::print
void print() const
Definition: Hector.cc:488
Hector::~Hector
virtual ~Hector()
Definition: Hector.cc:139
Hector::m_beamlineFP4201
H_BeamLine * m_beamlineFP4201
Definition: Hector.h:101
part
part
Definition: HCALResponse.h:20
Hector::m_smearE
bool m_smearE
Definition: Hector.h:89
Hector::etacut
double etacut
Definition: Hector.h:86
Hector::beam1filename
string beam1filename
Definition: Hector.h:128
Hector::add
void add(const HepMC::GenEvent *ev, const edm::EventSetup &es)
Definition: Hector.cc:183
Hector::m_smearAng
bool m_smearAng
Definition: Hector.h:87
Hector::m_beamlineD12
H_BeamLine * m_beamlineD12
Definition: Hector.h:106
Hector::m_beamlineFP4202
H_BeamLine * m_beamlineFP4202
Definition: Hector.h:102
Hector::m_rppzdc
float m_rppzdc
Definition: Hector.h:95
Hector::addPartToHepMC
HepMC::GenEvent * addPartToHepMC(HepMC::GenEvent *event)
Definition: Hector.cc:495
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
Hector::lengthzdc
double lengthzdc
Definition: Hector.h:83
Hector::m_yAtTrPoint
std::map< unsigned int, double > m_yAtTrPoint
Definition: Hector.h:118
pi
double pi
Definition: Vx3DHLTAnalyzer.h:63
Hector::m_pdg
std::map< unsigned int, int > m_pdg
Definition: Hector.h:124
Hector::m_TxAtTrPoint
std::map< unsigned int, double > m_TxAtTrPoint
Definition: Hector.h:119
Hector::m_xAtTrPoint
std::map< unsigned int, double > m_xAtTrPoint
Definition: Hector.h:117
Hector::lengthd1
double lengthd1
Definition: Hector.h:84
Hector::rootEngine_
TRandom3 * rootEngine_
Definition: Hector.h:137
Hector::filterD1
void filterD1()
Definition: Hector.cc:404
Hector::m_FP420Transport
bool m_FP420Transport
Definition: Hector.h:132
Hector::filterZDC
void filterZDC()
Definition: Hector.cc:333
Hector::clearApertureFlags
void clearApertureFlags()
Definition: Hector.cc:154
Generated for CMSSW Reference Manual by doxygen 1.8.5