BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
BesMcTruthWriter.cc
Go to the documentation of this file.
1//---------------------------------------------------------------------------//
2//// BOOST --- BESIII Object_Oriented Simulation Tool //
3////---------------------------------------------------------------------------//
4////Description:
5////Author: Dengzy
6////Created: Mar, 2004
7////Modified:
8////Comment:
9//
10#include "EmcSim/BesEmcHit.hh"
11#include "G4DigiManager.hh"
12#include "MdcSim/BesMdcHit.hh"
13#include "MucSim/BesMucHit.hh"
14#include "TofSim/BesTofHit.hh"
15#include "TruSim/BesSensitiveManager.hh"
16#include "TruSim/BesTruthTrack.hh"
17#include "TruSim/BesTruthVertex.hh"
18// #include "G4HCofThisEvent.hh"
19// #include "G4SDManager.hh"
20// #include "G4PrimaryVertex.hh"
21// #include "G4PrimaryParticle.hh"
22
23#include "CLHEP/Geometry/Point3D.h"
24#include "CLHEP/Vector/LorentzVector.h"
25#include "GaudiKernel/Bootstrap.h"
26#include "GaudiKernel/IDataProviderSvc.h"
27#include "GaudiKernel/ISvcLocator.h"
28#include "GaudiKernel/MsgStream.h"
29#include "GaudiKernel/RegistryEntry.h"
30
31#include "McTruth/DecayMode.h"
32#include "McTruth/EmcMcHit.h"
33#include "McTruth/MdcMcHit.h"
34#include "McTruth/MucMcHit.h"
35#include "McTruth/TofMcHit.h"
36
37#include "Identifier/EmcID.h"
38#include "Identifier/Identifier.h"
39#include "Identifier/MdcID.h"
40#include "Identifier/MucID.h"
41#include "Identifier/TofID.h"
42
43#include "EventModel/EventModel.h"
44#include "McTruth/McEvent.h"
45
46#include "BesMcTruthWriter.hh"
47#include "GaudiKernel/SmartDataPtr.h"
48
49BesMcTruthWriter::BesMcTruthWriter() {
50 m_DigiMan = G4DigiManager::GetDMpointer();
51 // mdcGeoPointer=BesMdcGeoParameter::GetGeo();
52}
53
54BesMcTruthWriter::~BesMcTruthWriter() {}
55
56void BesMcTruthWriter::SaveMcTruth() {
57 // interface to event data service
58 ISvcLocator* svcLocator = Gaudi::svcLocator();
59 StatusCode sc = svcLocator->service( "EventDataSvc", m_evtSvc );
60 if ( sc.isFailure() ) G4cout << "Could not accesss EventDataSvc!" << G4endl;
61
62 DataObject* aMcEvent;
63 m_evtSvc->findObject( "/Event/MC", aMcEvent );
64 if ( aMcEvent == NULL )
65 {
66 McEvent* aMcEvent = new McEvent;
67 sc = m_evtSvc->registerObject( "/Event/MC", aMcEvent );
68 if ( sc != StatusCode::SUCCESS ) G4cout << "Could not register McEvent" << G4endl;
69 }
70
71 SaveMcParticle();
72 SaveDecayMode();
73 SaveMdcTruth();
74 SaveTofTruth();
75 SaveEmcTruth();
76 SaveMucTruth();
77}
78
79void BesMcTruthWriter::SaveMcParticle() {
80 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
81 vector<BesTruthTrack*>* trackList = sensitiveManager->GetTrackList();
82 G4int nTrack = trackList->size();
83 BesTruthTrack* track;
84
85 vector<BesTruthVertex*>* vertexList = sensitiveManager->GetVertexList();
86 G4int nVertex = vertexList->size();
87 BesTruthVertex* vertex;
88
89 // arrange TruthTrack in trackList in order of trackIndex
90 for ( int i = 0; i < nTrack - 1; i++ )
91 for ( int j = i + 1; j < nTrack; j++ )
92 if ( ( *trackList )[i]->GetIndex() > ( *trackList )[j]->GetIndex() )
93 {
94 track = ( *trackList )[i];
95 ( *trackList )[i] = ( *trackList )[j];
96 ( *trackList )[j] = track;
97 }
98
99 Event::McParticleCol* particleCol = new Event::McParticleCol;
100
101 // loop over tracks
102 for ( int i = 0; i < nTrack; i++ )
103 {
104 track = ( *trackList )[i];
105
106 // find out the start point
107 bool isPrimary = false;
108 bool startPointFound = false;
109 BesTruthVertex* startPoint = nullptr;
110
111 for ( int i = 0; i < nVertex; i++ )
112 {
113 vertex = ( *vertexList )[i];
114 if ( track->GetVertex()->GetIndex() == vertex->GetIndex() )
115 {
116 if ( !vertex->GetParentTrack() ) isPrimary = true;
117 startPointFound = true;
118 startPoint = vertex;
119 break;
120 }
121 }
122
123 if ( !startPointFound )
124 std::cout << "error in finding the start point of a track" << std::endl;
125
126 bool endPointFound = false;
127 // find out the end point
128 for ( int i = 0; i < nVertex; i++ )
129 {
130 vertex = ( *vertexList )[i];
131 if ( track->GetTerminalVertex() )
132 {
133 if ( track->GetTerminalVertex()->GetIndex() == vertex->GetIndex() )
134 {
135 // create the mc particle
136 Event::McParticle* mcParticle = new Event::McParticle;
137
138 // initialization
139 HepLorentzVector initialMomentum(
140 track->GetP4().x() / 1000., track->GetP4().y() / 1000.,
141 track->GetP4().z() / 1000., track->GetP4().t() / 1000. );
142
143 HepLorentzVector initialPosition(
144 startPoint->GetPosition().x() / 10., startPoint->GetPosition().y() / 10.,
145 startPoint->GetPosition().z() / 10., startPoint->GetTime() );
146
147 mcParticle->initialize( track->GetPDGCode(), 0, initialMomentum, initialPosition );
148
149 // Set track index
150 mcParticle->setTrackIndex( track->GetIndex() );
151
152 // Adding status flag
153 if ( isPrimary ) mcParticle->addStatusFlag( Event::McParticle::PRIMARY );
154
155 if ( track->GetDaughterIndexes().size() == 0 )
156 mcParticle->addStatusFlag( Event::McParticle::LEAF );
157
158 // std::cout<<"pdg:"<<track->GetPDGCode()<<" "
159 // <<"source:"<<track->GetSource()<<" "
160 // <<std::endl;
161 if ( track->GetSource() == "FromGenerator" )
162 mcParticle->addStatusFlag( Event::McParticle::DECAYED );
163 else if ( track->GetSource() == "FromG4" )
164 mcParticle->addStatusFlag( Event::McParticle::DECAYFLT );
165 else mcParticle->addStatusFlag( Event::McParticle::ERROR );
166
167 // std::cout<<"status:"<<mcParticle->statusFlags()<<" "
168 // <<"FGN:"<<mcParticle->decayFromGenerator()<<" "
169 // <<"FG4:"<<mcParticle->decayInFlight()<<std::endl;
170 // Adding vertex index
171 mcParticle->setVertexIndex0( startPoint->GetIndex() );
172 mcParticle->setVertexIndex1( vertex->GetIndex() );
173
174 // Set the final position
175 HepLorentzVector finalPosition( vertex->GetPosition().x() / 10.,
176 vertex->GetPosition().y() / 10.,
177 vertex->GetPosition().z() / 10., vertex->GetTime() );
178 mcParticle->finalize( finalPosition );
179
180 // push back
181 particleCol->push_back( mcParticle );
182
183 endPointFound = true;
184 }
185 }
186 }
187
188 if ( !endPointFound )
189 {
190 // std::cout << "warning in finding the end point of a track" <<std::endl;
191 if ( track->GetDaughterIndexes().size() == 0 )
192 {
193 // create the mc particle
194 Event::McParticle* mcParticle = new Event::McParticle;
195 // initialization
196 HepLorentzVector initialMomentum(
197 track->GetP4().x() / 1000., track->GetP4().y() / 1000., track->GetP4().z() / 1000.,
198 track->GetP4().t() / 1000. );
199 HepLorentzVector initialPosition(
200 startPoint->GetPosition().x() / 10., startPoint->GetPosition().y() / 10.,
201 startPoint->GetPosition().z() / 10., startPoint->GetTime() );
202
203 mcParticle->initialize( track->GetPDGCode(), 0, initialMomentum, initialPosition );
204
205 // Set track index
206 mcParticle->setTrackIndex( track->GetIndex() );
207
208 // Adding status flag
209 mcParticle->addStatusFlag( Event::McParticle::LEAF );
210 if ( isPrimary ) mcParticle->addStatusFlag( Event::McParticle::PRIMARY );
211
212 // Adding vertex index
213 mcParticle->setVertexIndex0( startPoint->GetIndex() );
214 mcParticle->setVertexIndex1( -99 );
215
216 // std::cout<<"pdg:"<<track->GetPDGCode()<<" "
217 // <<"source:"<<track->GetSource()<<" "
218 // <<std::endl;
219 if ( track->GetSource() == "FromGenerator" )
220 mcParticle->addStatusFlag( Event::McParticle::DECAYED );
221 else if ( track->GetSource() == "FromG4" )
222 mcParticle->addStatusFlag( Event::McParticle::DECAYFLT );
223 else mcParticle->addStatusFlag( Event::McParticle::ERROR );
224
225 // std::cout<<"status:"<<mcParticle->statusFlags()<<" "
226 // <<"FGN:"<<mcParticle->decayFromGenerator()<<" "
227 // <<"FG4:"<<mcParticle->decayInFlight()<<std::endl;
228
229 // push back
230 particleCol->push_back( mcParticle );
231 continue;
232 }
233 }
234 }
235
236 // Get primary McParticles
237 SmartRefVector<Event::McParticle> primaryParticleCol;
238 Event::McParticleCol::iterator iter_particle;
239 for ( iter_particle = particleCol->begin(); iter_particle != particleCol->end();
240 iter_particle++ )
241 {
242
243 if ( ( *iter_particle )->primaryParticle() )
244 {
245 Event::McParticle* mcPart = (Event::McParticle*)( *iter_particle );
246 primaryParticleCol.push_back( mcPart );
247 }
248 }
249
250 if ( primaryParticleCol.empty() ) std::cout << "no primary particle found!" << std::endl;
251
252 // Add mother particle recursively
253 SmartRefVector<Event::McParticle>::iterator iter_primary;
254 for ( iter_primary = primaryParticleCol.begin(); iter_primary != primaryParticleCol.end();
255 iter_primary++ )
256 {
257 if ( !( ( *iter_primary )->leafParticle() ) ) AddMother( ( *iter_primary ), particleCol );
258 }
259
260 // register McParticle collection to TDS
261 StatusCode sc = m_evtSvc->registerObject( "/Event/MC/McParticleCol", particleCol );
262 if ( sc != StatusCode::SUCCESS )
263 G4cout << "Could not register McParticle collection" << G4endl;
264
265 // retrive McParticle from TDS
266 /*SmartDataPtr<Event::McParticleCol> mcParticleCol(m_evtSvc,"/Event/MC/McParticleCol");
267 if(!mcParticleCol)
268 G4cout<<"Could not retrieve McParticelCol"<<G4endl;
269 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
270 for (;iter_mc != mcParticleCol->end(); iter_mc++)
271 {
272 G4cout<<(*iter_mc)->getTrackIndex();
273 G4cout<<" "<<(*iter_mc)->particleProperty();
274 G4cout<<" "<<(*iter_mc)->getVertexIndex0();
275 G4cout<<" "<<(*iter_mc)->getVertexIndex1();
276 G4cout<<" "<<(*iter_mc)->initialFourMomentum().x();
277 G4cout<<" "<<(*iter_mc)->initialFourMomentum().y();
278 G4cout<<" "<<(*iter_mc)->initialFourMomentum().z();
279 G4cout<<" "<<(*iter_mc)->initialFourMomentum().t();
280 G4cout<<" "<<(*iter_mc)->initialPosition().x();
281 G4cout<<" "<<(*iter_mc)->initialPosition().y();
282 G4cout<<" "<<(*iter_mc)->initialPosition().z();
283 G4cout<<G4endl;
284 }
285 G4cout<<"end of retrieve McParticle from TDS"<<G4endl;
286 */
287}
288
289void BesMcTruthWriter::AddMother( Event::McParticle* currentParticle,
290 Event::McParticleCol* particleCol ) {
291 if ( currentParticle->leafParticle() ) { return; }
292 Event::McParticleCol::iterator iter;
293 bool found = false;
294 for ( iter = particleCol->begin(); iter != particleCol->end(); iter++ )
295 {
296 if ( currentParticle->vertexIndex1() == ( *iter )->vertexIndex0() )
297 {
298 found = true;
299 ( *iter )->setMother( currentParticle );
300 currentParticle->addDaughter( *iter );
301 AddMother( ( *iter ), particleCol );
302 }
303 }
304 if ( !found ) std::cout << "AddMother: inconsistency was found!" << std::endl;
305}
306
307void BesMcTruthWriter::SaveDecayMode() {
308 SmartDataPtr<DecayMode> decayMode( m_evtSvc, "/Event/MC/DecayMode" );
309 if ( !decayMode )
310 {
311 // G4cout<<"Could not retrieve DecayMode"<<G4endl;
312 DecayMode* aDecayMode = new DecayMode;
313 // register Decay Mode to TDS
314 int dm[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
315 aDecayMode->putData( dm, 10 );
316 StatusCode scDM = m_evtSvc->registerObject( "/Event/MC/DecayMode", aDecayMode );
317 if ( scDM != StatusCode::SUCCESS ) G4cout << "Could not register DecayMode" << G4endl;
318 }
319}
320
321void BesMcTruthWriter::SaveMdcTruth() {
322 // Mdc McHit collection defined in BOSS
323 Event::MdcMcHitCol* aMdcMcHitCol = new Event::MdcMcHitCol;
324
325 G4int HCID = -1;
326 HCID = m_DigiMan->GetHitsCollectionID( "BesMdcTruthCollection" );
327 if ( HCID > 0 )
328 {
329 BesMdcHitsCollection* HC = 0;
330 HC = (BesMdcHitsCollection*)( m_DigiMan->GetHitsCollection( HCID ) );
331 G4int n_hit = HC->entries();
332 if ( n_hit > 0 )
333 {
334 // arrange hits in hits collection in order of trackIndex
335 BesMdcHit* hit;
336 vector<BesMdcHit*>* vecHC = HC->GetVector();
337 for ( int i = 0; i < n_hit - 1; i++ )
338 for ( int j = i + 1; j < n_hit; j++ )
339 if ( ( *vecHC )[i]->GetTrackID() > ( *vecHC )[j]->GetTrackID() )
340 {
341 hit = ( *vecHC )[i];
342 ( *vecHC )[i] = ( *vecHC )[j];
343 ( *vecHC )[j] = hit;
344 }
345
346 // push back MdcMcHits to MdcMcHitCol in TDS
347 for ( G4int i = 0; i < n_hit; i++ )
348 {
349 hit = ( *HC )[i];
350 const Identifier ident = MdcID::wire_id( hit->GetLayerNo(), hit->GetCellNo() );
351 Event::MdcMcHit* mdcMcHit = new Event::MdcMcHit(
352 ident, hit->GetTrackID(), hit->GetPos().x(), hit->GetPos().y(), hit->GetPos().z(),
353 hit->GetDriftD(), hit->GetEdep(), hit->GetPosFlag() );
354 aMdcMcHitCol->push_back( mdcMcHit );
355 }
356 }
357 }
358
359 // register MDC McTruth collection to TDS
360 StatusCode scMdc = m_evtSvc->registerObject( "/Event/MC/MdcMcHitCol", aMdcMcHitCol );
361 if ( scMdc != StatusCode::SUCCESS )
362 G4cout << "Could not register MDC McTruth collection" << G4endl;
363
364 // retrieve MDC McTruth from TDS
365 /*SmartDataPtr<Event::MdcMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/MdcMcHitCol");
366 if(!aMcHitCol)
367 G4cout<<"Could not retrieve MDC McTruth collection"<<G4endl;
368
369 Event::MdcMcHitCol::iterator iMcHitCol;
370 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
371 {
372 const Identifier ident = (*iMcHitCol)->identify();
373 G4cout<<(*iMcHitCol)->getTrackIndex();
374 G4cout<<" "<<MdcID::layer(ident);
375 G4cout<<" "<<MdcID::wire(ident);
376 G4cout<<" "<<(*iMcHitCol)->getDepositEnergy();
377 G4cout<<" "<<(*iMcHitCol)->getDriftDistance();
378 G4cout<<" "<<(*iMcHitCol)->getPositionX();
379 G4cout<<" "<<(*iMcHitCol)->getPositionY();
380 G4cout<<" "<<(*iMcHitCol)->getPositionZ();
381 G4cout<<G4endl;
382 }
383 G4cout<<"end of retrieve MDC McTruth collection"<<G4endl;
384 */
385}
386
387void BesMcTruthWriter::SaveTofTruth() {
388 // Tof McHit collection defined in BOSS
389 Event::TofMcHitCol* aTofMcHitCol = new Event::TofMcHitCol;
390
391 G4int HCID = -1;
392 HCID = m_DigiMan->GetHitsCollectionID( "BesTofHitsList" );
393 if ( HCID > 0 )
394 {
395 BesTofHitsCollection* HC = 0;
396 HC = (BesTofHitsCollection*)( m_DigiMan->GetHitsCollection( HCID ) );
397 G4int n_hit = HC->entries();
398 if ( n_hit > 0 )
399 {
400 // arrange hits in hits collection in order of trackIndex
401 BesTofHit* hit;
402 vector<BesTofHit*>* vecHC = HC->GetVector();
403 for ( int i = 0; i < n_hit - 1; i++ )
404 {
405 for ( int j = i + 1; j < n_hit; j++ )
406 {
407 if ( ( *vecHC )[i]->GetTrackIndex() > ( *vecHC )[j]->GetTrackIndex() )
408 {
409 hit = ( *vecHC )[i];
410 ( *vecHC )[i] = ( *vecHC )[j];
411 ( *vecHC )[j] = hit;
412 }
413 }
414 }
415
416 // push back TofMcHit collection to TDS
417 for ( G4int i = 0; i < n_hit; i++ )
418 {
419 hit = ( *HC )[i];
420 Identifier ident;
421 unsigned int barrel_ec = hit->GetPartId();
422 // for Scintillator
423 if ( TofID::is_scin( barrel_ec ) )
424 {
425 unsigned int scinNum = hit->GetScinNb();
426 unsigned int layer = 0;
427 if ( TofID::is_barrel( barrel_ec ) && scinNum > TofID::getPHI_BARREL_MAX() )
428 {
429 layer = 1;
430 scinNum = scinNum - TofID::getPHI_BARREL_MAX() - 1;
431 }
432 ident = TofID::cell_id( barrel_ec, layer, scinNum, 0 );
433 }
434 // for ETF(MRPC)
435 else
436 { // for ETF(MRPC)
437 if ( barrel_ec == 3 || barrel_ec == 4 )
438 {
439 unsigned int endcap = 0;
440 unsigned int module = hit->GetModule();
441 unsigned int strip = hit->GetStrip();
442 if ( barrel_ec == 4 )
443 { // west end cap
444 endcap = 1;
445 }
446 ident = TofID::cell_id( 3, endcap, module, strip, 0 );
447 }
448 }
449 if ( barrel_ec >= 0 && barrel_ec <= 4 )
450 {
451 Event::TofMcHit* tofMcHit = new Event::TofMcHit(
452 ident, hit->GetTrackIndex(), hit->GetPos().x(), hit->GetPos().y(),
453 hit->GetPos().z(), hit->GetMomentum().x(), hit->GetMomentum().y(),
454 hit->GetMomentum().z(), hit->GetTrackL(), hit->GetTime() );
455 aTofMcHitCol->push_back( tofMcHit );
456 }
457 }
458 }
459 }
460
461 // register TOF McTruth collection to TDS
462 StatusCode scTof = m_evtSvc->registerObject( "/Event/MC/TofMcHitCol", aTofMcHitCol );
463 if ( scTof != StatusCode::SUCCESS )
464 G4cout << "Could not register TOF McTruth collection" << G4endl;
465
466 // retrieve TOF McTruth from TDS
467 /*SmartDataPtr<Event::TofMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/TofMcHitCol");
468 if(!aMcHitCol)
469 G4cout<<"Could not retrieve TOF McTruth collection"<<G4endl;
470
471 Event::TofMcHitCol::iterator iMcHitCol;
472 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
473 {
474 const Identifier ident = (*iMcHitCol)->identify();
475 G4cout<<(*iMcHitCol)->getTrackIndex();
476 G4cout<<" "<<TofID::barrel_ec(ident);;
477 G4cout<<" "<<TofID::layer(ident);
478 G4cout<<" "<<TofID::phi_module(ident);
479 G4cout<<" "<<(*iMcHitCol)->getPositionX();
480 G4cout<<" "<<(*iMcHitCol)->getPositionY();
481 G4cout<<" "<<(*iMcHitCol)->getPositionZ();
482 G4cout<<" "<<(*iMcHitCol)->getPx();
483 G4cout<<" "<<(*iMcHitCol)->getPy();
484 G4cout<<" "<<(*iMcHitCol)->getPz();
485 G4cout<<" "<<(*iMcHitCol)->getTrackLength();
486 G4cout<<" "<<(*iMcHitCol)->getFlightTime();
487 G4cout<<G4endl;
488 }
489 G4cout<<"end of retrieve TOF McTruth"<<G4endl;
490 */
491}
492
493void BesMcTruthWriter::SaveEmcTruth() {
494 // Emc McHit collection defined in BOSS
495 Event::EmcMcHitCol* aEmcMcHitCol = new Event::EmcMcHitCol;
496
497 G4int fullMc = 1;
498 if ( fullMc == 1 )
499 { // full Mc Truth
500 G4int HCID = -1;
501 HCID = m_DigiMan->GetHitsCollectionID( "BesEmcTruthHitsList" );
502 if ( HCID > 0 )
503 {
504 BesEmcTruthHitsCollection* HC = 0;
505 HC = (BesEmcTruthHitsCollection*)( m_DigiMan->GetHitsCollection( HCID ) );
506 G4int n_hit = HC->entries();
507 if ( n_hit > 0 )
508 {
509 // arrange hits in hits collection in order of trackIndex
510 BesEmcTruthHit* hit;
511 vector<BesEmcTruthHit*>* vecHC = HC->GetVector();
512 for ( int i = 0; i < n_hit - 1; i++ )
513 {
514 for ( int j = i + 1; j < n_hit; j++ )
515 {
516 if ( ( *vecHC )[i]->GetTrackIndex() > ( *vecHC )[j]->GetTrackIndex() )
517 {
518 hit = ( *vecHC )[i];
519 ( *vecHC )[i] = ( *vecHC )[j];
520 ( *vecHC )[j] = hit;
521 }
522 }
523 }
524
525 for ( G4int i = 0; i < n_hit; i++ )
526 {
527 BesEmcTruthHit* hit;
528 hit = ( *HC )[i];
529
530 std::map<Identifier, double> hitMap;
531 std::map<Identifier, double>::const_iterator iHitMap;
532 hitMap.clear();
533 for ( iHitMap = hit->Begin(); iHitMap != hit->End(); iHitMap++ )
534 { hitMap[iHitMap->first] = iHitMap->second; }
535
536 Event::EmcMcHit* emcMcHit = new Event::EmcMcHit(
537 hit->GetIdentify(), hit->GetTrackIndex(), hit->GetPosition().x(),
538 hit->GetPosition().y(), hit->GetPosition().z(), hit->GetMomentum().x(),
539 hit->GetMomentum().y(), hit->GetMomentum().z(), hit->GetEDep() );
540
541 emcMcHit->setHitEmc( hit->GetHitEmc() );
542 emcMcHit->setPDGCode( hit->GetPDGCode() );
543 emcMcHit->setPDGCharge( hit->GetPDGCharge() );
544 emcMcHit->setTime( hit->GetTime() );
545 emcMcHit->setHitMap( hitMap );
546
547 aEmcMcHitCol->push_back( emcMcHit );
548 }
549 }
550 }
551 }
552 else
553 { // simple Mc Truth
554 G4int HCID = -1;
555 HCID = m_DigiMan->GetHitsCollectionID( "BesEmcHitsList" );
556 if ( HCID > 0 )
557 {
558 BesEmcHitsCollection* HC = 0;
559 HC = (BesEmcHitsCollection*)( m_DigiMan->GetHitsCollection( HCID ) );
560 G4int n_hit = HC->entries();
561 if ( n_hit > 0 )
562 {
563 // arrange hits in hits collection in order of trackIndex
564 BesEmcHit* hit;
565 vector<BesEmcHit*>* vecHC = HC->GetVector();
566 for ( int i = 0; i < n_hit - 1; i++ )
567 for ( int j = i + 1; j < n_hit; j++ )
568 if ( ( *vecHC )[i]->GetTrackIndex() > ( *vecHC )[j]->GetTrackIndex() )
569 {
570 hit = ( *vecHC )[i];
571 ( *vecHC )[i] = ( *vecHC )[j];
572 ( *vecHC )[j] = hit;
573 }
574
575 // Emc McHit collection defined in BOSS
576 Event::EmcMcHitCol* aEmcMcHitCol = new Event::EmcMcHitCol;
577
578 for ( G4int i = 0; i < n_hit; i++ )
579 {
580 hit = ( *HC )[i];
581 Identifier ident = EmcID::crystal_id( hit->GetPartId(), hit->GetNumThetaCrystal(),
582 hit->GetNumPhiCrystal() );
583
584 std::map<Identifier, double> hitMap;
585 Event::EmcMcHit* emcMcHit = new Event::EmcMcHit(
586 ident, hit->GetTrackIndex(), hit->GetPosCrystal().x(), hit->GetPosCrystal().y(),
587 hit->GetPosCrystal().z(), hit->GetMomentum().x(), hit->GetMomentum().y(),
588 hit->GetMomentum().z(), hit->GetEdepCrystal() );
589
590 aEmcMcHitCol->push_back( emcMcHit );
591 }
592 }
593 }
594 }
595
596 // register EMC McTruth collection to TDS
597 StatusCode scEmc = m_evtSvc->registerObject( "/Event/MC/EmcMcHitCol", aEmcMcHitCol );
598 if ( scEmc != StatusCode::SUCCESS )
599 G4cout << "Could not register EMC McTruth collection" << G4endl;
600
601 // retrieve EMC McTruth from TDS
602 /*SmartDataPtr<Event::EmcMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/EmcMcHitCol");
603 if(!aMcHitCol)
604 G4cout<<"Could not retrieve EMC McTruth collection"<<G4endl;
605
606 Event::EmcMcHitCol::iterator iMcHitCol;
607 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
608 {
609 const Identifier ident = (*iMcHitCol)->identify();
610 G4cout<<(*iMcHitCol)->getTrackIndex();
611 G4cout<<" "<<EmcID::barrel_ec(ident);
612 G4cout<<" "<<EmcID::theta_module(ident);
613 G4cout<<" "<<EmcID::phi_module(ident);
614 G4cout<<" "<<(*iMcHitCol)->getPositionX();
615 G4cout<<" "<<(*iMcHitCol)->getPositionY();
616 G4cout<<" "<<(*iMcHitCol)->getPositionZ();
617 G4cout<<" "<<(*iMcHitCol)->getPx();
618 G4cout<<" "<<(*iMcHitCol)->getPy();
619 G4cout<<" "<<(*iMcHitCol)->getPz();
620 G4cout<<" "<<(*iMcHitCol)->getDepositEnergy();
621 G4cout<<G4endl;
622 }
623 G4cout<<"end of retrieve EMC McTruth"<<G4endl;
624 */
625}
626
627void BesMcTruthWriter::SaveMucTruth() {
628 // Muc McHit collection defined in BOSS
629 Event::MucMcHitCol* aMucMcHitCol = new Event::MucMcHitCol;
630
631 G4int HCID = -1;
632 HCID = m_DigiMan->GetHitsCollectionID( "BesMucHitsList" );
633 if ( HCID > 0 )
634 {
635 BesMucHitsCollection* HC = 0;
636 HC = (BesMucHitsCollection*)( m_DigiMan->GetHitsCollection( HCID ) );
637 G4int n_hit = HC->entries();
638 if ( n_hit > 0 )
639 {
640 // arrange hits in hits collection in order of trackIndex
641 BesMucHit* hit;
642 vector<BesMucHit*>* vecHC = HC->GetVector();
643 for ( int i = 0; i < n_hit - 1; i++ )
644 for ( int j = i + 1; j < n_hit; j++ )
645 if ( ( *vecHC )[i]->GetTrackIndex() > ( *vecHC )[j]->GetTrackIndex() )
646 {
647 hit = ( *vecHC )[i];
648 ( *vecHC )[i] = ( *vecHC )[j];
649 ( *vecHC )[j] = hit;
650 }
651
652 for ( G4int i = 0; i < n_hit; i++ )
653 {
654 hit = ( *HC )[i];
655 Identifier ident =
656 MucID::channel_id( hit->GetPart(), hit->GetSeg(), hit->GetGap(), hit->GetStrip() );
657 Event::MucMcHit* mucMcHit =
658 new Event::MucMcHit( ident, hit->GetTrackIndex(), hit->GetPos().x(),
659 hit->GetPos().y(), hit->GetPos().z(), hit->GetMomentum().x(),
660 hit->GetMomentum().y(), hit->GetMomentum().z() );
661 aMucMcHitCol->push_back( mucMcHit );
662 }
663 }
664 }
665
666 // register MUC McTruth collection to TDS
667 StatusCode scMuc = m_evtSvc->registerObject( "/Event/MC/MucMcHitCol", aMucMcHitCol );
668 if ( scMuc != StatusCode::SUCCESS )
669 G4cout << "Could not register MUC McTruth collection" << G4endl;
670
671 // retrieve MUC McTruth from TDS
672 /*SmartDataPtr<Event::MucMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/MucMcHitCol");
673 if(!aMcHitCol)
674 G4cout<<"Could not retrieve MUC McTruth collection"<<G4endl;
675
676 Event::MucMcHitCol::iterator iMcHitCol;
677 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
678 {
679 const Identifier ident = (*iMcHitCol)->identify();
680 G4cout<<(*iMcHitCol)->getTrackIndex();
681 G4cout<<" "<<MucID::part(ident);
682 G4cout<<" "<<MucID::seg(ident);
683 G4cout<<" "<<MucID::gap(ident);
684 G4cout<<" "<<MucID::strip(ident);
685 G4cout<<" "<<(*iMcHitCol)->getPositionX();
686 G4cout<<" "<<(*iMcHitCol)->getPositionY();
687 G4cout<<" "<<(*iMcHitCol)->getPositionZ();
688 G4cout<<" "<<(*iMcHitCol)->getPx();
689 G4cout<<" "<<(*iMcHitCol)->getPy();
690 G4cout<<" "<<(*iMcHitCol)->getPz();
691 G4cout<<G4endl;
692 }
693 G4cout<<"end of retrieve MUC McTruth"<<G4endl;
694 */
695}
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
BesEmcTruthHitsCollection
G4THitsCollection< BesEmcTruthHit > BesEmcTruthHitsCollection
Definition InstallArea/x86_64-el9-gcc13-dbg/include/EmcSim/BesEmcHit.hh:177
BesEmcHitsCollection
G4THitsCollection< BesEmcHit > BesEmcHitsCollection
Definition InstallArea/x86_64-el9-gcc13-dbg/include/EmcSim/BesEmcHit.hh:86
BesMdcHitsCollection
G4THitsCollection< BesMdcHit > BesMdcHitsCollection
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcSim/BesMdcHit.hh:76
BesMucHitsCollection
G4THitsCollection< BesMucHit > BesMucHitsCollection
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MucSim/BesMucHit.hh:100
BesTofHitsCollection
G4THitsCollection< BesTofHit > BesTofHitsCollection
Definition InstallArea/x86_64-el9-gcc13-dbg/include/TofSim/BesTofHit.hh:113
BesEmcTruthHit::End
std::map< Identifier, G4double >::const_iterator End() const
Definition BesEmcHit.cc:163
BesEmcTruthHit::Begin
std::map< Identifier, G4double >::const_iterator Begin() const
Definition BesEmcHit.cc:159
BesMcTruthWriter::AddMother
void AddMother(Event::McParticle *, Event::McParticleCol *)
Definition BesMcTruthWriter.cc:289
DecayMode::putData
void putData(int *data, unsigned int size)
Definition Event/McTruth/src/DecayMode.cxx:9
EmcID::crystal_id
static Identifier crystal_id(const unsigned int barrel_ec, const unsigned int theta_module, const unsigned int phi_module)
For a single crystal.
Definition EmcID.cxx:63
Event::EmcMcHit::setHitMap
void setHitMap(std::map< Identifier, double > &hitMap)
Definition Event/McTruth/include/McTruth/EmcMcHit.h:88
Event::McParticle::setVertexIndex0
void setVertexIndex0(int index0)
methods for setting and getting vertex indexes
Definition Event/McTruth/include/McTruth/McParticle.h:104
Event::McParticle::ERROR
@ ERROR
this particle is a leaf in the particle tree
Definition Event/McTruth/include/McTruth/McParticle.h:38
Event::McParticle::initialize
void initialize(StdHepId id, unsigned int statusBits, const HepLorentzVector &initialMomentum, const HepLorentzVector &initialPosition, const std::string process="")
Set the initial attributes of the McParticle.
Definition McParticle.cxx:33
Event::McParticle::leafParticle
bool leafParticle() const
Retrieve whether this is a leaf particle.
Definition McParticle.cxx:13
Event::McParticle::addStatusFlag
void addStatusFlag(unsigned int flag)
add a new flag to the status flags
Definition Event/McTruth/include/McTruth/McParticle.h:88
Event::McParticle::finalize
void finalize(const HepLorentzVector &finalPosition)
Set the final attributes of the McParticle.
Definition McParticle.cxx:64
Event::McParticle::addDaughter
void addDaughter(const SmartRef< McParticle > d)
add a daugther particle to this particle
Definition Event/McTruth/include/McTruth/McParticle.h:132
MdcID::wire_id
static Identifier wire_id(int wireType, int layer, int wire)
For a single wire.
Definition MdcID.cxx:69
MucID::channel_id
static Identifier channel_id(int barrel_ec, int segment, int layer, int channel)
For a single crystal.
Definition MucID.cxx:114
TofID::cell_id
static Identifier cell_id(int barrel_ec, int layer, int phi_module, int end)
For a single crystal.
Definition TofID.cxx:126
TofID::is_scin
static bool is_scin(const Identifier &id)
Definition TofID.cxx:88
TofID::getPHI_BARREL_MAX
static value_type getPHI_BARREL_MAX()
Definition TofID.cxx:167
TofID::is_barrel
static bool is_barrel(const Identifier &id)
Test for barrel.
Definition TofID.cxx:40
Event::MucMcHitCol
ObjectVector< MucMcHit > MucMcHitCol
Definition Event/McTruth/include/McTruth/MucMcHit.h:82
Event::McParticleCol
ObjectList< McParticle > McParticleCol
Definition Event/McTruth/include/McTruth/McParticle.h:188
Event::EmcMcHitCol
ObjectVector< EmcMcHit > EmcMcHitCol
Definition Event/McTruth/include/McTruth/EmcMcHit.h:126
Event::TofMcHitCol
ObjectVector< TofMcHit > TofMcHitCol
Definition Event/McTruth/include/McTruth/TofMcHit.h:93
Event::MdcMcHitCol
ObjectVector< MdcMcHit > MdcMcHitCol
Definition Event/McTruth/include/McTruth/MdcMcHit.h:138