#include <BesMcTruthWriter.hh>
Public Member Functions | |
| void | AddMother (Event::McParticle *, Event::McParticleCol *) |
| void | AddMother (Event::McParticle *, Event::McParticleCol *) |
| BesMcTruthWriter () | |
| BesMcTruthWriter () | |
| void | SaveDecayMode () |
| void | SaveDecayMode () |
| void | SaveEmcTruth () |
| void | SaveEmcTruth () |
| void | SaveMcParticle () |
| void | SaveMcParticle () |
| void | SaveMcTruth () |
| void | SaveMcTruth () |
| void | SaveMdcTruth () |
| void | SaveMdcTruth () |
| void | SaveMucTruth () |
| void | SaveMucTruth () |
| void | SaveTofTruth () |
| void | SaveTofTruth () |
| ~BesMcTruthWriter () | |
| ~BesMcTruthWriter () | |
Private Attributes | |
| G4DigiManager * | m_DigiMan |
| G4DigiManager * | m_DigiMan |
| IDataProviderSvc * | m_evtSvc |
| IDataProviderSvc * | m_evtSvc |
| BesMdcGeoParameter * | mdcGeoPointer |
| BesMdcGeoParameter * | mdcGeoPointer |
|
|
00050 {
00051 m_DigiMan = G4DigiManager::GetDMpointer();
00052 //mdcGeoPointer=BesMdcGeoParameter::GetGeo();
00053 }
|
|
|
00056 {
00057 }
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
||||||||||||
|
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 }
|
|
|
|
|
|
00308 {
00309 SmartDataPtr<DecayMode> decayMode(m_evtSvc,"/Event/MC/DecayMode");
00310 if(!decayMode)
00311 {
00312 //G4cout<<"Could not retrieve DecayMode"<<G4endl;
00313 DecayMode* aDecayMode = new DecayMode;
00314 //register Decay Mode to TDS
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 }
|
|
|
|
|
|
00471 {
00472 //Emc McHit collection defined in BOSS
00473 Event::EmcMcHitCol* aEmcMcHitCol = new Event::EmcMcHitCol;
00474
00475 G4int fullMc = 1;
00476 if(fullMc==1) { //full Mc Truth
00477 G4int HCID = -1;
00478 HCID = m_DigiMan->GetHitsCollectionID("BesEmcTruthHitsList");
00479 if(HCID>0)
00480 {
00481 BesEmcTruthHitsCollection* HC = 0;
00482 HC = (BesEmcTruthHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00483 G4int n_hit = HC->entries();
00484 if(n_hit>0)
00485 {
00486 //arrange hits in hits collection in order of trackIndex
00487 BesEmcTruthHit* hit;
00488 vector<BesEmcTruthHit*>* vecHC = HC->GetVector();
00489 for(int i=0;i<n_hit-1;i++) {
00490 for(int j=i+1;j<n_hit;j++) {
00491 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex()) {
00492 hit = (*vecHC)[i];
00493 (*vecHC)[i] = (*vecHC)[j];
00494 (*vecHC)[j] = hit;
00495 }
00496 }
00497 }
00498
00499 for(G4int i=0;i<n_hit;i++)
00500 {
00501 BesEmcTruthHit* hit;
00502 hit = (*HC)[i];
00503
00504 std::map<Identifier,double> hitMap;
00505 std::map<Identifier,double>::const_iterator iHitMap;
00506 hitMap.clear();
00507 for(iHitMap=hit->Begin();iHitMap!=hit->End();iHitMap++) {
00508 hitMap[iHitMap->first]=iHitMap->second;
00509 }
00510
00511 Event::EmcMcHit* emcMcHit = new Event::EmcMcHit(hit->GetIdentify(), hit->GetTrackIndex(),
00512 hit->GetPosition().x(), hit->GetPosition().y(), hit->GetPosition().z(),
00513 hit->GetMomentum().x(), hit->GetMomentum().y(), hit->GetMomentum().z(),
00514 hit->GetEDep() );
00515
00516 emcMcHit->setHitEmc(hit->GetHitEmc());
00517 emcMcHit->setPDGCode(hit->GetPDGCode());
00518 emcMcHit->setPDGCharge(hit->GetPDGCharge());
00519 emcMcHit->setTime(hit->GetTime());
00520 emcMcHit->setHitMap(hitMap);
00521
00522 aEmcMcHitCol->push_back(emcMcHit);
00523 }
00524 }
00525 }
00526 } else { //simple Mc Truth
00527 G4int HCID = -1;
00528 HCID = m_DigiMan->GetHitsCollectionID("BesEmcHitsList");
00529 if(HCID>0)
00530 {
00531 BesEmcHitsCollection* HC = 0;
00532 HC = (BesEmcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00533 G4int n_hit = HC->entries();
00534 if(n_hit>0)
00535 {
00536 //arrange hits in hits collection in order of trackIndex
00537 BesEmcHit* hit;
00538 vector<BesEmcHit*>* vecHC = HC->GetVector();
00539 for(int i=0;i<n_hit-1;i++)
00540 for(int j=i+1;j<n_hit;j++)
00541 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
00542 {
00543 hit = (*vecHC)[i];
00544 (*vecHC)[i] = (*vecHC)[j];
00545 (*vecHC)[j] = hit;
00546 }
00547
00548 //Emc McHit collection defined in BOSS
00549 Event::EmcMcHitCol* aEmcMcHitCol = new Event::EmcMcHitCol;
00550
00551 for(G4int i=0;i<n_hit;i++)
00552 {
00553 hit = (*HC)[i];
00554 Identifier ident = EmcID::crystal_id ( hit->GetPartId(), hit->GetNumThetaCrystal(), hit->GetNumPhiCrystal() );
00555
00556 std::map<Identifier,double> hitMap;
00557 Event::EmcMcHit* emcMcHit = new Event::EmcMcHit(ident, hit->GetTrackIndex(),
00558 hit->GetPosCrystal().x(), hit->GetPosCrystal().y(), hit->GetPosCrystal().z(),
00559 hit->GetMomentum().x(), hit->GetMomentum().y(), hit->GetMomentum().z(),
00560 hit->GetEdepCrystal() );
00561
00562 aEmcMcHitCol->push_back(emcMcHit);
00563 }
00564 }
00565 }
00566 }
00567
00568 //register EMC McTruth collection to TDS
00569 StatusCode scEmc = m_evtSvc->registerObject("/Event/MC/EmcMcHitCol", aEmcMcHitCol);
00570 if(scEmc!=StatusCode::SUCCESS)
00571 G4cout<< "Could not register EMC McTruth collection" <<G4endl;
00572
00573 //retrieve EMC McTruth from TDS
00574 /*SmartDataPtr<Event::EmcMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/EmcMcHitCol");
00575 if(!aMcHitCol)
00576 G4cout<<"Could not retrieve EMC McTruth collection"<<G4endl;
00577
00578 Event::EmcMcHitCol::iterator iMcHitCol;
00579 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
00580 {
00581 const Identifier ident = (*iMcHitCol)->identify();
00582 G4cout<<(*iMcHitCol)->getTrackIndex();
00583 G4cout<<" "<<EmcID::barrel_ec(ident);
00584 G4cout<<" "<<EmcID::theta_module(ident);
00585 G4cout<<" "<<EmcID::phi_module(ident);
00586 G4cout<<" "<<(*iMcHitCol)->getPositionX();
00587 G4cout<<" "<<(*iMcHitCol)->getPositionY();
00588 G4cout<<" "<<(*iMcHitCol)->getPositionZ();
00589 G4cout<<" "<<(*iMcHitCol)->getPx();
00590 G4cout<<" "<<(*iMcHitCol)->getPy();
00591 G4cout<<" "<<(*iMcHitCol)->getPz();
00592 G4cout<<" "<<(*iMcHitCol)->getDepositEnergy();
00593 G4cout<<G4endl;
00594 }
00595 G4cout<<"end of retrieve EMC McTruth"<<G4endl;
00596 */
00597 }
|
|
|
|
|
|
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 //arrange TruthTrack in trackList in order of trackIndex
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 //loop over tracks
00108 for(int i=0;i<nTrack;i++)
00109 {
00110 track = (*trackList) [i];
00111
00112 // find out the start point
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 // find out the end point
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 //create the mc particle
00143 Event::McParticle* mcParticle = new Event::McParticle;
00144
00145 // initialization
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 // Set track index
00153 mcParticle->setTrackIndex( track->GetIndex() );
00154
00155 // Adding status flag
00156 if (isPrimary ) mcParticle->addStatusFlag(Event::McParticle::PRIMARY);
00157
00158 if (track->GetDaughterIndexes().size()==0)
00159 mcParticle->addStatusFlag(Event::McParticle::LEAF);
00160
00161 //std::cout<<"pdg:"<<track->GetPDGCode()<<" "
00162 // <<"source:"<<track->GetSource()<<" "
00163 // <<std::endl;
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 //std::cout<<"status:"<<mcParticle->statusFlags()<<" "
00172 // <<"FGN:"<<mcParticle->decayFromGenerator()<<" "
00173 // <<"FG4:"<<mcParticle->decayInFlight()<<std::endl;
00174 // Adding vertex index
00175 mcParticle->setVertexIndex0(startPoint->GetIndex());
00176 mcParticle->setVertexIndex1(vertex->GetIndex());
00177
00178 // Set the final position
00179 HepLorentzVector finalPosition(vertex->GetPosition().x()/10., vertex->GetPosition().y()/10., vertex->GetPosition().z()/10., vertex->GetTime());
00180 mcParticle->finalize(finalPosition);
00181
00182 // push back
00183 particleCol->push_back(mcParticle);
00184
00185 endPointFound = true;
00186 }
00187 }
00188 }
00189
00190 if (!endPointFound)
00191 {
00192 //std::cout << "warning in finding the end point of a track" <<std::endl;
00193 if (track->GetDaughterIndexes().size()==0)
00194 {
00195 // create the mc particle
00196 Event::McParticle* mcParticle = new Event::McParticle;
00197 // initialization
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 // Set track index
00204 mcParticle->setTrackIndex( track->GetIndex() );
00205
00206 // Adding status flag
00207 mcParticle->addStatusFlag(Event::McParticle::LEAF);
00208 if (isPrimary ) mcParticle->addStatusFlag(Event::McParticle::PRIMARY);
00209
00210 // Adding vertex index
00211 mcParticle->setVertexIndex0( startPoint->GetIndex() );
00212 mcParticle->setVertexIndex1( -99 );
00213
00214 //std::cout<<"pdg:"<<track->GetPDGCode()<<" "
00215 // <<"source:"<<track->GetSource()<<" "
00216 // <<std::endl;
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 //std::cout<<"status:"<<mcParticle->statusFlags()<<" "
00225 // <<"FGN:"<<mcParticle->decayFromGenerator()<<" "
00226 // <<"FG4:"<<mcParticle->decayInFlight()<<std::endl;
00227
00228 // push back
00229 particleCol->push_back(mcParticle);
00230 continue;
00231 }
00232 }
00233 }
00234
00235 // Get primary McParticles
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 // Add mother particle recursively
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 //register McParticle collection to TDS
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 //retrive McParticle from TDS
00264 /*SmartDataPtr<Event::McParticleCol> mcParticleCol(m_evtSvc,"/Event/MC/McParticleCol");
00265 if(!mcParticleCol)
00266 G4cout<<"Could not retrieve McParticelCol"<<G4endl;
00267 Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
00268 for (;iter_mc != mcParticleCol->end(); iter_mc++)
00269 {
00270 G4cout<<(*iter_mc)->getTrackIndex();
00271 G4cout<<" "<<(*iter_mc)->particleProperty();
00272 G4cout<<" "<<(*iter_mc)->getVertexIndex0();
00273 G4cout<<" "<<(*iter_mc)->getVertexIndex1();
00274 G4cout<<" "<<(*iter_mc)->initialFourMomentum().x();
00275 G4cout<<" "<<(*iter_mc)->initialFourMomentum().y();
00276 G4cout<<" "<<(*iter_mc)->initialFourMomentum().z();
00277 G4cout<<" "<<(*iter_mc)->initialFourMomentum().t();
00278 G4cout<<" "<<(*iter_mc)->initialPosition().x();
00279 G4cout<<" "<<(*iter_mc)->initialPosition().y();
00280 G4cout<<" "<<(*iter_mc)->initialPosition().z();
00281 G4cout<<G4endl;
00282 }
00283 G4cout<<"end of retrieve McParticle from TDS"<<G4endl;
00284 */
00285
00286 }
|
|
|
|
|
|
00060 {
00061 //interface to event data service
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 }
|
|
|
|
|
|
00324 {
00325 //Mdc McHit collection defined in BOSS
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 //arrange hits in hits collection in order of trackIndex
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 //push back MdcMcHits to MdcMcHitCol in TDS
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 //register MDC McTruth collection to TDS
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 //retrieve MDC McTruth from TDS
00369 /*SmartDataPtr<Event::MdcMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/MdcMcHitCol");
00370 if(!aMcHitCol)
00371 G4cout<<"Could not retrieve MDC McTruth collection"<<G4endl;
00372
00373 Event::MdcMcHitCol::iterator iMcHitCol;
00374 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
00375 {
00376 const Identifier ident = (*iMcHitCol)->identify();
00377 G4cout<<(*iMcHitCol)->getTrackIndex();
00378 G4cout<<" "<<MdcID::layer(ident);
00379 G4cout<<" "<<MdcID::wire(ident);
00380 G4cout<<" "<<(*iMcHitCol)->getDepositEnergy();
00381 G4cout<<" "<<(*iMcHitCol)->getDriftDistance();
00382 G4cout<<" "<<(*iMcHitCol)->getPositionX();
00383 G4cout<<" "<<(*iMcHitCol)->getPositionY();
00384 G4cout<<" "<<(*iMcHitCol)->getPositionZ();
00385 G4cout<<G4endl;
00386 }
00387 G4cout<<"end of retrieve MDC McTruth collection"<<G4endl;
00388 */
00389 }
|
|
|
|
|
|
00600 {
00601 //Muc McHit collection defined in BOSS
00602 Event::MucMcHitCol* aMucMcHitCol = new Event::MucMcHitCol;
00603
00604 G4int HCID = -1;
00605 HCID = m_DigiMan->GetHitsCollectionID("BesMucHitsList");
00606 if(HCID>0)
00607 {
00608 BesMucHitsCollection* HC = 0;
00609 HC = (BesMucHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00610 G4int n_hit = HC->entries();
00611 if(n_hit>0)
00612 {
00613 //arrange hits in hits collection in order of trackIndex
00614 BesMucHit* hit;
00615 vector<BesMucHit*>* vecHC = HC->GetVector();
00616 for(int i=0;i<n_hit-1;i++)
00617 for(int j=i+1;j<n_hit;j++)
00618 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
00619 {
00620 hit = (*vecHC)[i];
00621 (*vecHC)[i] = (*vecHC)[j];
00622 (*vecHC)[j] = hit;
00623 }
00624
00625 for(G4int i=0;i<n_hit;i++)
00626 {
00627 hit = (*HC)[i];
00628 Identifier ident = MucID::channel_id ( hit->GetPart(), hit->GetSeg(), hit->GetGap(),
00629 hit->GetStrip() );
00630 Event::MucMcHit* mucMcHit = new Event::MucMcHit(ident, hit->GetTrackIndex(),
00631 hit->GetPos().x(), hit->GetPos().y(), hit->GetPos().z(), hit->GetMomentum().x(),
00632 hit->GetMomentum().y(), hit->GetMomentum().z() );
00633 aMucMcHitCol->push_back(mucMcHit);
00634 }
00635
00636 }
00637 }
00638
00639 //register MUC McTruth collection to TDS
00640 StatusCode scMuc = m_evtSvc->registerObject("/Event/MC/MucMcHitCol", aMucMcHitCol);
00641 if(scMuc!=StatusCode::SUCCESS)
00642 G4cout<< "Could not register MUC McTruth collection" <<G4endl;
00643
00644 //retrieve MUC McTruth from TDS
00645 /*SmartDataPtr<Event::MucMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/MucMcHitCol");
00646 if(!aMcHitCol)
00647 G4cout<<"Could not retrieve MUC McTruth collection"<<G4endl;
00648
00649 Event::MucMcHitCol::iterator iMcHitCol;
00650 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
00651 {
00652 const Identifier ident = (*iMcHitCol)->identify();
00653 G4cout<<(*iMcHitCol)->getTrackIndex();
00654 G4cout<<" "<<MucID::part(ident);
00655 G4cout<<" "<<MucID::seg(ident);
00656 G4cout<<" "<<MucID::gap(ident);
00657 G4cout<<" "<<MucID::strip(ident);
00658 G4cout<<" "<<(*iMcHitCol)->getPositionX();
00659 G4cout<<" "<<(*iMcHitCol)->getPositionY();
00660 G4cout<<" "<<(*iMcHitCol)->getPositionZ();
00661 G4cout<<" "<<(*iMcHitCol)->getPx();
00662 G4cout<<" "<<(*iMcHitCol)->getPy();
00663 G4cout<<" "<<(*iMcHitCol)->getPz();
00664 G4cout<<G4endl;
00665 }
00666 G4cout<<"end of retrieve MUC McTruth"<<G4endl;
00667 */
00668 }
|
|
|
|
|
|
00392 {
00393 //Tof McHit collection defined in BOSS
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 //arrange hits in hits collection in order of trackIndex
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 {
00412 hit = (*vecHC)[i];
00413 (*vecHC)[i] = (*vecHC)[j];
00414 (*vecHC)[j] = hit;
00415 }
00416
00417 //push back TofMcHit collection to TDS
00418 for(G4int i=0;i<n_hit;i++)
00419 {
00420 hit = (*HC)[i];
00421 unsigned int scinNum, barrel_ec, layer = 0;
00422 scinNum = hit->GetScinNb();
00423 barrel_ec = hit->GetPartId();
00424 if (TofID::is_barrel(barrel_ec) && scinNum > TofID::getPHI_BARREL_MAX()) {
00425 layer = 1;
00426 scinNum = scinNum - TofID::getPHI_BARREL_MAX() - 1;
00427 }
00428 Identifier ident = TofID::cell_id ( barrel_ec, layer, scinNum, 0);
00429 Event::TofMcHit* tofMcHit = new Event::TofMcHit( ident, hit->GetTrackIndex(),
00430 hit->GetPos().x(), hit->GetPos().y(), hit->GetPos().z(), hit->GetMomentum().x(),
00431 hit->GetMomentum().y(), hit->GetMomentum().z(), hit->GetTrackL(), hit->GetTime() );
00432
00433 aTofMcHitCol->push_back(tofMcHit);
00434 }
00435 }
00436 }
00437
00438 //register TOF McTruth collection to TDS
00439 StatusCode scTof = m_evtSvc->registerObject("/Event/MC/TofMcHitCol", aTofMcHitCol);
00440 if(scTof!=StatusCode::SUCCESS)
00441 G4cout<< "Could not register TOF McTruth collection" <<G4endl;
00442
00443 //retrieve TOF McTruth from TDS
00444 /*SmartDataPtr<Event::TofMcHitCol> aMcHitCol(m_evtSvc,"/Event/MC/TofMcHitCol");
00445 if(!aMcHitCol)
00446 G4cout<<"Could not retrieve TOF McTruth collection"<<G4endl;
00447
00448 Event::TofMcHitCol::iterator iMcHitCol;
00449 for(iMcHitCol=aMcHitCol->begin(); iMcHitCol!=aMcHitCol->end(); iMcHitCol++)
00450 {
00451 const Identifier ident = (*iMcHitCol)->identify();
00452 G4cout<<(*iMcHitCol)->getTrackIndex();
00453 G4cout<<" "<<TofID::barrel_ec(ident);;
00454 G4cout<<" "<<TofID::layer(ident);
00455 G4cout<<" "<<TofID::phi_module(ident);
00456 G4cout<<" "<<(*iMcHitCol)->getPositionX();
00457 G4cout<<" "<<(*iMcHitCol)->getPositionY();
00458 G4cout<<" "<<(*iMcHitCol)->getPositionZ();
00459 G4cout<<" "<<(*iMcHitCol)->getPx();
00460 G4cout<<" "<<(*iMcHitCol)->getPy();
00461 G4cout<<" "<<(*iMcHitCol)->getPz();
00462 G4cout<<" "<<(*iMcHitCol)->getTrackLength();
00463 G4cout<<" "<<(*iMcHitCol)->getFlightTime();
00464 G4cout<<G4endl;
00465 }
00466 G4cout<<"end of retrieve TOF McTruth"<<G4endl;
00467 */
00468 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1