85{
87 vector<BesTruthTrack*>* trackList = sensitiveManager->
GetTrackList();
88 G4int nTrack = trackList->size();
90
91 vector<BesTruthVertex*>* vertexList = sensitiveManager->
GetVertexList();
92 G4int nVertex = vertexList->size();
94
95
96 for(int i=0;i<nTrack-1;i++)
97 for(int j=i+1;j<nTrack;j++)
98 if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
99 {
100 track=(*trackList)[i];
101 (*trackList)[i]=(*trackList)[j];
102 (*trackList)[j]=track;
103 }
104
106
107
108 for(int i=0;i<nTrack;i++)
109 {
110 track = (*trackList) [i];
111
112
113 bool isPrimary = false;
114 bool startPointFound = false;
116
117 for (int i=0;i<nVertex;i++)
118 {
119 vertex = (*vertexList) [i];
121 {
123 startPointFound = true;
124 startPoint = vertex;
125 break;
126 }
127 }
128
129 if (!startPointFound)
130 std::cout << "error in finding the start point of a track" <<std::endl;
131
132
133 bool endPointFound = false;
134
135 for (int i=0;i<nVertex;i++)
136 {
137 vertex = (*vertexList) [i];
139 {
141 {
142
144
145
146 HepLorentzVector initialMomentum(track->
GetP4().x()/1000., track->
GetP4().y()/1000., track->
GetP4().z()/1000., track->
GetP4().t()/1000.);
147
149
151
152
154
155
157
160
161
162
163
168 else
170
171
172
173
174
177
178
180 mcParticle->
finalize(finalPosition);
181
182
183 particleCol->push_back(mcParticle);
184
185 endPointFound = true;
186 }
187 }
188 }
189
190 if (!endPointFound)
191 {
192
194 {
195
197
198 HepLorentzVector initialMomentum( track->
GetP4().x()/1000., track->
GetP4().y()/1000., track->
GetP4().z()/1000., track->
GetP4().t()/1000.);
200
202
203
205
206
209
210
213
214
215
216
221 else
223
224
225
226
227
228
229 particleCol->push_back(mcParticle);
230 continue;
231 }
232 }
233 }
234
235
236 SmartRefVector<Event::McParticle> primaryParticleCol;
237 Event::McParticleCol::iterator iter_particle;
238 for ( iter_particle = particleCol->begin();
239 iter_particle != particleCol->end(); iter_particle++) {
240
241 if ((*iter_particle)->primaryParticle()) {
243 primaryParticleCol.push_back(mcPart);
244 }
245 }
246
247 if (primaryParticleCol.empty())
248 std::cout << "no primary particle found!" << std::endl;
249
250
251 SmartRefVector<Event::McParticle>::iterator iter_primary;
252 for ( iter_primary = primaryParticleCol.begin();
253 iter_primary != primaryParticleCol.end(); iter_primary++) {
254 if ( !((*iter_primary)->leafParticle()) )
256 }
257
258
259 StatusCode sc = m_evtSvc->registerObject("/Event/MC/McParticleCol", particleCol);
260 if(sc!=StatusCode::SUCCESS)
261 G4cout<< "Could not register McParticle collection" <<G4endl;
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286}
std::vector< BesTruthVertex * > * GetVertexList()
std::vector< BesTruthTrack * > * GetTrackList()
static BesSensitiveManager * GetSensitiveManager()
HepLorentzVector GetP4() const
BesTruthVertex * GetTerminalVertex() const
BesTruthVertex * GetVertex() const
vector< int > GetDaughterIndexes() const
BesTruthTrack * GetParentTrack() const
G4ThreeVector GetPosition() const
void setVertexIndex0(int index0)
methods for setting and getting vertex indexes
@ PRIMARY
Decayed in flight.
@ ERROR
this particle is a leaf in the particle tree
@ DECAYFLT
Decayed by generator.
void initialize(StdHepId id, unsigned int statusBits, const HepLorentzVector &initialMomentum, const HepLorentzVector &initialPosition, const std::string process="")
Set the initial attributes of the McParticle.
void addStatusFlag(unsigned int flag)
add a new flag to the status flags
void finalize(const HepLorentzVector &finalPosition)
Set the final attributes of the McParticle.
void setVertexIndex1(int index1)
void setTrackIndex(int trackIndex)
ObjectList< McParticle > McParticleCol