00001
00009
00010 #include "BesMdcHit.hh"
00011 #include "BesMdcDigi.hh"
00012 #include "BesTofHit.hh"
00013 #include "BesTofDigi.hh"
00014 #include "BesEmcHit.hh"
00015 #include "BesEmcDigi.hh"
00016 #include "BesMucHit.hh"
00017 #include "BesMucDigi.hh"
00018 #include "BesEventAction.hh"
00019 #include "G4RunManager.hh"
00020 #include "BesTruthTrack.hh"
00021 #include "BesTruthVertex.hh"
00022 #include "BesSensitiveManager.hh"
00023 #include "G4HCofThisEvent.hh"
00024 #include "G4SDManager.hh"
00025 #include "G4PrimaryVertex.hh"
00026 #include "G4PrimaryParticle.hh"
00027
00028 #include "BesAsciiIO.hh"
00029 #include "AsciiDmp/dmplib.hh"
00030 #include "G4DigiManager.hh"
00031 #include <iostream>
00032 using namespace std;
00033
00034 #include "GaudiKernel/IDataProviderSvc.h"
00035 #include "GaudiKernel/ISvcLocator.h"
00036 #include "GaudiKernel/Bootstrap.h"
00037 #include "GaudiKernel/RegistryEntry.h"
00038 #include "GaudiKernel/MsgStream.h"
00039 #include "GaudiKernel/SmartDataPtr.h"
00040 #include "McTruth/DecayMode.h"
00041
00042 BesAsciiIO::BesAsciiIO(
00043 G4int mdcTruFlag,G4int mdcDigiFlag,
00044 G4int tofTruFlag, G4int tofDigiFlag,
00045 G4int emcTruFlag, G4int emcDigiFlag,
00046 G4int mucTruFlag, G4int mucDigiFlag,
00047 G4String name)
00048 :m_mdcTruFlag(mdcTruFlag),m_mdcDigiFlag(mdcDigiFlag),
00049 m_tofTruFlag(tofTruFlag),m_tofDigiFlag(tofDigiFlag),
00050 m_emcTruFlag(emcTruFlag),m_emcDigiFlag(emcDigiFlag),
00051 m_mucTruFlag(mucTruFlag),m_mucDigiFlag(mucDigiFlag),
00052 m_asciiFile(name)
00053 {
00054 m_DigiMan = G4DigiManager::GetDMpointer();
00055 }
00056
00057 BesAsciiIO::~BesAsciiIO()
00058 {
00059 }
00060
00061 void BesAsciiIO::SaveAsciiEvents(G4int runId, const G4Event* evt)
00062 {
00063 EVENT asciiEvt;
00064
00065 asciiEvt.set_initialized();
00066 asciiEvt.header.set_initialized();
00067 asciiEvt.header.eventNo=evt->GetEventID();
00068 asciiEvt.header.runNo= runId;
00069 asciiEvt.decayMode.set_initialized();
00070 SaveDecayMode(asciiEvt);
00071
00072 asciiEvt.trackTruth.set_initialized();
00073 asciiEvt.vertexTruth.set_initialized();
00074 SaveTrackTruth(asciiEvt);
00075 SaveVertexTruth(asciiEvt);
00076
00077 if(m_mdcTruFlag)
00078 {
00079 asciiEvt.mdcTruth.set_initialized();
00080 SaveMdcTruth(asciiEvt);
00081 }
00082
00083 if(m_mdcDigiFlag)
00084 {
00085 asciiEvt.mdcDigi.set_initialized();
00086 SaveMdcDigits(asciiEvt);
00087 }
00088
00089 if(m_tofTruFlag)
00090 {
00091 asciiEvt.tofTruth.set_initialized();
00092 SaveTofTruth(asciiEvt);
00093 }
00094
00095 if(m_tofDigiFlag)
00096 {
00097 asciiEvt.tofDigi.set_initialized();
00098 SaveTofDigits(asciiEvt);
00099 }
00100
00101 if(m_emcTruFlag)
00102 {
00103 asciiEvt.emcTruth.set_initialized();
00104 SaveEmcTruth(asciiEvt);
00105 }
00106
00107 if(m_emcDigiFlag)
00108 {
00109 asciiEvt.emcDigi.set_initialized();
00110 SaveEmcDigits(asciiEvt);
00111 }
00112
00113 if(m_mucTruFlag)
00114 {
00115 asciiEvt.mucTruth.set_initialized();
00116 SaveMucTruth(asciiEvt);
00117 }
00118
00119 if(m_mucDigiFlag)
00120 {
00121 asciiEvt.mucDigi.set_initialized();
00122 SaveMucDigits(asciiEvt);
00123 }
00124
00125 ofstream os;
00126 if(evt->GetEventID()==0)
00127 {
00128 os.open(m_asciiFile);
00129 FRMTVERSION version;
00130 version.set_initialized();
00131 version.major = 1;
00132 version.minor = 0;
00133 os << version;
00134 }
00135 else
00136 os.open(m_asciiFile,ios::out|ios::app);
00137
00138 try {
00139 os << asciiEvt;
00140 }
00141 catch (AsciiWrongTag& ex) {
00142 std::cerr << "wrong tag, got " << ex.got()
00143 << " expected: " << ex.expected()
00144 << std::endl;
00145 } catch (AsciiDumpException& ) {
00146 std::cerr << "AsciiDumpException was caught!" << std::endl;
00147 }
00148 os.close();
00149 }
00150
00151 void BesAsciiIO::SaveDecayMode(EVENT& asciiEvt)
00152 {
00153 asciiEvt.decayMode.size=10;
00154
00155 ISvcLocator* svcLocator = Gaudi::svcLocator();
00156 IDataProviderSvc* evtSvc;
00157 StatusCode sc=svcLocator->service("EventDataSvc", evtSvc);
00158 if (sc.isFailure())
00159 G4cout<<"Could not accesss EventDataSvc!"<<G4endl;
00160
00161 SmartDataPtr<DecayMode> decayMode(evtSvc,"/Event/MC/DecayMode");
00162 if(!decayMode)
00163 {
00164 for(int i=0;i<10;i++)
00165 asciiEvt.decayMode.data[i]=0;
00166 }
00167 else
00168 {
00169 int dm[10]={0,0,0,0,0,0,0,0,0,0};
00170 decayMode->getData(dm,10);
00171 for(int i=0;i<10;i++)
00172 asciiEvt.decayMode.data[i]=dm[i];
00173 }
00174 }
00175
00176 void BesAsciiIO::SaveTrackTruth(EVENT& asciiEvt)
00177 {
00178
00179 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
00180
00181 vector<BesTruthTrack*>* trackList = sensitiveManager->GetTrackList();
00182
00183
00184 G4int nTrack = trackList->size();
00185 BesTruthTrack* track;
00186 for(int i=0;i<nTrack-1;i++)
00187 for(int j=i+1;j<nTrack;j++)
00188 if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
00189 {
00190 track=(*trackList)[i];
00191 (*trackList)[i]=(*trackList)[j];
00192 (*trackList)[j]=track;
00193 }
00194
00195 for(int i=0;i<nTrack;i++)
00196 {
00197 TrackTruthType trackTruth;
00198 track = (*(sensitiveManager->GetTrackList())) [i];
00199 trackTruth.trackIndex = track->GetIndex();
00200 if(track->GetIndex()<0)
00201 G4cout<<"trackIndex<0!"<<G4endl;
00202 trackTruth.PDGCode = track->GetPDGCode();
00203 trackTruth.PDGCharge = track->GetPDGCharge();
00204 trackTruth.v0Index = track->GetVertex()->GetIndex();
00205 if( track->GetTerminalVertex() )
00206 trackTruth.v1Index = track->GetTerminalVertex()->GetIndex();
00207 else
00208 trackTruth.v1Index = -99;
00209 trackTruth.px = track->GetP4().x()/1000.;
00210 trackTruth.py = track->GetP4().y()/1000.;
00211 trackTruth.pz = track->GetP4().z()/1000.;
00212 trackTruth.E = track->GetP4().t()/1000.;
00213
00214 G4int size = track->GetDaughterIndexes().size();
00215 if(size>0)
00216 {
00217 trackTruth.minDaughterIndex = (track->GetDaughterIndexes())[0];
00218 trackTruth.maxDaughterIndex = (track->GetDaughterIndexes())[size-1];
00219 }
00220 else
00221 {
00222 trackTruth.minDaughterIndex = -99;
00223 trackTruth.maxDaughterIndex = -99;
00224 }
00225
00226
00227 asciiEvt.trackTruth.truthCol.push_back(trackTruth);
00228 }
00229
00230 asciiEvt.trackTruth.nTruth=asciiEvt.trackTruth.truthCol.size();
00231 }
00232
00233 void BesAsciiIO::SaveVertexTruth(EVENT& asciiEvt)
00234 {
00235 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
00236
00237 vector<BesTruthVertex*>* vertexList = sensitiveManager->GetVertexList();
00238
00239 G4int nVertex = vertexList->size();
00240 BesTruthVertex* vertex;
00241 for(int i=0;i<nVertex;i++)
00242 {
00243 VertexTruthType vertexTruth;
00244 vertex = (*vertexList) [i];
00245 vertexTruth.vertexIndex = vertex->GetIndex();
00246 if(vertex->GetParentTrack())
00247 vertexTruth.parentTrackIndex = vertex->GetParentTrack()->GetIndex();
00248 else
00249 vertexTruth.parentTrackIndex = -99;
00250
00251 vertexTruth.posX = vertex->GetPosition().x()/10.;
00252 vertexTruth.posY = vertex->GetPosition().y()/10.;
00253 vertexTruth.posZ = vertex->GetPosition().z()/10.;
00254 vertexTruth.time = vertex->GetTime();
00255 asciiEvt.vertexTruth.truthCol.push_back(vertexTruth);
00256 }
00257 asciiEvt.vertexTruth.nTruth = asciiEvt.vertexTruth.truthCol.size();
00258 }
00259
00260 void BesAsciiIO::SaveMdcTruth(EVENT& asciiEvt)
00261 {
00262 G4int HCID = -1;
00263 HCID = m_DigiMan->GetHitsCollectionID("BesMdcTruthCollection");
00264 if(HCID>0)
00265 {
00266 BesMdcHitsCollection* HC = 0;
00267 HC = (BesMdcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00268 G4int n_hit = HC->entries();
00269 if(n_hit>0)
00270 {
00271
00272 BesMdcHit* hit;
00273 vector<BesMdcHit*>* vecHC = HC->GetVector();
00274 for(int i=0;i<n_hit-1;i++)
00275 for(int j=i+1;j<n_hit;j++)
00276 if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
00277 {
00278 hit = (*vecHC)[i];
00279 (*vecHC)[i] = (*vecHC)[j];
00280 (*vecHC)[j] = hit;
00281 }
00282
00283 for(G4int i=0;i<n_hit;i++)
00284 {
00285 hit = (*HC)[i];
00286 MdcTruthType mdcTruth;
00287 mdcTruth.trackIndex = hit->GetTrackID();
00288 mdcTruth.layerNo = hit->GetLayerNo();
00289 mdcTruth.cellNo = hit->GetCellNo();
00290 mdcTruth.edep = hit->GetEdep();
00291 mdcTruth.driftD = hit->GetDriftD();
00292 mdcTruth.posX = hit->GetPos().x();
00293 mdcTruth.posY = hit->GetPos().y();
00294 mdcTruth.posZ = hit->GetPos().z();
00295 mdcTruth.posFlag = hit->GetPosFlag();
00296
00297 asciiEvt.mdcTruth.truthCol.push_back(mdcTruth);
00298 }
00299 }
00300 }
00301 asciiEvt.mdcTruth.nTruth = asciiEvt.mdcTruth.truthCol.size();
00302 }
00303
00304 void BesAsciiIO::SaveTofTruth(EVENT& asciiEvt)
00305 {
00306 G4int HCID = -1;
00307 HCID = m_DigiMan->GetHitsCollectionID("BesTofHitsList");
00308 if(HCID>0)
00309 {
00310 BesTofHitsCollection* HC = 0;
00311 HC = (BesTofHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00312 G4int n_hit = HC->entries();
00313 if(n_hit>0)
00314 {
00315
00316 BesTofHit* hit;
00317 vector<BesTofHit*>* vecHC = HC->GetVector();
00318 for(int i=0;i<n_hit-1;i++)
00319 for(int j=i+1;j<n_hit;j++)
00320 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
00321 {
00322 hit = (*vecHC)[i];
00323 (*vecHC)[i] = (*vecHC)[j];
00324 (*vecHC)[j] = hit;
00325 }
00326
00327 for(G4int i=0;i<n_hit;i++)
00328 {
00329 hit = (*HC)[i];
00330 TofTruthType tofTruth;
00331 tofTruth.trackIndex = hit->GetTrackIndex();
00332 tofTruth.partId = hit->GetPartId();
00333 tofTruth.scinNb = hit->GetScinNb();
00334 tofTruth.posX = hit->GetPos().x();
00335 tofTruth.posY = hit->GetPos().y();
00336 tofTruth.posZ = hit->GetPos().z();
00337 tofTruth.px = hit->GetMomentum().x();
00338 tofTruth.py = hit->GetMomentum().y();
00339 tofTruth.pz = hit->GetMomentum().z();
00340 tofTruth.trackL = hit->GetTrackL();
00341 tofTruth.time = hit->GetTime();
00342
00343 asciiEvt.tofTruth.truthCol.push_back(tofTruth);
00344 }
00345 }
00346 }
00347 asciiEvt.tofTruth.nTruth=asciiEvt.tofTruth.truthCol.size();
00348 }
00349
00350 void BesAsciiIO::SaveEmcTruth(EVENT& asciiEvt)
00351 {
00352 G4int HCID = -1;
00353 HCID = m_DigiMan->GetHitsCollectionID("BesEmcHitsList");
00354 if(HCID>0)
00355 {
00356 BesEmcHitsCollection* HC = 0;
00357 HC = (BesEmcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00358 G4int n_hit = HC->entries();
00359 if(n_hit>0)
00360 {
00361
00362 BesEmcHit* hit;
00363 vector<BesEmcHit*>* vecHC = HC->GetVector();
00364 for(int i=0;i<n_hit-1;i++)
00365 for(int j=i+1;j<n_hit;j++)
00366 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
00367 {
00368 hit = (*vecHC)[i];
00369 (*vecHC)[i] = (*vecHC)[j];
00370 (*vecHC)[j] = hit;
00371 }
00372
00373 for(G4int i=0;i<n_hit;i++)
00374 {
00375 hit = (*HC)[i];
00376 EmcTruthType emcTruth;
00377 emcTruth.trackIndex = hit->GetTrackIndex();
00378 emcTruth.partId = hit->GetPartId();
00379 emcTruth.numTheta = hit->GetNumThetaCrystal();
00380 emcTruth.numPhi = hit->GetNumPhiCrystal();
00381 emcTruth.posX = hit->GetPosCrystal().x();
00382 emcTruth.posY = hit->GetPosCrystal().y();
00383 emcTruth.posZ = hit->GetPosCrystal().z();
00384 emcTruth.px = hit->GetMomentum().x();
00385 emcTruth.py = hit->GetMomentum().y();
00386 emcTruth.pz = hit->GetMomentum().z();
00387 emcTruth.totalEdep = hit->GetEdepCrystal();
00388
00389 asciiEvt.emcTruth.truthCol.push_back(emcTruth);
00390 }
00391 }
00392 }
00393 asciiEvt.emcTruth.nTruth=asciiEvt.emcTruth.truthCol.size();
00394 }
00395
00396 void BesAsciiIO::SaveMucTruth(EVENT& asciiEvt)
00397 {
00398 G4int HCID = -1;
00399 HCID = m_DigiMan->GetHitsCollectionID("BesMucHitsList");
00400 if(HCID>0)
00401 {
00402 BesMucHitsCollection* HC = 0;
00403 HC = (BesMucHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00404 G4int n_hit = HC->entries();
00405 if(n_hit>0)
00406 {
00407
00408 BesMucHit* hit;
00409 vector<BesMucHit*>* vecHC = HC->GetVector();
00410 for(int i=0;i<n_hit-1;i++)
00411 for(int j=i+1;j<n_hit;j++)
00412 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
00413 {
00414 hit = (*vecHC)[i];
00415 (*vecHC)[i] = (*vecHC)[j];
00416 (*vecHC)[j] = hit;
00417 }
00418
00419 for(G4int i=0;i<n_hit;i++)
00420 {
00421 hit = (*HC)[i];
00422 MucTruthType mucTruth;
00423 mucTruth.trackIndex = hit->GetTrackIndex();
00424 mucTruth.partId = hit->GetPart();
00425 mucTruth.segId = hit->GetSeg();
00426 mucTruth.gapId = hit->GetGap();
00427 mucTruth.stripId = hit->GetStrip();
00428 mucTruth.posX = hit->GetPos().x();
00429 mucTruth.posY = hit->GetPos().y();
00430 mucTruth.posZ = hit->GetPos().z();
00431 mucTruth.px = hit->GetMomentum().x();
00432 mucTruth.py = hit->GetMomentum().y();
00433 mucTruth.pz = hit->GetMomentum().z();
00434 asciiEvt.mucTruth.truthCol.push_back(mucTruth);
00435 }
00436 }
00437 }
00438 asciiEvt.mucTruth.nTruth=asciiEvt.mucTruth.truthCol.size();
00439 }
00440
00441 void BesAsciiIO::SaveMdcDigits(EVENT& asciiEvt)
00442 {
00443 G4int mdcDigiCollID = -1;
00444 mdcDigiCollID = m_DigiMan->GetDigiCollectionID("BesMdcDigisCollection");
00445 if(mdcDigiCollID>=0){
00446 BesMdcDigisCollection* mdcDC = (BesMdcDigisCollection*)m_DigiMan->GetDigiCollection(mdcDigiCollID);
00447 G4int nDigi = mdcDC->entries();
00448 if(nDigi>0)
00449 {
00450 BesMdcDigi* mdcDigi;
00451 for(int i=0;i<nDigi;i++)
00452 {
00453 MdcDigiType mdcData;
00454 mdcDigi=(*mdcDC)[i];
00455 mdcData.trackIndex = mdcDigi->GetTrackID();
00456 mdcData.layerNo = mdcDigi->GetLayerNo();
00457 mdcData.cellNo = mdcDigi->GetCellNo();
00458 mdcData.energyDeposit = mdcDigi->GetEdep();
00459 mdcData.driftTime = mdcDigi->GetDriftT();
00460 asciiEvt.mdcDigi.digiCol.push_back(mdcData);
00461 }
00462 }
00463 }
00464 asciiEvt.mdcDigi.nDigi=asciiEvt.mdcDigi.digiCol.size();
00465 }
00466
00467
00468 void BesAsciiIO::SaveTofDigits(EVENT& asciiEvt)
00469 {
00470
00471 G4int tofDigiCollID = -1;
00472 tofDigiCollID = m_DigiMan->GetDigiCollectionID("BesTofDigitsCollection");
00473
00474 if(tofDigiCollID>=0)
00475 {
00476 BesTofDigitsCollection* tofDC = (BesTofDigitsCollection*)m_DigiMan->GetDigiCollection(tofDigiCollID);
00477 G4int nDigi = tofDC->entries();
00478 if(nDigi>0)
00479 {
00480
00481 BesTofDigi* digi;
00482 vector<BesTofDigi*>* vecDC = tofDC->GetVector();
00483 for(int i=0;i<nDigi-1;i++)
00484 for(int j=i+1;j<nDigi;j++)
00485 if((*vecDC)[i]->GetTrackIndex()>(*vecDC)[j]->GetTrackIndex())
00486 {
00487 digi = (*vecDC)[i];
00488 (*vecDC)[i] = (*vecDC)[j];
00489 (*vecDC)[j] = digi;
00490 }
00491
00492 for(int i=0;i<nDigi;i++)
00493 {
00494 TofDigiType tofData;
00495 digi = (*tofDC)[i];
00496 tofData.trackIndex = digi->GetTrackIndex();
00497 tofData.partId = digi->GetPartId();
00498 tofData.scinNb = digi->GetScinNb();
00499 tofData.forwADC = digi->GetForwADC();
00500 tofData.forwTDC = digi->GetForwTDC();
00501 tofData.backADC = digi->GetBackADC();
00502 tofData.backTDC = digi->GetBackTDC();
00503
00504 asciiEvt.tofDigi.digiCol.push_back(tofData);
00505 }
00506 }
00507 }
00508 asciiEvt.tofDigi.nDigi=asciiEvt.tofDigi.digiCol.size();
00509 }
00510
00511 void BesAsciiIO::SaveEmcDigits(EVENT& asciiEvt)
00512 {
00513 G4int emcDigiCollID = -1;
00514 emcDigiCollID = m_DigiMan->GetDigiCollectionID("BesEmcDigitsCollection");
00515 if(emcDigiCollID>=0)
00516 {
00517 BesEmcDigitsCollection* emcDC = (BesEmcDigitsCollection*)m_DigiMan->GetDigiCollection(emcDigiCollID);
00518 G4int nDigi = emcDC->entries();
00519 if(nDigi>0)
00520 {
00521
00522 BesEmcDigi* digi;
00523 vector<BesEmcDigi*>* vecDC = emcDC->GetVector();
00524 for(int i=0;i<nDigi-1;i++)
00525 for(int j=i+1;j<nDigi;j++)
00526 if((*vecDC)[i]->GetTrackIndex()>(*vecDC)[j]->GetTrackIndex())
00527 {
00528 digi = (*vecDC)[i];
00529 (*vecDC)[i] = (*vecDC)[j];
00530 (*vecDC)[j] = digi;
00531 }
00532 for(int i=0;i<nDigi;i++)
00533 {
00534 EmcDigiType emcData;
00535 digi = (*emcDC)[i];
00536 emcData.trackIndex = digi->GetTrackIndex();
00537 emcData.partId = digi->GetPartId();
00538 emcData.numTheta = digi->GetThetaNb();
00539 emcData.numPhi = digi->GetPhiNb();
00540 emcData.energyDeposit = digi->GetEnergy();
00541 emcData.hitTime = (G4double)digi->GetTime();
00542 asciiEvt.emcDigi.digiCol.push_back(emcData);
00543 }
00544 }
00545 }
00546 asciiEvt.emcDigi.nDigi=asciiEvt.emcDigi.digiCol.size();
00547 }
00548
00549 void BesAsciiIO::SaveMucDigits(EVENT& asciiEvt)
00550 {
00551 G4int mucDigiCollID =-1;
00552 mucDigiCollID = m_DigiMan->GetDigiCollectionID("BesMucDigisCollection");
00553 if(mucDigiCollID>=0)
00554 {
00555 BesMucDigisCollection* mucDC = (BesMucDigisCollection*)m_DigiMan->GetDigiCollection(mucDigiCollID);
00556 G4int nDigi = mucDC->entries();
00557 if(nDigi > 0) {
00558 BesMucDigi* mucDigi;
00559 for(int i = 0; i < nDigi; i++)
00560 {
00561 MucDigiType mucData;
00562 mucDigi = (*mucDC)[i];
00563 mucData.trackIndex = mucDigi->GetTrackIndex();
00564 mucData.partNo = mucDigi->GetPartId();
00565 mucData.segNo = mucDigi->GetSegId();
00566 mucData.gapNo = mucDigi->GetGapId();
00567 mucData.stripNo = mucDigi->GetStripId();
00568
00569 asciiEvt.mucDigi.digiCol.push_back(mucData);
00570 }
00571 }
00572 }
00573 asciiEvt.mucDigi.nDigi=asciiEvt.mucDigi.digiCol.size();
00574 }
00575
00576
00577 void BesAsciiIO::SaveHitAsciiEvents(G4int runId, const G4Event* evt){
00578 HitEVENT asciiEvt;
00579
00580 asciiEvt.set_initialized();
00581 asciiEvt.header.set_initialized();
00582 asciiEvt.header.eventNo=evt->GetEventID();
00583 asciiEvt.header.runNo= runId;
00584 asciiEvt.decayMode.set_initialized();
00585 SaveDecayMode(asciiEvt);
00586
00587 asciiEvt.trackTruth.set_initialized();
00588 asciiEvt.vertexTruth.set_initialized();
00589 SaveTrackTruth(asciiEvt);
00590 SaveVertexTruth(asciiEvt);
00591
00592 if(m_mdcTruFlag)
00593 {
00594 asciiEvt.mdcTruth.set_initialized();
00595 SaveMdcTruth(asciiEvt);
00596 }
00597
00598 if(m_mdcDigiFlag)
00599 {
00600 asciiEvt.mdcHit.set_initialized();
00601 SaveMdcHits(asciiEvt);
00602 }
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640 ofstream os;
00641 if(evt->GetEventID()==0){
00642 os.open(m_asciiFile);
00643 FRMTVERSION version;
00644 version.set_initialized();
00645 version.major = 1;
00646 version.minor = 0;
00647 os << version;
00648 }
00649 else
00650 os.open(m_asciiFile,ios::out|ios::app);
00651
00652 try {
00653 os << asciiEvt;
00654 }catch (AsciiWrongTag& ex) {
00655 std::cerr << "wrong tag, got " << ex.got()
00656 << " expected: " << ex.expected()
00657 << std::endl;
00658 } catch (AsciiDumpException& ) {
00659 std::cerr << "AsciiDumpException was caught!" << std::endl;
00660 }
00661
00662 os.close();
00663 }
00664
00665 void BesAsciiIO::SaveDecayMode(HitEVENT& asciiEvt)
00666 {
00667 asciiEvt.decayMode.size=10;
00668
00669 ISvcLocator* svcLocator = Gaudi::svcLocator();
00670 IDataProviderSvc* evtSvc;
00671 StatusCode sc=svcLocator->service("EventDataSvc", evtSvc);
00672 if (sc.isFailure())
00673 G4cout<<"Could not accesss EventDataSvc!"<<G4endl;
00674
00675 SmartDataPtr<DecayMode> decayMode(evtSvc,"/Event/MC/DecayMode");
00676 if(!decayMode)
00677 {
00678 for(int i=0;i<10;i++)
00679 asciiEvt.decayMode.data[i]=0;
00680 }
00681 else
00682 {
00683 int dm[10]={0,0,0,0,0,0,0,0,0,0};
00684 decayMode->getData(dm,10);
00685 for(int i=0;i<10;i++)
00686 asciiEvt.decayMode.data[i]=dm[i];
00687 }
00688 }
00689
00690 void BesAsciiIO::SaveTrackTruth(HitEVENT& asciiEvt)
00691 {
00692
00693 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
00694
00695 vector<BesTruthTrack*>* trackList = sensitiveManager->GetTrackList();
00696
00697
00698 G4int nTrack = trackList->size();
00699 BesTruthTrack* track;
00700 for(int i=0;i<nTrack-1;i++)
00701 for(int j=i+1;j<nTrack;j++)
00702 if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
00703 {
00704 track=(*trackList)[i];
00705 (*trackList)[i]=(*trackList)[j];
00706 (*trackList)[j]=track;
00707 }
00708
00709 for(int i=0;i<nTrack;i++)
00710 {
00711 TrackTruthType trackTruth;
00712 track = (*(sensitiveManager->GetTrackList())) [i];
00713 trackTruth.trackIndex = track->GetIndex();
00714 trackTruth.PDGCode = track->GetPDGCode();
00715 trackTruth.PDGCharge = track->GetPDGCharge();
00716 trackTruth.v0Index = track->GetVertex()->GetIndex();
00717 if( track->GetTerminalVertex() )
00718 trackTruth.v1Index = track->GetTerminalVertex()->GetIndex();
00719 else
00720 trackTruth.v1Index = -99;
00721 trackTruth.px = track->GetP4().x();
00722 trackTruth.py = track->GetP4().y();
00723 trackTruth.pz = track->GetP4().z();
00724 trackTruth.E = track->GetP4().t();
00725
00726 G4int size = track->GetDaughterIndexes().size();
00727 if(size>0)
00728 {
00729 trackTruth.minDaughterIndex = (track->GetDaughterIndexes())[0];
00730 trackTruth.maxDaughterIndex = (track->GetDaughterIndexes())[size-1];
00731 }
00732 else
00733 {
00734 trackTruth.minDaughterIndex = -99;
00735 trackTruth.maxDaughterIndex = -99;
00736 }
00737
00738
00739 asciiEvt.trackTruth.truthCol.push_back(trackTruth);
00740 }
00741
00742 asciiEvt.trackTruth.nTruth=asciiEvt.trackTruth.truthCol.size();
00743 }
00744
00745 void BesAsciiIO::SaveVertexTruth(HitEVENT& asciiEvt)
00746 {
00747 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
00748
00749 vector<BesTruthVertex*>* vertexList = sensitiveManager->GetVertexList();
00750
00751 G4int nVertex = vertexList->size();
00752 BesTruthVertex* vertex;
00753 for(int i=0;i<nVertex;i++)
00754 {
00755 VertexTruthType vertexTruth;
00756 vertex = (*vertexList) [i];
00757 vertexTruth.vertexIndex = vertex->GetIndex();
00758 if(vertex->GetParentTrack())
00759 vertexTruth.parentTrackIndex = vertex->GetParentTrack()->GetIndex();
00760 else
00761 vertexTruth.parentTrackIndex = -99;
00762
00763 vertexTruth.posX = vertex->GetPosition().x();
00764 vertexTruth.posY = vertex->GetPosition().y();
00765 vertexTruth.posZ = vertex->GetPosition().z();
00766 vertexTruth.time = vertex->GetTime();
00767 asciiEvt.vertexTruth.truthCol.push_back(vertexTruth);
00768 }
00769 asciiEvt.vertexTruth.nTruth = asciiEvt.vertexTruth.truthCol.size();
00770 }
00771
00772 void BesAsciiIO::SaveMdcTruth(HitEVENT& asciiEvt)
00773 {
00774 G4int HCID = -1;
00775 HCID = m_DigiMan->GetHitsCollectionID("BesMdcTruthCollection");
00776 if(HCID>0)
00777 {
00778 BesMdcHitsCollection* HC = 0;
00779 HC = (BesMdcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
00780 G4int n_hit = HC->entries();
00781 if(n_hit>0)
00782 {
00783
00784 BesMdcHit* hit;
00785 vector<BesMdcHit*>* vecHC = HC->GetVector();
00786 for(int i=0;i<n_hit-1;i++)
00787 for(int j=i+1;j<n_hit;j++)
00788 if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
00789 {
00790 hit = (*vecHC)[i];
00791 (*vecHC)[i] = (*vecHC)[j];
00792 (*vecHC)[j] = hit;
00793 }
00794
00795 for(G4int i=0;i<n_hit;i++)
00796 {
00797 hit = (*HC)[i];
00798 MdcTruthType mdcTruth;
00799 mdcTruth.trackIndex = hit->GetTrackID();
00800 mdcTruth.layerNo = hit->GetLayerNo();
00801 mdcTruth.cellNo = hit->GetCellNo();
00802 mdcTruth.edep = hit->GetEdep();
00803 mdcTruth.driftD = hit->GetDriftD();
00804 mdcTruth.posX = hit->GetPos().x();
00805 mdcTruth.posY = hit->GetPos().y();
00806 mdcTruth.posZ = hit->GetPos().z();
00807 mdcTruth.posFlag = hit->GetPosFlag();
00808
00809 asciiEvt.mdcTruth.truthCol.push_back(mdcTruth);
00810 }
00811 }
00812 }
00813 asciiEvt.mdcTruth.nTruth = asciiEvt.mdcTruth.truthCol.size();
00814 }
00815
00816 void BesAsciiIO::SaveMdcHits(HitEVENT& asciiEvt)
00817 {
00818 G4int mdcHitCollID = -1;
00819 mdcHitCollID = m_DigiMan->GetHitsCollectionID("BesMdcHitsCollection");
00820 if(mdcHitCollID>=0)
00821 {
00822 BesMdcHitsCollection* mdcDC = (BesMdcHitsCollection*)m_DigiMan->GetHitsCollection(mdcHitCollID);
00823 G4int nHit = mdcDC->entries();
00824 if(nHit>0)
00825 {
00826 BesMdcHit* mdcHit;
00827 for(int i=0;i<nHit;i++)
00828 {
00829 MdcHitType mdcData;
00830 mdcHit=(*mdcDC)[i];
00831 mdcData.trackIndex = mdcHit->GetTrackID();
00832 mdcData.layerNo = mdcHit->GetLayerNo();
00833 mdcData.cellNo = mdcHit->GetCellNo();
00834 mdcData.posX = mdcHit->GetPos().x();
00835 mdcData.posY = mdcHit->GetPos().y();
00836 mdcData.posZ = mdcHit->GetPos().z();
00837 mdcData.energyDeposit = mdcHit->GetEdep();
00838 mdcData.driftDistance = mdcHit->GetDriftD();
00839 mdcData.globalT = mdcHit->GetGlobalT();
00840 mdcData.theta = mdcHit->GetTheta();
00841 mdcData.enterAngle = mdcHit->GetEnterAngle();
00842 mdcData.posFlag = mdcHit->GetPosFlag();
00843
00844 asciiEvt.mdcHit.hitCol.push_back(mdcData);
00845 }
00846 }
00847 }
00848 asciiEvt.mdcHit.nHit=asciiEvt.mdcHit.hitCol.size();
00849 }
00850