00001
00009
00010 #include "G4DigiManager.hh"
00011 #include "BesMdcHit.hh"
00012 #include "BesTofHit.hh"
00013 #include "BesEmcHit.hh"
00014 #include "BesMucHit.hh"
00015 #include "BesTruthTrack.hh"
00016 #include "BesTruthVertex.hh"
00017 #include "BesSensitiveManager.hh"
00018
00019
00020
00021
00022
00023 #include "GaudiKernel/IDataProviderSvc.h"
00024 #include "GaudiKernel/ISvcLocator.h"
00025 #include "GaudiKernel/Bootstrap.h"
00026 #include "GaudiKernel/RegistryEntry.h"
00027 #include "GaudiKernel/MsgStream.h"
00028 #include "CLHEP/Vector/LorentzVector.h"
00029 #include "CLHEP/Geometry/Point3D.h"
00030
00031 #include "McTruth/DecayMode.h"
00032 #include "McTruth/MdcMcHit.h"
00033 #include "McTruth/TofMcHit.h"
00034 #include "McTruth/EmcMcHit.h"
00035 #include "McTruth/MucMcHit.h"
00036
00037 #include "Identifier/Identifier.h"
00038 #include "Identifier/MdcID.h"
00039 #include "Identifier/TofID.h"
00040 #include "Identifier/EmcID.h"
00041 #include "Identifier/MucID.h"
00042
00043 #include "McTruth/McEvent.h"
00044 #include "EventModel/EventModel.h"
00045
00046 #include "GaudiKernel/SmartDataPtr.h"
00047 #include "BesMcTruthWriter.hh"
00048
00049 BesMcTruthWriter::BesMcTruthWriter()
00050 {
00051 m_DigiMan = G4DigiManager::GetDMpointer();
00052
00053 }
00054
00055 BesMcTruthWriter::~BesMcTruthWriter()
00056 {
00057 }
00058
00059 void BesMcTruthWriter::SaveMcTruth()
00060 {
00061
00062 ISvcLocator* svcLocator = Gaudi::svcLocator();
00063 StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
00064 if (sc.isFailure())
00065 G4cout<<"Could not accesss EventDataSvc!"<<G4endl;
00066
00067 DataObject *aMcEvent;
00068 m_evtSvc->findObject("/Event/MC",aMcEvent);
00069 if(aMcEvent==NULL) {
00070 McEvent* aMcEvent = new McEvent;
00071 sc = m_evtSvc->registerObject("/Event/MC",aMcEvent);
00072 if(sc!=StatusCode::SUCCESS)
00073 G4cout<< "Could not register McEvent" <<G4endl;
00074 }
00075
00076 SaveMcParticle();
00077 SaveDecayMode();
00078 SaveMdcTruth();
00079 SaveTofTruth();
00080 SaveEmcTruth();
00081 SaveMucTruth();
00082 }
00083
00084 void BesMcTruthWriter::SaveMcParticle()
00085 {
00086 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
00087 vector<BesTruthTrack*>* trackList = sensitiveManager->GetTrackList();
00088 G4int nTrack = trackList->size();
00089 BesTruthTrack* track;
00090
00091 vector<BesTruthVertex*>* vertexList = sensitiveManager->GetVertexList();
00092 G4int nVertex = vertexList->size();
00093 BesTruthVertex* vertex;
00094
00095
00096 for(int i=0;i<nTrack-1;i++)
00097 for(int j=i+1;j<nTrack;j++)
00098 if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
00099 {
00100 track=(*trackList)[i];
00101 (*trackList)[i]=(*trackList)[j];
00102 (*trackList)[j]=track;
00103 }
00104
00105 Event::McParticleCol *particleCol = new Event::McParticleCol;
00106
00107
00108 for(int i=0;i<nTrack;i++)
00109 {
00110 track = (*trackList) [i];
00111
00112
00113 bool isPrimary = false;
00114 bool startPointFound = false;
00115 BesTruthVertex* startPoint;
00116
00117 for (int i=0;i<nVertex;i++)
00118 {
00119 vertex = (*vertexList) [i];
00120 if ( track->GetVertex()->GetIndex() == vertex->GetIndex() )
00121 {
00122 if (!vertex->GetParentTrack()) isPrimary = true;
00123 startPointFound = true;
00124 startPoint = vertex;
00125 break;
00126 }
00127 }
00128
00129 if (!startPointFound)
00130 std::cout << "error in finding the start point of a track" <<std::endl;
00131
00132
00133 bool endPointFound = false;
00134
00135 for (int i=0;i<nVertex;i++)
00136 {
00137 vertex = (*vertexList) [i];
00138 if( track->GetTerminalVertex() )
00139 {
00140 if (track->GetTerminalVertex()->GetIndex() == vertex->GetIndex() )
00141 {
00142
00143 Event::McParticle* mcParticle = new Event::McParticle;
00144
00145
00146 HepLorentzVector initialMomentum(track->GetP4().x()/1000., track->GetP4().y()/1000., track->GetP4().z()/1000., track->GetP4().t()/1000.);
00147
00148 HepLorentzVector initialPosition(startPoint->GetPosition().x()/10., startPoint->GetPosition().y()/10., startPoint->GetPosition().z()/10.,startPoint->GetTime());
00149
00150 mcParticle->initialize(track->GetPDGCode(), 0, initialMomentum, initialPosition);
00151
00152
00153 mcParticle->setTrackIndex( track->GetIndex() );
00154
00155
00156 if (isPrimary ) mcParticle->addStatusFlag(Event::McParticle::PRIMARY);
00157
00158 if (track->GetDaughterIndexes().size()==0)
00159 mcParticle->addStatusFlag(Event::McParticle::LEAF);
00160
00161
00162
00163
00164 if (track->GetSource()=="FromGenerator")
00165 mcParticle->addStatusFlag(Event::McParticle::DECAYED);
00166 else if(track->GetSource()=="FromG4")
00167 mcParticle->addStatusFlag(Event::McParticle::DECAYFLT);
00168 else
00169 mcParticle->addStatusFlag(Event::McParticle::ERROR);
00170
00171
00172
00173
00174
00175 mcParticle->setVertexIndex0(startPoint->GetIndex());
00176 mcParticle->setVertexIndex1(vertex->GetIndex());
00177
00178
00179 HepLorentzVector finalPosition(vertex->GetPosition().x()/10., vertex->GetPosition().y()/10., vertex->GetPosition().z()/10., vertex->GetTime());
00180 mcParticle->finalize(finalPosition);
00181
00182
00183 particleCol->push_back(mcParticle);
00184
00185 endPointFound = true;
00186 }
00187 }
00188 }
00189
00190 if (!endPointFound)
00191 {
00192
00193 if (track->GetDaughterIndexes().size()==0)
00194 {
00195
00196 Event::McParticle* mcParticle = new Event::McParticle;
00197
00198 HepLorentzVector initialMomentum( track->GetP4().x()/1000., track->GetP4().y()/1000., track->GetP4().z()/1000., track->GetP4().t()/1000.);
00199 HepLorentzVector initialPosition(startPoint->GetPosition().x()/10., startPoint->GetPosition().y()/10., startPoint->GetPosition().z()/10., startPoint->GetTime());
00200
00201 mcParticle->initialize(track->GetPDGCode(), 0, initialMomentum, initialPosition);
00202
00203
00204 mcParticle->setTrackIndex( track->GetIndex() );
00205
00206
00207 mcParticle->addStatusFlag(Event::McParticle::LEAF);
00208 if (isPrimary ) mcParticle->addStatusFlag(Event::McParticle::PRIMARY);
00209
00210
00211 mcParticle->setVertexIndex0( startPoint->GetIndex() );
00212 mcParticle->setVertexIndex1( -99 );
00213
00214
00215
00216
00217 if (track->GetSource()=="FromGenerator")
00218 mcParticle->addStatusFlag(Event::McParticle::DECAYED);
00219 else if(track->GetSource()=="FromG4")
00220 mcParticle->addStatusFlag(Event::McParticle::DECAYFLT);
00221 else
00222 mcParticle->addStatusFlag(Event::McParticle::ERROR);
00223
00224
00225
00226
00227
00228
00229 particleCol->push_back(mcParticle);
00230 continue;
00231 }
00232 }
00233 }
00234
00235
00236 SmartRefVector<Event::McParticle> primaryParticleCol;
00237 Event::McParticleCol::iterator iter_particle;
00238 for ( iter_particle = particleCol->begin();
00239 iter_particle != particleCol->end(); iter_particle++) {
00240
00241 if ((*iter_particle)->primaryParticle()) {
00242 Event::McParticle* mcPart = (Event::McParticle*)(*iter_particle);
00243 primaryParticleCol.push_back(mcPart);
00244 }
00245 }
00246
00247 if (primaryParticleCol.empty())
00248 std::cout << "no primary particle found!" << std::endl;
00249
00250
00251 SmartRefVector<Event::McParticle>::iterator iter_primary;
00252 for ( iter_primary = primaryParticleCol.begin();
00253 iter_primary != primaryParticleCol.end(); iter_primary++) {
00254 if ( !((*iter_primary)->leafParticle()) )
00255 AddMother((*iter_primary), particleCol);
00256 }
00257
00258
00259 StatusCode sc = m_evtSvc->registerObject("/Event/MC/McParticleCol", particleCol);
00260 if(sc!=StatusCode::SUCCESS)
00261 G4cout<< "Could not register McParticle collection" <<G4endl;
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286 }
00287
00288 void BesMcTruthWriter::AddMother(Event::McParticle* currentParticle, Event::McParticleCol *particleCol)
00289 {
00290 if (currentParticle->leafParticle()) {
00291 return;
00292 }
00293 Event::McParticleCol::iterator iter;
00294 bool found = false;
00295 for ( iter = particleCol->begin(); iter != particleCol->end(); iter++) {
00296 if (currentParticle->vertexIndex1() == (*iter)->vertexIndex0()) {
00297 found = true;
00298 (*iter)->setMother(currentParticle);
00299 currentParticle->addDaughter(*iter);
00300 AddMother((*iter), particleCol);
00301 }
00302 }
00303 if (!found) std::cout << "AddMother: inconsistency was found!" << std::endl;
00304
00305 }
00306
00307 void BesMcTruthWriter::SaveDecayMode()
00308 {
00309 SmartDataPtr<DecayMode> decayMode(m_evtSvc,"/Event/MC/DecayMode");
00310 if(!decayMode)
00311 {
00312
00313 DecayMode* aDecayMode = new DecayMode;
00314
00315 int dm[10]={0,0,0,0,0,0,0,0,0,0};
00316 aDecayMode->putData(dm,10);
00317 StatusCode scDM = m_evtSvc->registerObject("/Event/MC/DecayMode", aDecayMode);
00318 if(scDM!=StatusCode::SUCCESS)
00319 G4cout<< "Could not register DecayMode" <<G4endl;
00320 }
00321 }
00322
00323 void BesMcTruthWriter::SaveMdcTruth()
00324 {
00325
00326 Event::MdcMcHitCol* aMdcMcHitCol = new Event::MdcMcHitCol;
00327
00328 G4int HCID = -1;
00329 HCID = m_DigiMan->GetHitsCollectionID("BesMdcTruthCollection");
00330 if(HCID>0)
00331 {
00332 BesMdcHitsCollection* HC = 0;
00333 HC = (BesMdcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00334 G4int n_hit = HC->entries();
00335 if(n_hit>0)
00336 {
00337
00338 BesMdcHit* hit;
00339 vector<BesMdcHit*>* vecHC = HC->GetVector();
00340 for(int i=0;i<n_hit-1;i++)
00341 for(int j=i+1;j<n_hit;j++)
00342 if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
00343 {
00344 hit = (*vecHC)[i];
00345 (*vecHC)[i] = (*vecHC)[j];
00346 (*vecHC)[j] = hit;
00347 }
00348
00349
00350 for(G4int i=0;i<n_hit;i++)
00351 {
00352 hit = (*HC)[i];
00353 const Identifier ident = MdcID::wire_id ( hit->GetLayerNo(), hit->GetCellNo() );
00354 Event::MdcMcHit* mdcMcHit = new Event::MdcMcHit(ident, hit->GetTrackID(),
00355 hit->GetPos().x(), hit->GetPos().y(), hit->GetPos().z(),
00356 hit->GetDriftD(), hit->GetEdep(),hit->GetPosFlag() );
00357 aMdcMcHitCol->push_back(mdcMcHit);
00358 }
00359
00360 }
00361 }
00362
00363
00364 StatusCode scMdc = m_evtSvc->registerObject("/Event/MC/MdcMcHitCol", aMdcMcHitCol);
00365 if(scMdc!=StatusCode::SUCCESS)
00366 G4cout<< "Could not register MDC McTruth collection" <<G4endl;
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389 }
00390
00391 void BesMcTruthWriter::SaveTofTruth()
00392 {
00393
00394 Event::TofMcHitCol* aTofMcHitCol = new Event::TofMcHitCol;
00395
00396 G4int HCID = -1;
00397 HCID = m_DigiMan->GetHitsCollectionID("BesTofHitsList");
00398 if(HCID>0)
00399 {
00400 BesTofHitsCollection* HC = 0;
00401 HC = (BesTofHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00402 G4int n_hit = HC->entries();
00403 if(n_hit>0)
00404 {
00405
00406 BesTofHit* hit;
00407 vector<BesTofHit*>* vecHC = HC->GetVector();
00408 for( int i=0; i<n_hit-1; i++ ) {
00409 for( int j=i+1; j<n_hit; j++ ) {
00410 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex()) {
00411 hit = (*vecHC)[i];
00412 (*vecHC)[i] = (*vecHC)[j];
00413 (*vecHC)[j] = hit;
00414 }
00415 }
00416 }
00417
00418
00419 for( G4int i=0; i<n_hit; i++ ) {
00420 hit = (*HC)[i];
00421 Identifier ident;
00422 unsigned int barrel_ec = hit->GetPartId();
00423
00424 if( TofID::is_scin( barrel_ec ) ) {
00425 unsigned int scinNum = hit->GetScinNb();
00426 unsigned int layer = 0;
00427 if( TofID::is_barrel(barrel_ec) && scinNum>TofID::getPHI_BARREL_MAX()) {
00428 layer = 1;
00429 scinNum = scinNum - TofID::getPHI_BARREL_MAX() - 1;
00430 }
00431 ident = TofID::cell_id ( barrel_ec, layer, scinNum, 0);
00432 }
00433
00434 else {
00435 if( barrel_ec==3 || barrel_ec==4 ) {
00436 unsigned int endcap = 0;
00437 unsigned int module = hit->GetModule();
00438 unsigned int strip = hit->GetStrip();
00439 if( barrel_ec==4 ) {
00440 endcap = 1;
00441 }
00442 ident = TofID::cell_id( 3, endcap, module, strip, 0 );
00443 }
00444 }
00445 if( barrel_ec>=0 && barrel_ec<=4 ) {
00446 Event::TofMcHit* tofMcHit = new Event::TofMcHit( ident, hit->GetTrackIndex(),
00447 hit->GetPos().x(), hit->GetPos().y(), hit->GetPos().z(), hit->GetMomentum().x(),
00448 hit->GetMomentum().y(), hit->GetMomentum().z(), hit->GetTrackL(), hit->GetTime() );
00449 aTofMcHitCol->push_back(tofMcHit);
00450 }
00451 }
00452 }
00453 }
00454
00455
00456 StatusCode scTof = m_evtSvc->registerObject("/Event/MC/TofMcHitCol", aTofMcHitCol);
00457 if(scTof!=StatusCode::SUCCESS)
00458 G4cout<< "Could not register TOF McTruth collection" <<G4endl;
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485 }
00486
00487 void BesMcTruthWriter::SaveEmcTruth()
00488 {
00489
00490 Event::EmcMcHitCol* aEmcMcHitCol = new Event::EmcMcHitCol;
00491
00492 G4int fullMc = 1;
00493 if(fullMc==1) {
00494 G4int HCID = -1;
00495 HCID = m_DigiMan->GetHitsCollectionID("BesEmcTruthHitsList");
00496 if(HCID>0)
00497 {
00498 BesEmcTruthHitsCollection* HC = 0;
00499 HC = (BesEmcTruthHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00500 G4int n_hit = HC->entries();
00501 if(n_hit>0)
00502 {
00503
00504 BesEmcTruthHit* hit;
00505 vector<BesEmcTruthHit*>* vecHC = HC->GetVector();
00506 for(int i=0;i<n_hit-1;i++) {
00507 for(int j=i+1;j<n_hit;j++) {
00508 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex()) {
00509 hit = (*vecHC)[i];
00510 (*vecHC)[i] = (*vecHC)[j];
00511 (*vecHC)[j] = hit;
00512 }
00513 }
00514 }
00515
00516 for(G4int i=0;i<n_hit;i++)
00517 {
00518 BesEmcTruthHit* hit;
00519 hit = (*HC)[i];
00520
00521 std::map<Identifier,double> hitMap;
00522 std::map<Identifier,double>::const_iterator iHitMap;
00523 hitMap.clear();
00524 for(iHitMap=hit->Begin();iHitMap!=hit->End();iHitMap++) {
00525 hitMap[iHitMap->first]=iHitMap->second;
00526 }
00527
00528 Event::EmcMcHit* emcMcHit = new Event::EmcMcHit(hit->GetIdentify(), hit->GetTrackIndex(),
00529 hit->GetPosition().x(), hit->GetPosition().y(), hit->GetPosition().z(),
00530 hit->GetMomentum().x(), hit->GetMomentum().y(), hit->GetMomentum().z(),
00531 hit->GetEDep() );
00532
00533 emcMcHit->setHitEmc(hit->GetHitEmc());
00534 emcMcHit->setPDGCode(hit->GetPDGCode());
00535 emcMcHit->setPDGCharge(hit->GetPDGCharge());
00536 emcMcHit->setTime(hit->GetTime());
00537 emcMcHit->setHitMap(hitMap);
00538
00539 aEmcMcHitCol->push_back(emcMcHit);
00540 }
00541 }
00542 }
00543 } else {
00544 G4int HCID = -1;
00545 HCID = m_DigiMan->GetHitsCollectionID("BesEmcHitsList");
00546 if(HCID>0)
00547 {
00548 BesEmcHitsCollection* HC = 0;
00549 HC = (BesEmcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00550 G4int n_hit = HC->entries();
00551 if(n_hit>0)
00552 {
00553
00554 BesEmcHit* hit;
00555 vector<BesEmcHit*>* vecHC = HC->GetVector();
00556 for(int i=0;i<n_hit-1;i++)
00557 for(int j=i+1;j<n_hit;j++)
00558 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
00559 {
00560 hit = (*vecHC)[i];
00561 (*vecHC)[i] = (*vecHC)[j];
00562 (*vecHC)[j] = hit;
00563 }
00564
00565
00566 Event::EmcMcHitCol* aEmcMcHitCol = new Event::EmcMcHitCol;
00567
00568 for(G4int i=0;i<n_hit;i++)
00569 {
00570 hit = (*HC)[i];
00571 Identifier ident = EmcID::crystal_id ( hit->GetPartId(), hit->GetNumThetaCrystal(), hit->GetNumPhiCrystal() );
00572
00573 std::map<Identifier,double> hitMap;
00574 Event::EmcMcHit* emcMcHit = new Event::EmcMcHit(ident, hit->GetTrackIndex(),
00575 hit->GetPosCrystal().x(), hit->GetPosCrystal().y(), hit->GetPosCrystal().z(),
00576 hit->GetMomentum().x(), hit->GetMomentum().y(), hit->GetMomentum().z(),
00577 hit->GetEdepCrystal() );
00578
00579 aEmcMcHitCol->push_back(emcMcHit);
00580 }
00581 }
00582 }
00583 }
00584
00585
00586 StatusCode scEmc = m_evtSvc->registerObject("/Event/MC/EmcMcHitCol", aEmcMcHitCol);
00587 if(scEmc!=StatusCode::SUCCESS)
00588 G4cout<< "Could not register EMC McTruth collection" <<G4endl;
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614 }
00615
00616 void BesMcTruthWriter::SaveMucTruth()
00617 {
00618
00619 Event::MucMcHitCol* aMucMcHitCol = new Event::MucMcHitCol;
00620
00621 G4int HCID = -1;
00622 HCID = m_DigiMan->GetHitsCollectionID("BesMucHitsList");
00623 if(HCID>0)
00624 {
00625 BesMucHitsCollection* HC = 0;
00626 HC = (BesMucHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00627 G4int n_hit = HC->entries();
00628 if(n_hit>0)
00629 {
00630
00631 BesMucHit* hit;
00632 vector<BesMucHit*>* vecHC = HC->GetVector();
00633 for(int i=0;i<n_hit-1;i++)
00634 for(int j=i+1;j<n_hit;j++)
00635 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
00636 {
00637 hit = (*vecHC)[i];
00638 (*vecHC)[i] = (*vecHC)[j];
00639 (*vecHC)[j] = hit;
00640 }
00641
00642 for(G4int i=0;i<n_hit;i++)
00643 {
00644 hit = (*HC)[i];
00645 Identifier ident = MucID::channel_id ( hit->GetPart(), hit->GetSeg(), hit->GetGap(),
00646 hit->GetStrip() );
00647 Event::MucMcHit* mucMcHit = new Event::MucMcHit(ident, hit->GetTrackIndex(),
00648 hit->GetPos().x(), hit->GetPos().y(), hit->GetPos().z(), hit->GetMomentum().x(),
00649 hit->GetMomentum().y(), hit->GetMomentum().z() );
00650 aMucMcHitCol->push_back(mucMcHit);
00651 }
00652
00653 }
00654 }
00655
00656
00657 StatusCode scMuc = m_evtSvc->registerObject("/Event/MC/MucMcHitCol", aMucMcHitCol);
00658 if(scMuc!=StatusCode::SUCCESS)
00659 G4cout<< "Could not register MUC McTruth collection" <<G4endl;
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685 }
00686