00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00015
00016
00017 #include "BesSensitiveManager.hh"
00018 #include "BesTruthTrack.hh"
00019 #include "BesTruthVertex.hh"
00020
00021 #include "G4TrackingManager.hh"
00022 #include "G4Track.hh"
00023 #include "G4VProcess.hh"
00024 #include "G4Event.hh"
00025
00026 #include "HepMC/GenEvent.h"
00027
00028 #include "GaudiKernel/ISvcLocator.h"
00029 #include "GaudiKernel/Bootstrap.h"
00030 #include "GaudiKernel/MsgStream.h"
00031 #include "GaudiKernel/IMessageSvc.h"
00032 #include "GaudiKernel/IDataProviderSvc.h"
00033
00034 #include "GeneratorObject/McGenEvent.h"
00035 #include "GaudiKernel/SmartDataPtr.h"
00036
00037 BesSensitiveManager* BesSensitiveManager::m_sensitiveManager=0;
00038
00039
00040
00041
00042
00043 void BesSensitiveManager::BeginOfTruthEvent(const G4Event* evt)
00044 {
00045
00046
00047 m_trackList = new std::vector<BesTruthTrack*>;
00048 m_vertexList = new std::vector<BesTruthVertex*>;
00049
00050 SetVertex0(evt);
00051
00052 m_count = 0;
00053 SaveParticlesFromGenerator();
00054 m_count += m_trackList->size();
00055
00056
00057 iter=clients.begin();
00058 iter_end=clients.end();
00059 while( iter != iter_end )
00060 {
00061 (*iter)->BeginOfTruthEvent(evt);
00062 iter++;
00063 }
00064 }
00065
00066 void BesSensitiveManager::SetVertex0(const G4Event* anEvent)
00067 {
00068
00069 for( G4int i=0; i<1; i++ )
00070 {
00071 G4PrimaryVertex* primaryVertex = anEvent->GetPrimaryVertex(i);
00072 m_pos0 = primaryVertex->GetPosition();
00073 m_t0 = primaryVertex->GetT0();
00074 }
00075 if(m_logLevel<=4)
00076 G4cout<<"pos0:"<<m_pos0<<" t0:"<<m_t0<<G4endl;
00077 }
00078
00079 G4int BesSensitiveManager::CheckType(const HepMC::GenEvent* hepmcevt)
00080 {
00081 G4int flag =0;
00082 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00083 vitr != hepmcevt->vertices_end(); ++vitr ) {
00084 for (HepMC::GenVertex::particle_iterator
00085 pitr= (*vitr)->particles_begin(HepMC::children);
00086 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
00087 if((*pitr)->end_vertex())
00088 {flag = 1; break;}
00089 }
00090 if(flag) break;
00091 }
00092 if(m_logLevel <= 4)
00093 {
00094 if(flag == 0)
00095 G4cout<<G4endl<<"generator is GENBES!"<<G4endl;
00096 else
00097 G4cout<<G4endl<<"generator is not GENBES!"<<G4endl;
00098 }
00099 return flag;
00100 }
00101
00102 void BesSensitiveManager::SaveParticlesFromGenerator()
00103 {
00104 IDataProviderSvc* p_evtSvc=0;
00105 ISvcLocator* svcLocator = Gaudi::svcLocator();
00106 StatusCode sc=svcLocator->service("EventDataSvc", p_evtSvc);
00107 if (sc.isFailure())
00108 std::cout<<"BesHepMCInterface could not access EventDataSvc!!"<<std::endl;
00109
00110 SmartDataPtr<McGenEventCol> mcCollptr( p_evtSvc, "/Event/Gen");
00111
00112 if ( mcCollptr == 0 )
00113 std::cout << "no McGenEventCollection found." << std::endl;
00114
00115 else
00116 {
00117 McGenEventCol::const_iterator it = mcCollptr->begin();
00118 McGenEvent* mcEvent = (McGenEvent* ) (*it);
00119 m_hepmcevt = mcEvent->getGenEvt();
00120 G4int typeFlag = CheckType(m_hepmcevt);
00121
00122
00123 for(HepMC::GenEvent::vertex_const_iterator vitr= m_hepmcevt->vertices_begin();
00124 vitr != m_hepmcevt->vertices_end(); ++vitr )
00125 {
00126 G4int barcodeVtx = (*vitr)->barcode();
00127 if(m_logLevel<=4)
00128 G4cout<<G4endl<<"barcodeVtx:"<<barcodeVtx<<" ";
00129
00130 G4LorentzVector xvtx((*vitr)->position().x(),(*vitr)->position().y(),(*vitr)->position().z(),(*vitr)->position().t());
00131 G4ThreeVector v(xvtx.x(), xvtx.y(), xvtx.z());
00132 BesTruthVertex* newVertex = new BesTruthVertex();
00133 newVertex->SetPosition(v);
00134 newVertex->SetTime(xvtx.t()*mm/c_light);
00135 if(m_logLevel<=4)
00136 G4cout<<"xyzt:"<<xvtx.x()<<" "<<xvtx.y()<<" "
00137 <<xvtx.z()<<" "<<xvtx.t()*mm/c_light<<" ";
00138 G4int nTrack = m_trackList->size();
00139 BesTruthTrack* parentTrack=0;
00140 G4int parentTrackIndex= -99;
00141 for(int i=0;i<nTrack;i++)
00142 {
00143 parentTrack = (*m_trackList)[i];
00144 if(parentTrack->GetBarcodeEndVtx() == barcodeVtx)
00145 {
00146 newVertex->SetParentTrack(parentTrack);
00147 parentTrackIndex = parentTrack->GetIndex();
00148 if(m_logLevel<=4)
00149 G4cout<<"parentTrack:"<<parentTrackIndex<<" ";
00150 parentTrack->SetTerminalVertex(newVertex);
00151
00152 }
00153 }
00154
00155 G4int vtxFlag=0;
00156 G4int pOut = (*vitr)->particles_out_size();
00157 HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
00158 G4int pdgcode= (*pitr)-> pdg_id();
00159 if(pOut == 1 && pdgcode == -11 && typeFlag==1)
00160 vtxFlag=1;
00161
00162 if(vtxFlag==0)
00163 {
00164 m_vertexList->push_back(newVertex);
00165 newVertex->SetIndex(m_vertexList->size()-1);
00166 if(m_logLevel<=4)
00167 G4cout<<"vindex:"<<m_vertexList->size()-1<<G4endl;
00168
00169 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children); pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
00170 {
00171 G4LorentzVector p((*pitr)->momentum().x(), (*pitr)->momentum().y(),(*pitr)->momentum().z(),(*pitr)->momentum().t());
00172 G4LorentzVector pGeV(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
00173 G4int pdgcode = (*pitr)->pdg_id();
00174
00175 BesTruthTrack* newTrack = new BesTruthTrack;
00176 newTrack->SetP4(pGeV);
00177 newTrack->SetPDGCode(pdgcode);
00178 if(m_logLevel<=4)
00179 {
00180 G4cout<<"pdg:"<<pdgcode<<" ";
00181 G4cout<<"p:"<<p.x()<<" "<<p.y()<<" "<<p.z()<<" "<<p.t()<<" ";
00182 G4double mass = sqrt(p.t()*p.t()-p.x()*p.x()-p.y()*p.y()-p.z()*p.z());
00183 G4cout<<mass<<" ";
00184 }
00185 newTrack->SetPDGCharge(99);
00186 newTrack->SetParticleName("unknown");
00187 newTrack->SetVertex(newVertex);
00188 newTrack->SetTerminalVertex(0);
00189 newTrack->SetSource("FromGenerator");
00190
00191 G4int barcodeEndVtx=0;
00192 if((*pitr)->end_vertex())
00193 {
00194 barcodeEndVtx = (*pitr)->end_vertex()->barcode();
00195 if(m_logLevel<=4)
00196 G4cout<<"endVtx:"<<barcodeEndVtx<<" ";
00197 }
00198 newTrack->SetBarcodeEndVtx(barcodeEndVtx);
00199
00200 m_trackList->push_back(newTrack);
00201 newTrack->SetIndex(m_trackList->size()-1);
00202 if(m_logLevel<=4)
00203 G4cout<<"index:"<<m_trackList->size()-1<<" ";
00204
00205
00206 if(parentTrackIndex>=0)
00207 {
00208 if(m_logLevel<=4)
00209 G4cout<<"mother:"<<parentTrackIndex;
00210 (*m_trackList)[parentTrackIndex]->AddDaughterIndex(m_trackList->size()-1);
00211 }
00212 if(m_logLevel<=4)
00213 G4cout<<G4endl;
00214 }
00215 }
00216 }
00217 }
00218 }
00219
00220
00221
00222 void BesSensitiveManager::EndOfTruthEvent(const G4Event* evt)
00223 {
00224
00225 iter=clients.begin();
00226 iter_end=clients.end();
00227 while( iter != iter_end )
00228 {
00229 (*iter)->EndOfTruthEvent(evt);
00230 iter++;
00231 }
00232 }
00233
00234 void BesSensitiveManager::ClearEvent()
00235 {
00236 if(m_trackList)
00237 {
00238 for(size_t i=0;i<m_trackList->size();i++)
00239 delete (*m_trackList)[i];
00240 m_trackList->clear();
00241 delete m_trackList;
00242 }
00243 if(m_vertexList)
00244 {
00245 for(size_t i=0;i<m_vertexList->size();i++)
00246 delete (*m_vertexList)[i];
00247 m_vertexList->clear();
00248 delete m_vertexList;
00249 }
00250 }
00251
00252 void BesSensitiveManager::BeginOfTrack( const G4Track *track )
00253 {
00254 m_trackFlag = 0;
00255
00256
00257 chain.push_back(FollowTrack(track));
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270
00271 iter=clients.begin();
00272 iter_end=clients.end();
00273 while(iter!=iter_end)
00274 {
00275 (*iter)->BeginOfTrack( track );
00276 iter++;
00277 }
00278 }
00279 void BesSensitiveManager::EndOfTrack( const G4Track *track , G4TrackingManager* trackingManager )
00280 {
00281 if(chain.back().savedByDefault==true)
00282 {
00283 BesTStats &stat = chain.back();
00284 G4int endVtxFlag = 0;
00285 if(m_trackFlag==1)
00286 {
00287 BesTruthTrack* truthTrack = (*m_trackList)[m_trackIndex];
00288 if(truthTrack->GetTerminalVertex())
00289 {
00290 UpdateVertex(stat,track); endVtxFlag = 1;}
00291 }
00292 if(endVtxFlag == 0)
00293 {
00294
00295 G4StepPoint *finalStep = track->GetStep()->GetPostStepPoint();
00296 G4StepStatus stepStatus = finalStep->GetStepStatus();
00297 if (track->GetNextVolume() != 0 ||
00298 (stepStatus != fGeomBoundary && stepStatus != fWorldBoundary) )
00299 {
00300 if(m_logLevel<=4)
00301 G4cout<<"Particle died. make a terminal vertex: ";
00302
00303 const G4ThreeVector pos(finalStep->GetPosition());
00304 G4int vindex = m_vertexList->size();
00305 if(m_logLevel<=4)
00306 G4cout<<vindex<<G4endl;
00307 stat.vertices.push_back(vindex);
00308 BesTruthVertex *newVertex = new BesTruthVertex();
00309 newVertex->SetPosition(pos);
00310 newVertex->SetTime( finalStep->GetGlobalTime());
00311 if(m_trackFlag==1)
00312 newVertex->SetParentTrack( (*m_trackList)[m_trackIndex] );
00313 else
00314 newVertex->SetParentTrack( m_trackList->back() );
00315 newVertex->SetTerminal( true );
00316 newVertex->SetIndex( vindex );
00317
00318
00319 newVertex->SetMinDau(m_count);
00320
00321
00322 G4TrackVector* trackVector = trackingManager->GimmeSecondaries();
00323 G4int nSecon = trackVector->size();
00324 G4Track* seconTrack;
00325 G4int nDau=0;
00326 if(nSecon>0)
00327 {
00328 for(G4int i=0;i<nSecon;i++)
00329 {
00330 seconTrack = (*trackVector)[i];
00331 G4String processName = seconTrack->GetCreatorProcess()->GetProcessName();
00332 if(processName == "Decay")
00333 nDau++;
00334 }
00335 }
00336 m_vertexList->push_back( newVertex );
00337
00338 m_count += nDau;
00339
00340
00341
00342
00343
00344
00345 if(m_trackFlag==1)
00346 (*m_trackList)[m_trackIndex]->SetTerminalVertex( m_vertexList->back() );
00347 else
00348 m_trackList->back()->SetTerminalVertex( m_vertexList->back() );
00349 }
00350 }
00351 }
00352
00353 iter=clients.begin();
00354 iter_end=clients.end();
00355 while( iter!=iter_end) {
00356 (*iter)->EndOfTrack( track );
00357 iter++;
00358 }
00359 }
00360
00361 void BesSensitiveManager::UpdateVertex(BesTStats stat, const G4Track *track)
00362 {
00363 BesTruthTrack* truthTrack = (*m_trackList)[m_trackIndex];
00364 G4int terminalIndex = truthTrack->GetTerminalVertex()->GetIndex();
00365 if(m_logLevel<=4)
00366 G4cout<<"updateVertex:"<<terminalIndex<<" ";
00367 BesTruthVertex* vertex = (*m_vertexList)[terminalIndex];
00368 G4StepPoint *finalStep = track->GetStep()->GetPostStepPoint();
00369 const G4ThreeVector pos(finalStep->GetPosition());
00370 G4double time = finalStep->GetGlobalTime();
00371 if(m_logLevel<=4)
00372 G4cout<<"newPos:"<<pos<<" "<<"newTime:"<<time<<G4endl;
00373 vertex->SetPosition(pos);
00374 vertex->SetTime(time);
00375 vertex->SetTerminal( true );
00376
00377
00378
00379 }
00380
00381 void BesSensitiveManager::MakeNewTrack( BesTStats &stat, const G4Track *track)
00382 {
00383 if (stat.originVertex < 0 && chain.size() > 0)
00384 {
00385 if(m_logLevel<=4)
00386 G4cout<<"MakeNewTrack for decayed particles,";
00387
00388
00389
00390 G4int parentVstat = chain.back().vertices.size();
00391 while( --parentVstat >= 0)
00392 {
00393 G4int vindex = chain.back().vertices[parentVstat];
00394
00395 G4ThreeVector diff((*m_vertexList)[vindex]->GetPosition());
00396 diff = diff - track->GetPosition();
00397 if (diff.mag2() < 1E-9)
00398 {
00399 stat.originVertex = vindex;
00400 if(m_logLevel<=4)
00401 G4cout<<"find a vertex:"<<vindex;
00402 (*m_vertexList)[vindex]->AddCurrentDau();
00403 break;
00404 }
00405 }
00406
00407 }
00408
00409 if(stat.originVertex < 0 && chain.size() == 0)
00410 {
00411
00412 G4int nVertex = m_vertexList->size();
00413 for(G4int i=0;i<nVertex;i++)
00414 {
00415 G4ThreeVector diff((*m_vertexList)[i]->GetPosition());
00416 diff = diff -track->GetPosition();
00417 if(diff.mag2() < 1E-9)
00418 {
00419 stat.originVertex = i;
00420 if(m_logLevel<=4)
00421 G4cout<<"MakeNewTrack for primary particles,find a vertex:"
00422 <<i;
00423 break;
00424 }
00425 }
00426 }
00427
00428 if (stat.originVertex < 0)
00429 {
00430 if(m_logLevel<=4)
00431 G4cout<<"MakeNewTrack for primary particles,make new vertex:"
00432 <<m_vertexList->size();
00433
00434
00435 const G4ThreeVector v(track->GetPosition());
00436
00437 stat.originVertex = m_vertexList->size();
00438
00439 BesTruthVertex *newVertex = new BesTruthVertex();
00440 newVertex->SetPosition(v);
00441 newVertex->SetTime( track->GetGlobalTime());
00442
00443
00444
00445
00446
00447
00448
00449
00450 newVertex->SetParentTrack( 0 );
00451 newVertex->SetTerminal( false );
00452 newVertex->SetIndex( stat.originVertex );
00453 if(track->GetCreatorProcess())
00454 newVertex->SetProcessName(track->GetCreatorProcess()->GetProcessName());
00455 else
00456 newVertex->SetProcessName("generator");
00457
00458 m_vertexList->push_back( newVertex );
00459 }
00460
00461
00462
00463
00464 BesTruthTrack *newTrack = new BesTruthTrack();
00465 newTrack->SetP4( stat.p4 );
00466 newTrack->SetPDGCode(track->GetDefinition()->GetPDGEncoding());
00467 newTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
00468 newTrack->SetParticleName(track->GetDefinition()->GetParticleName());
00469 newTrack->SetSource("FromG4");
00470 BesTruthVertex* vertex = (*m_vertexList)[stat.originVertex];
00471 newTrack->SetVertex(vertex);
00472
00473 if(track->GetParentID()==0)
00474 {
00475 stat.trackIndex = m_count;
00476 m_count++;
00477 }
00478 else
00479 stat.trackIndex = vertex->GetMinDau() + vertex->GetCurrentDau()-1;
00480
00481 newTrack->SetIndex( stat.trackIndex );
00482 m_trackIndex = stat.trackIndex;
00483 if(m_logLevel<=4)
00484 G4cout<<" m_trackIndex:"<<m_trackIndex<<G4endl;
00485
00486 newTrack->SetG4TrackId( track->GetTrackID());
00487 m_trackList->push_back( newTrack );
00488
00489 BesTruthTrack* parent = newTrack->GetVertex()->GetParentTrack();
00490 if(parent)
00491 {
00492 parent->AddDaughterIndex(stat.trackIndex);
00493 if(m_logLevel<=4)
00494 G4cout<<"add this daughter to:"<<parent->GetIndex()<<G4endl;
00495 }
00496 }
00497
00498
00499
00500
00501
00502
00503
00504
00505 BesTStats BesSensitiveManager::FollowTrack( const G4Track *track )
00506 {
00507
00508 BesTStats stat( track->GetTrackID(),
00509 HepLorentzVector( track->GetMomentum(), track->GetTotalEnergy() ),
00510 track->GetGlobalTime());
00511
00512
00513 G4int parentId = track->GetParentID();
00514 G4int pdg = track->GetDefinition()->GetPDGEncoding();
00515 G4String particleName = track->GetDefinition()->GetParticleName();
00516 G4ThreeVector trackPos = track->GetPosition();
00517 G4double diffT = m_t0-track->GetGlobalTime();
00518 G4ThreeVector diffPos = (m_pos0-trackPos).mag2();
00519
00520 if (parentId == 0)
00521 {
00522 if(m_logLevel<=4)
00523 {
00524 G4cout<<G4endl;
00525 G4cout<<"trackId:"<<track->GetTrackID()<<" ";
00526 G4cout<<"parentId:"<<parentId<<" ";
00527 G4cout<<"pdg:"<<pdg<<" ";
00528 G4cout<<"name:"<<particleName<<" ";
00529 G4cout<<"timeDiff:"<<diffT<<" ";
00530 G4cout<<"posDiff:"<<diffPos<<G4endl;
00531 }
00532
00533
00534 chain.clear();
00535
00536 stat.savedByDefault = true;
00537
00538 if(m_hepmcevt==0)
00539 MakeNewTrack( stat, track );
00540 else
00541 {
00542 UpdatePrimaryTrack(track);
00543 m_trackFlag = 1;
00544 }
00545 return stat;
00546 }
00547
00548
00549
00550 for(;;)
00551 {
00552 if (chain.back().G4index == parentId) break;
00553 chain.pop_back();
00554 }
00555
00556
00557
00558 if (chain.back().savedByDefault)
00559 {
00560
00561
00562
00563
00564 G4String processName = track->GetCreatorProcess()->GetProcessName();
00565 if (processName=="Decay")
00566 {
00567 if(m_logLevel<=4)
00568 {
00569 G4cout<<G4endl;
00570 G4cout<<"trackId: "<<track->GetTrackID()<<" ";
00571 G4cout<<"parentId: "<<parentId<<" ";
00572 G4cout<<"pdg: "<<pdg<<" ";
00573 G4cout<<"pname: "<<particleName<<G4endl;
00574 }
00575 stat.savedByDefault = true;
00576
00577 if(CheckDecayTrack(track))
00578 m_trackFlag = 1;
00579 else
00580 MakeNewTrack( stat, track);
00581 return stat;
00582 }
00583 }
00584
00585
00586 stat.savedByDefault = false;
00587 return stat;
00588 }
00589
00590 G4bool BesSensitiveManager::CheckDecayTrack(const G4Track* track)
00591 {
00592 G4int flag = 0;
00593 G4int trackID = track->GetTrackID();
00594 G4int parentID = track->GetParentID();
00595 G4int pdgcode = track->GetDefinition()->GetPDGEncoding();
00596 G4ThreeVector p3 = track->GetMomentum();
00597 if(m_logLevel<=4)
00598 G4cout<<"CheckDecayTrack p3: "<<p3<<" "<<track->GetTotalEnergy()<<G4endl;
00599
00600
00601
00602 BesTruthTrack* truthTrack=0;
00603 G4int nTrack = m_trackList->size();
00604 G4int parentTrackIndex=-99;
00605 G4int terminalVertexIndex=-99;
00606 BesTruthTrack* pTrack;
00607 for(int i=0;i<nTrack;i++)
00608 {
00609 truthTrack = (*m_trackList)[i];
00610 if(truthTrack->GetG4TrackId() == parentID)
00611 {
00612
00613 parentTrackIndex = truthTrack->GetIndex();
00614 if(truthTrack->GetTerminalVertex())
00615 {
00616 terminalVertexIndex = truthTrack->GetTerminalVertex()->GetIndex();
00617 pTrack = truthTrack;
00618 if(m_logLevel<=4)
00619 {
00620 G4cout<<"parentTrackIndex:"<<parentTrackIndex<<" "
00621 <<"parent terminate at vertex: "<<terminalVertexIndex<<G4endl;
00622 if(parentTrackIndex != i)
00623 G4cout<<"i: "<<i<<std::endl;
00624 }
00625 break;
00626 }
00627 }
00628 }
00629 if(parentTrackIndex==-99 || terminalVertexIndex==-99)
00630 {
00631 G4cout<<"MatchDecayError!"<<G4endl;
00632
00633 exit(1);
00634 }
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644 G4double minDiffP=1e99;
00645 G4int truthI=-9;
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
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
00687
00688
00689
00690
00691
00692
00693
00694 if(m_logLevel<=4)
00695 std::cout<<"chain.back().vertices.size(): "<<chain.back().vertices.size()<<std::endl;
00696
00697
00698 if(chain.back().vertices.size()<1)
00699 {
00700 if(m_logLevel<=4)
00701 std::cout<<"trackList size: "<<m_trackList->size()<<std::endl;
00702 std::vector<int>* vDau = new std::vector<int>;
00703 GetDaughterVertexes(pTrack, vDau);
00704 G4int sizev = vDau->size();
00705 if(sizev>0)
00706 {
00707 for(int i=0;i<nTrack;i++)
00708 {
00709 truthTrack = (*m_trackList)[i];
00710 G4int vIndex = truthTrack->GetVertex()->GetIndex();
00711 G4int pdgT = truthTrack->GetPDGCode();
00712 if(pdgT==-22) pdgT = 22;
00713 if(pdgT == pdgcode)
00714 {
00715 G4bool matchFlag = MatchVertex(vIndex, vDau);
00716 if(matchFlag)
00717 {
00718 HepLorentzVector tp4 = truthTrack->GetP4();
00719 G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
00720 G4double diffP = (p3-tp3).mag2();
00721 if(m_logLevel<=4)
00722 {
00723 G4cout<<"index: "<<truthTrack->GetIndex()<<"tp3: "<<tp3<<G4endl;
00724 G4cout<<"diffP: "<<diffP<<G4endl;
00725 }
00726 if(diffP<minDiffP)
00727 {
00728 minDiffP = diffP; truthI = i;
00729 if(m_logLevel<=4)
00730 G4cout<<"truthI: "<<i<<G4endl;
00731 }
00732 }
00733 }
00734 }
00735 if(vDau)
00736 delete vDau;
00737 }
00738 }
00739
00740 if(truthI!=-9)
00741 {
00742 m_trackIndex = truthI;
00743 if(m_logLevel<=4)
00744 {
00745 G4cout<<"MatchDecayedTrackWithTrueMother"<<G4endl;
00746 G4cout<<"MatchWithTrueMother m_trackIndex: "<<m_trackIndex<<G4endl;
00747 if(minDiffP>1e-9)
00748 G4cout<<"large minDiffP: "<<minDiffP<<G4endl;
00749 }
00750 truthTrack=(*m_trackList)[truthI];
00751 G4double mass = truthTrack->GetP4().m();
00752 G4double E = sqrt(mass*mass+p3.x()*p3.x()+p3.y()*p3.y()+p3.z()*p3.z());
00753 truthTrack->GetP4().setX(p3.x());
00754 truthTrack->GetP4().setY(p3.y());
00755 truthTrack->GetP4().setZ(p3.z());
00756 truthTrack->GetP4().setE(E);
00757 HepLorentzVector p4 = truthTrack->GetP4();
00758 if(m_logLevel<=4)
00759 {
00760 G4cout<<"primary p: "<<p4.x()<<" "<<p4.y()<<" "<<p4.z()<<" "<<p4.m()<<G4endl;
00761 }
00762 truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
00763 truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
00764 truthTrack->SetG4TrackId(track->GetTrackID());
00765 truthTrack->Found();
00766
00767
00768 return true;
00769 }
00770 return false;
00771 }
00772
00773 void BesSensitiveManager::GetDaughterVertexes(BesTruthTrack* pTrack, std::vector<int>* vDau)
00774 {
00775 G4int size = pTrack->GetDaughterIndexes().size();
00776 if(size>0)
00777 {
00778 G4int minDau = (pTrack->GetDaughterIndexes())[0];
00779 G4int maxDau = (pTrack->GetDaughterIndexes())[size-1];
00780
00781 for(G4int dau=minDau;dau<=maxDau;dau++)
00782 {
00783 if(dau<m_trackList->size())
00784 {
00785 BesTruthTrack* truthTrack = (*m_trackList)[dau];
00786 if(truthTrack->GetVertex())
00787 {
00788 vDau->push_back(truthTrack->GetVertex()->GetIndex());
00789 GetDaughterVertexes(truthTrack,vDau);
00790 }
00791 }
00792 }
00793 }
00794 }
00795
00796 G4bool BesSensitiveManager::MatchVertex(G4int vIndex, std::vector<int>* vDau)
00797 {
00798 G4int size = vDau->size();
00799 if(size>0)
00800 {
00801 for(G4int i=0;i<size;i++)
00802 {
00803 if(vIndex==(*vDau)[i])
00804 return true;
00805 }
00806 }
00807 return false;
00808 }
00809
00810 G4bool BesSensitiveManager::MatchDaughterTrack(const G4Track* track)
00811 {
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855 return true;
00856 }
00857
00858 void BesSensitiveManager::UpdatePrimaryTrack(const G4Track* track)
00859 {
00860 if(m_logLevel<=4)
00861 {
00862 G4cout<<"UpdatePrimaryTrack! ";
00863 G4cout<<"position: "<<track->GetPosition()<<G4endl;
00864 G4cout<<"time: "<<track->GetGlobalTime()<<G4endl;
00865 G4cout<<"momentum: "<<track->GetMomentum()<<" "<<track->GetTotalEnergy()<<G4endl;
00866 }
00867 G4int pdgcode = track->GetDefinition()->GetPDGEncoding();
00868 G4ThreeVector p = track->GetMomentum();
00869
00870 G4int nTrack = m_trackList->size();
00871 BesTruthTrack* truthTrack;
00872 G4int flag = 0;
00873 G4double minDiffP=1e99;
00874 for(int i=0;i<nTrack;i++)
00875 {
00876 truthTrack = (*m_trackList)[i];
00877 HepLorentzVector tp4 = truthTrack->GetP4();
00878 G4ThreeVector tp3(tp4.x(), tp4.y(), tp4.z());
00879 G4double diffP=(p-tp3).mag2();
00880 G4int pdgT = truthTrack->GetPDGCode();
00881 if(pdgT==-22) pdgT=22;
00882 if(pdgcode == pdgT && diffP<minDiffP)
00883 minDiffP = diffP;
00884 if(pdgcode == pdgT && diffP < 1E-9 &&
00885 truthTrack->NotFound())
00886 {
00887 m_trackIndex = truthTrack->GetIndex();
00888 if(m_logLevel<=4)
00889 G4cout<<"m_trackIndex: "<<m_trackIndex<<" "<<G4endl;
00890 truthTrack->SetPDGCharge(track->GetDefinition()->GetPDGCharge());
00891 truthTrack->SetParticleName(track->GetDefinition()->GetParticleName());
00892 truthTrack->SetG4TrackId(track->GetTrackID());
00893 truthTrack->Found();
00894 G4double mass = truthTrack->GetP4().m();
00895 G4double E = sqrt(mass*mass+p.x()*p.x()+p.y()*p.y()+p.z()*p.z());
00896 truthTrack->GetP4().setX(p.x());
00897 truthTrack->GetP4().setY(p.y());
00898 truthTrack->GetP4().setZ(p.z());
00899 truthTrack->GetP4().setT(E);
00900 if(m_logLevel<=4)
00901 {
00902 HepLorentzVector p4 = truthTrack->GetP4();
00903 G4cout<<"primary p: "<<p4.x()<<" "<<p4.y()<<" "<<p4.z()<<" "<<p4.m()<<G4endl;
00904 }
00905 G4int size = truthTrack->GetDaughterIndexes().size();
00906 if(size>0)
00907 {
00908 G4double minDau = (truthTrack->GetDaughterIndexes())[0];
00909 G4double maxDau = (truthTrack->GetDaughterIndexes())[size-1];
00910 if(m_logLevel<=4)
00911 G4cout<<"daughters: "<<minDau<<" "<<maxDau<<G4endl;
00912 }
00913 else
00914 {
00915 if(m_logLevel<=4)
00916 G4cout<<G4endl;
00917 }
00918 flag = 1; break;
00919 }
00920 }
00921 if(flag==0)
00922 {
00923 G4cout<<"MatchError! parentID=0, but match no track in trackList"<<G4endl;
00924 G4cout<<"pdgcode: "<<pdgcode<<" min p diff: "<<minDiffP<<G4endl;
00925 exit(1);
00926 }
00927 }
00928
00929