CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SimTrackManager.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: Application
4 // Class : SimTrackManager
5 //
6 // Implementation:
7 // <Notes on implementation>
8 //
9 // Original Author:
10 // Created: Fri Nov 25 17:44:19 EST 2005
11 //
12 
13 // system include files
14 #include <iostream>
15 
16 // user include files
20 
22 
23 #include "G4VProcess.hh"
24 
25 //#define DebugLog
26 //using namespace std;
27 
28 //
29 // constants, enums and typedefs
30 //
31 
32 //
33 // static data member definitions
34 //
35 
36 //
37 // constructors and destructor
38 //
39 SimTrackManager::SimTrackManager(bool iCollapsePrimaryVertices) :
40  m_trksForThisEvent(0),m_nVertices(0),
41  m_collapsePrimaryVertices(iCollapsePrimaryVertices),
42  lastTrack(0),lastHist(0),theLHCTlink(0){}
43 
44 
46 {
47  if ( m_trksForThisEvent != 0 ) deleteTracks() ;
48 }
49 
50 //
51 // member functions
52 //
54 {
56  else
57  {
58  for (unsigned int i = 0; i < m_trksForThisEvent->size(); i++) {
59  delete (*m_trksForThisEvent)[i];
60  }
61  delete m_trksForThisEvent;
63  }
66  std::vector<std::pair <int, int> >().swap(idsave);
67  ancestorList.clear();
68  lastTrack=0;
69  lastHist=0;
70 }
71 
73 {
74  for (unsigned int i = 0; i < m_trksForThisEvent->size(); i++) {
75  delete (*m_trksForThisEvent)[i];
76  }
77  delete m_trksForThisEvent;
79 }
80 
83 {
84  using namespace std;
85  TrackWithHistory * trkH = trkWHist;
86  if (trkH == 0)
87  {
88  edm::LogError("SimTrackManager")
89  << " SimTrackManager::saveTrackAndItsBranch got 0 pointer ";
90  abort();
91  }
92  trkH->save();
93  unsigned int parent = trkH->parentID();
94 
95  TrackContainer::const_iterator tk_itr =
96  std::lower_bound((*m_trksForThisEvent).begin(),(*m_trksForThisEvent).end(),
98 
99  TrackWithHistory * tempTk = *tk_itr;
100 
101  if (tk_itr!=m_trksForThisEvent->end() && (*tk_itr)->trackID()==parent) {
102  saveTrackAndItsBranch(tempTk);
103  }
104 }
105 
107 {
109 
110  // fill the map with the final mother-daughter relationship
111  idsave.swap(ancestorList);
112  stable_sort(idsave.begin(),idsave.end());
113 
114  std::vector<std::pair<int,int> >().swap(ancestorList);
115 
116  // to get a backward compatible order
117  stable_sort(m_trksForThisEvent->begin(),m_trksForThisEvent->end(),trkIDLess());
118 
119  // to reset the GenParticle ID of a SimTrack to its pre-LHCTransport value
120  resetGenID();
121 
122  reallyStoreTracks(simEvent);
123 }
124 
126 {
127  // loop over the (now ordered) vector and really save the tracks
128 #ifdef DebugLog
129  LogDebug("SimTrackManager")
130  << "Inside the reallyStoreTracks method object to be stored = "
131  << m_trksForThisEvent->size();
132 #endif
133 
134  for (unsigned int it = 0; it < m_trksForThisEvent->size(); it++)
135  {
136  TrackWithHistory * trkH = (*m_trksForThisEvent)[it];
137  // at this stage there is one vertex per track,
138  // so the vertex id of track N is also N
139  int ivertex = -1;
140  int ig;
141 
142  math::XYZVectorD pm(0.,0.,0.);
143  unsigned int iParentID = trkH->parentID();
144  for(unsigned int iit = 0; iit < m_trksForThisEvent->size(); ++iit)
145  {
146  if((*m_trksForThisEvent)[iit]->trackID()==iParentID){
147  pm = (*m_trksForThisEvent)[iit]->momentum();
148  break;
149  }
150  }
151  ig = trkH->genParticleID();
152  ivertex = getOrCreateVertex(trkH,iParentID,simEvent);
153  std::map<uint32_t,std::pair<math::XYZVectorD,math::XYZTLorentzVectorD> >::const_iterator cit =
154  mapTkCaloStateInfo.find(trkH->trackID());
155  std::pair<math::XYZVectorD,math::XYZTLorentzVectorD> tcinfo;
156  if (cit != mapTkCaloStateInfo.end()){
157  tcinfo = cit->second;
158  }
159  simEvent->add(new G4SimTrack(trkH->trackID(),trkH->particleID(),
160  trkH->momentum(),trkH->totalEnergy(),
161  ivertex,ig,pm,tcinfo.first,tcinfo.second));
162  }
163 }
164 
166  G4SimEvent * simEvent)
167 {
168  int parent = iParentID;
169  int check = -1;
170 
171  for( std::vector<TrackWithHistory*>::const_iterator it = (*m_trksForThisEvent).begin();
172  it!= (*m_trksForThisEvent).end();it++)
173  {
174  if ((*it)->trackID() == uint32_t(parent)) {
175  check = 0;
176  break;
177  }
178  }
179 
180  if(check==-1) { parent = -1; }
181 
182  VertexMap::const_iterator iterator = m_vertexMap.find(parent);
183  if (iterator != m_vertexMap.end()) {
184 
185  // loop over saved vertices
186  for(unsigned int k=0; k<m_vertexMap[parent].size(); k++){
187  if(sqrt((trkH->vertexPosition()-(((m_vertexMap[parent])[k]).second)).Mag2())<0.001)
188  { return (((m_vertexMap[parent])[k]).first); }
189  }
190  }
191 
192  unsigned int ptype = 0;
193  const G4VProcess* pr = trkH->creatorProcess();
194  if(pr) { ptype = pr->GetProcessSubType(); }
195  simEvent->add(new G4SimVertex(trkH->vertexPosition(),trkH->globalTime(),parent,ptype));
197  m_nVertices++;
198  return (m_nVertices-1);
199 
200 }
201 
203 {
204  m_vertexMap.clear();
206  m_nVertices=0;
207 }
208 
210 {
211  mapTkCaloStateInfo.clear();
212  std::map<uint32_t,std::pair<math::XYZVectorD,math::XYZTLorentzVectorD > >().swap(mapTkCaloStateInfo);
213 }
214 
216 {
217  int idMother = id;
218  if(id > 0) {
219  unsigned int n = idsave.size();
220  if(0 < n) {
221  int jmax = n - 1;
222  int j, id1;
223 
224  // first loop forward
225  bool notFound = true;
226  for(j=0; j<=jmax; ++j) {
227  if((idsave[j]).first == idMother) {
228  id1 = (idsave[j]).second;
229  if(0 == id1 || id1 == idMother) { return id1; }
230  jmax = j - 1;
231  idMother = id1;
232  notFound = false;
233  break;
234  }
235  }
236  if(notFound) { return 0; }
237 
238  // recursive loop
239  do {
240 
241  notFound = true;
242  // search ID scan backward
243  for(j=jmax; j>=0; --j) {
244  if((idsave[j]).first == idMother) {
245  id1 = (idsave[j]).second;
246  if(0 == id1 || id1 == idMother) { return id1; }
247  jmax = j - 1;
248  idMother = id1;
249  notFound = false;
250  break;
251  }
252  }
253  if(notFound) {
254  // ID not in the list of saved track - look into ancestors
255  jmax = ancestorList.size()-1;
256  for(j=jmax; j>=0; --j) {
257  if((ancestorList[j]).first == idMother) {
258  idMother = (ancestorList[j]).second;
259  return idMother;
260  }
261  }
262  return 0;
263  }
264  } while (!notFound);
265  }
266  }
267  return idMother;
268 }
269 
271 {
272  if ( ancestorList.size() > 0 && lastHist > ancestorList.size() ) {
273  lastHist = ancestorList.size();
274  edm::LogError("SimTrackManager")
275  << " SimTrackManager::fillMotherList track index corrupted";
276  }
277  /*
278  std::cout << "### SimTrackManager::fillMotherList: "
279  << idsave.size() << " saved; ancestor: " << lastHist
280  << " " << ancestorList.size() << std::endl;
281  for (unsigned int i = 0; i< idsave.size(); ++i) {
282  std::cout << " ISV: Track ID = " << (idsave[i]).first
283  << " Mother ID = " << (idsave[i]).second << std::endl;
284  }
285  */
286  for (unsigned int n = lastHist; n < ancestorList.size(); n++) {
287 
288  int theMotherId = idSavedTrack((ancestorList[n]).first);
289  ancestorList[n].second = theMotherId;
290  /*
291  std::cout << " ANC: Track ID = " << (ancestorList[n]).first
292  << " Mother ID = " << (ancestorList[n]).second << std::endl;
293  */
294 #ifdef DebugLog
295  LogDebug("SimTrackManager") << "Track ID = " << (ancestorList[n]).first
296  << " Mother ID = " << (ancestorList[n]).second;
297 #endif
298  }
299 
300  lastHist = ancestorList.size();
301 
302  idsave.clear();
303 
304 }
305 
307 
308  if ((*m_trksForThisEvent).size() == 0 && idsave.size() == 0) { return; }
309 
310 #ifdef DebugLog
311  LogDebug("SimTrackManager")
312  << "SimTrackManager::cleanTracksWithHistory has "
313  << idsave.size()
314  << " mother-daughter relationships stored with lastTrack = " << lastTrack;
315 #endif
316 
317  if ( lastTrack > 0 && lastTrack >= (*m_trksForThisEvent).size() ) {
318  lastTrack = 0;
319  edm::LogError("SimTrackManager")
320  << " SimTrackManager::cleanTracksWithHistory track index corrupted";
321  }
322 
323  stable_sort(m_trksForThisEvent->begin()+lastTrack,m_trksForThisEvent->end(),trkIDLess());
324 
325  stable_sort(idsave.begin(),idsave.end());
326 
327 #ifdef DebugLog
328  LogDebug("SimTrackManager")
329  << " SimTrackManager::cleanTracksWithHistory knows " << m_trksForThisEvent->size()
330  << " tracks with history before branching";
331  for (unsigned int it =0; it <(*m_trksForThisEvent).size(); it++) {
332  LogDebug("SimTrackManager")
333  << " 1 - Track in position " << it << " G4 track number "
334  << (*m_trksForThisEvent)[it]->trackID()
335  << " mother " << (*m_trksForThisEvent)[it]->parentID()
336  << " status " << (*m_trksForThisEvent)[it]->saved();
337  }
338 #endif
339 
340  for (unsigned int it = lastTrack; it < m_trksForThisEvent->size(); it++)
341  {
342  TrackWithHistory * t = (*m_trksForThisEvent)[it];
343  if (t->saved()) { saveTrackAndItsBranch(t); }
344  }
345  unsigned int num = lastTrack;
346  for (unsigned int it = lastTrack; it < m_trksForThisEvent->size(); it++)
347  {
348  TrackWithHistory * t = (*m_trksForThisEvent)[it];
349  int g4ID = t->trackID();
350  if (t->saved() == true)
351  {
352  if (it>num) (*m_trksForThisEvent)[num] = t;
353  num++;
354  for (unsigned int itr=0; itr<idsave.size(); itr++) {
355  if ((idsave[itr]).first == g4ID) {
356  (idsave[itr]).second = g4ID;
357  break;
358  }
359  }
360  }
361  else
362  {
363  delete t;
364  }
365  }
366 
367  (*m_trksForThisEvent).resize(num);
368 
369 #ifdef DebugLog
370  LogDebug("SimTrackManager")
371  << " AFTER CLEANING, I GET " << (*m_trksForThisEvent).size()
372  << " tracks to be saved persistently";
373  for (unsigned int it = 0; it < (*m_trksForThisEvent).size(); it++) {
374  LogDebug("SimTrackManager")
375  << " Track in position " << it
376  << " G4 track number " << (*m_trksForThisEvent)[it]->trackID()
377  << " mother " << (*m_trksForThisEvent)[it]->parentID()
378  << " Status " << (*m_trksForThisEvent)[it]->saved()
379  << " id " << (*m_trksForThisEvent)[it]->particleID()
380  << " E(MeV)= " << (*m_trksForThisEvent)[it]->totalEnergy();
381  }
382 #endif
383 
384  fillMotherList();
385 
386  lastTrack = (*m_trksForThisEvent).size();
387 }
388 
390 {
391  if ( theLHCTlink == 0 ) return;
392 
393  for (unsigned int it = 0; it < m_trksForThisEvent->size(); it++)
394  {
395  TrackWithHistory * trkH = (*m_trksForThisEvent)[it];
396  int genParticleID_ = trkH->genParticleID();
397  if ( genParticleID_ == -1 ) { continue; }
398  else {
399  for ( unsigned int itrlink = 0; itrlink < (*theLHCTlink).size(); itrlink++ ) {
400  if ( (*theLHCTlink)[itrlink].afterHector() == genParticleID_ ) {
401  trkH->setGenParticleID( (*theLHCTlink)[itrlink].beforeHector() );
402  continue;
403  }
404  }
405  }
406  }
407 
408  theLHCTlink = 0;
409 
410 }
#define LogDebug(id)
int i
Definition: DBlmapReader.cc:9
int getOrCreateVertex(TrackWithHistory *, int, G4SimEvent *simEvent)
bool saved() const
double totalEnergy() const
SimTrackManager(bool iCollapsePrimaryVertices=false)
void setGenParticleID(int i)
void cleanTkCaloStateInfoMap()
std::map< uint32_t, std::pair< math::XYZVectorD, math::XYZTLorentzVectorD > > mapTkCaloStateInfo
void swap(Association< C > &lhs, Association< C > &rhs)
Definition: Association.h:116
unsigned int lastHist
virtual ~SimTrackManager()
TrackContainer * m_trksForThisEvent
U second(std::pair< T, U > const &p)
void add(G4SimTrack *t)
Definition: G4SimEvent.h:35
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > XYZVectorD
spatial vector with cartesian internal representation
Definition: Vector3D.h:8
std::vector< TrackWithHistory * > TrackContainer
Definition: TrackContainer.h:8
T sqrt(T t)
Definition: SSEVec.h:18
std::vector< std::pair< int, int > > idsave
unsigned int trackID() const
void cleanTracksWithHistory()
int j
Definition: DBlmapReader.cc:9
std::pair< int, math::XYZVectorD > MapVertexPosition
this map contains association between vertex number and position
MotherParticleToVertexMap m_vertexMap
int particleID() const
void storeTracks(G4SimEvent *simEvent)
int idSavedTrack(int) const
int parentID() const
const G4VProcess * creatorProcess() const
std::map< int, MapVertexPositionVector > MotherParticleToVertexMap
const math::XYZVectorD & momentum() const
void saveTrackAndItsBranch(TrackWithHistory *)
this saves a track and all its parents looping over the non ordered vector
unsigned int lastTrack
double globalTime() const
static const GlobalPoint notFound(0, 0, 0)
const edm::LHCTransportLinkContainer * theLHCTlink
void reallyStoreTracks(G4SimEvent *simEvent)
std::vector< std::pair< int, int > > ancestorList
const math::XYZVectorD & vertexPosition() const
int genParticleID() const