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