CGEM BOSS 6.6.5.f
BESIII Offline Software System
Loading...
Searching...
No Matches
bak_BesSim-00-04-14/src/BesAsciiIO.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 "BesMdcHit.hh"
11#include "BesMdcDigi.hh"
12#include "BesCgemHit.hh"
13#include "BesCgemDigi.hh"
14#include "BesTofHit.hh"
15#include "BesTofDigi.hh"
16#include "BesEmcHit.hh"
17#include "BesEmcDigi.hh"
18#include "BesMucHit.hh"
19#include "BesMucDigi.hh"
20#include "BesEventAction.hh"
21#include "G4RunManager.hh"
22#include "BesTruthTrack.hh"
23#include "BesTruthVertex.hh"
24#include "BesSensitiveManager.hh"
25#include "G4HCofThisEvent.hh"
26#include "G4SDManager.hh"
27#include "G4PrimaryVertex.hh"
28#include "G4PrimaryParticle.hh"
29
30#include "BesAsciiIO.hh"
31#include "AsciiDmp/dmplib.hh"
32#include "G4DigiManager.hh"
33#include <iostream>
34using namespace std;
35
36#include "GaudiKernel/IDataProviderSvc.h"
37#include "GaudiKernel/ISvcLocator.h"
38#include "GaudiKernel/Bootstrap.h"
39#include "GaudiKernel/RegistryEntry.h"
40#include "GaudiKernel/MsgStream.h"
41#include "GaudiKernel/SmartDataPtr.h"
42#include "McTruth/DecayMode.h"
43
44BesAsciiIO::BesAsciiIO(
45 G4int mdcTruFlag,G4int mdcDigiFlag,
46 G4int cgemTruFlag,G4int cgemDigiFlag,
47 G4int tofTruFlag, G4int tofDigiFlag,
48 G4int emcTruFlag, G4int emcDigiFlag,
49 G4int mucTruFlag, G4int mucDigiFlag,
50 G4String name)
51:m_mdcTruFlag(mdcTruFlag),m_mdcDigiFlag(mdcDigiFlag),
52 m_cgemTruFlag(cgemTruFlag),m_cgemDigiFlag(cgemDigiFlag),
53 m_tofTruFlag(tofTruFlag),m_tofDigiFlag(tofDigiFlag),
54 m_emcTruFlag(emcTruFlag),m_emcDigiFlag(emcDigiFlag),
55 m_mucTruFlag(mucTruFlag),m_mucDigiFlag(mucDigiFlag),
56 m_asciiFile(name)
57{
58 m_DigiMan = G4DigiManager::GetDMpointer();
59}
60
61BesAsciiIO::~BesAsciiIO()
62{
63}
64
65void BesAsciiIO::SaveAsciiEvents(G4int runId, const G4Event* evt)
66{
67 EVENT asciiEvt;
68
69 asciiEvt.set_initialized();
70 asciiEvt.header.set_initialized();
71 asciiEvt.header.eventNo=evt->GetEventID();
72 asciiEvt.header.runNo= runId;
73 asciiEvt.decayMode.set_initialized();
74 SaveDecayMode(asciiEvt);
75
76 asciiEvt.trackTruth.set_initialized();
77 asciiEvt.vertexTruth.set_initialized();
78 SaveTrackTruth(asciiEvt);
79 SaveVertexTruth(asciiEvt);
80
81 if(m_mdcTruFlag)
82 {
83 asciiEvt.mdcTruth.set_initialized();
84 SaveMdcTruth(asciiEvt);
85 }
86
87 if(m_mdcDigiFlag)
88 {
89 asciiEvt.mdcDigi.set_initialized();
90 SaveMdcDigits(asciiEvt);
91 }
92
93 if(m_cgemTruFlag)
94 {
95 asciiEvt.cgemTruth.set_initialized();
96 SaveCgemTruth(asciiEvt);
97 }
98
99 if(m_cgemDigiFlag)
100 {
101 asciiEvt.cgemDigi.set_initialized();
102 SaveCgemDigits(asciiEvt);
103 }
104
105 if(m_tofTruFlag)
106 {
107 asciiEvt.tofTruth.set_initialized();
108 SaveTofTruth(asciiEvt);
109 }
110
111 if(m_tofDigiFlag)
112 {
113 asciiEvt.tofDigi.set_initialized();
114 SaveTofDigits(asciiEvt);
115 }
116
117 if(m_emcTruFlag)
118 {
119 asciiEvt.emcTruth.set_initialized();
120 SaveEmcTruth(asciiEvt);
121 }
122
123 if(m_emcDigiFlag)
124 {
125 asciiEvt.emcDigi.set_initialized();
126 SaveEmcDigits(asciiEvt);
127 }
128
129 if(m_mucTruFlag)
130 {
131 asciiEvt.mucTruth.set_initialized();
132 SaveMucTruth(asciiEvt);
133 }
134
135 if(m_mucDigiFlag)
136 {
137 asciiEvt.mucDigi.set_initialized();
138 SaveMucDigits(asciiEvt);
139 }
140
141 ofstream os;
142 if(evt->GetEventID()==0)
143 {
144 os.open(m_asciiFile);
145 FRMTVERSION version;
146 version.set_initialized();
147 version.major = 1;
148 version.minor = 0;
149 os << version;
150 }
151 else
152 os.open(m_asciiFile,ios::out|ios::app);
153
154 try {
155 os << asciiEvt;
156 }
157 catch (AsciiWrongTag& ex) {
158 std::cerr << "wrong tag, got " << ex.got()
159 << " expected: " << ex.expected()
160 << std::endl;
161 } catch (AsciiDumpException& ) {
162 std::cerr << "AsciiDumpException was caught!" << std::endl;
163 }
164 os.close();
165}
166
167void BesAsciiIO::SaveDecayMode(EVENT& asciiEvt)
168{
169 asciiEvt.decayMode.size=10;
170 //interface to event data service
171 ISvcLocator* svcLocator = Gaudi::svcLocator();
172 IDataProviderSvc* evtSvc;
173 StatusCode sc=svcLocator->service("EventDataSvc", evtSvc);
174 if (sc.isFailure())
175 G4cout<<"Could not accesss EventDataSvc!"<<G4endl;
176
177 SmartDataPtr<DecayMode> decayMode(evtSvc,"/Event/MC/DecayMode");
178 if(!decayMode)
179 {
180 for(int i=0;i<10;i++)
181 asciiEvt.decayMode.data[i]=0;
182 }
183 else
184 {
185 int dm[10]={0,0,0,0,0,0,0,0,0,0};
186 decayMode->getData(dm,10);
187 for(int i=0;i<10;i++)
188 asciiEvt.decayMode.data[i]=dm[i];
189 }
190}
191
192void BesAsciiIO::SaveTrackTruth(EVENT& asciiEvt)
193{
194
195 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
196
197 vector<BesTruthTrack*>* trackList = sensitiveManager->GetTrackList();
198
199 //arrange TruthTrack in trackList in order of trackIndex
200 G4int nTrack = trackList->size();
201 BesTruthTrack* track;
202 for(int i=0;i<nTrack-1;i++)
203 for(int j=i+1;j<nTrack;j++)
204 if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
205 {
206 track=(*trackList)[i];
207 (*trackList)[i]=(*trackList)[j];
208 (*trackList)[j]=track;
209 }
210
211 for(int i=0;i<nTrack;i++)
212 {
213 TrackTruthType trackTruth;
214 track = (*(sensitiveManager->GetTrackList())) [i];
215 trackTruth.trackIndex = track->GetIndex();
216 if(track->GetIndex()<0)
217 G4cout<<"trackIndex<0!"<<G4endl;
218 trackTruth.PDGCode = track->GetPDGCode();
219 trackTruth.PDGCharge = track->GetPDGCharge();
220 trackTruth.v0Index = track->GetVertex()->GetIndex();
221 if( track->GetTerminalVertex() )
222 trackTruth.v1Index = track->GetTerminalVertex()->GetIndex();
223 else
224 trackTruth.v1Index = -99;
225 trackTruth.px = track->GetP4().x()/1000.;
226 trackTruth.py = track->GetP4().y()/1000.;
227 trackTruth.pz = track->GetP4().z()/1000.;
228 trackTruth.E = track->GetP4().t()/1000.;
229
230 G4int size = track->GetDaughterIndexes().size();
231 if(size>0)
232 {
233 trackTruth.minDaughterIndex = (track->GetDaughterIndexes())[0];
234 trackTruth.maxDaughterIndex = (track->GetDaughterIndexes())[size-1];
235 }
236 else
237 {
238 trackTruth.minDaughterIndex = -99;
239 trackTruth.maxDaughterIndex = -99;
240 }
241
242
243 asciiEvt.trackTruth.truthCol.push_back(trackTruth);
244 }
245
246 asciiEvt.trackTruth.nTruth=asciiEvt.trackTruth.truthCol.size();
247}
248
249void BesAsciiIO::SaveVertexTruth(EVENT& asciiEvt)
250{
251 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
252
253 vector<BesTruthVertex*>* vertexList = sensitiveManager->GetVertexList();
254
255 G4int nVertex = vertexList->size();
256 BesTruthVertex* vertex;
257 for(int i=0;i<nVertex;i++)
258 {
259 VertexTruthType vertexTruth;
260 vertex = (*vertexList) [i];
261 vertexTruth.vertexIndex = vertex->GetIndex();
262 if(vertex->GetParentTrack())
263 vertexTruth.parentTrackIndex = vertex->GetParentTrack()->GetIndex();
264 else
265 vertexTruth.parentTrackIndex = -99;
266
267 vertexTruth.posX = vertex->GetPosition().x()/10.;
268 vertexTruth.posY = vertex->GetPosition().y()/10.;
269 vertexTruth.posZ = vertex->GetPosition().z()/10.;
270 vertexTruth.time = vertex->GetTime();
271 asciiEvt.vertexTruth.truthCol.push_back(vertexTruth);
272 }
273 asciiEvt.vertexTruth.nTruth = asciiEvt.vertexTruth.truthCol.size();
274}
275
276void BesAsciiIO::SaveMdcTruth(EVENT& asciiEvt)
277{
278 G4int HCID = -1;
279 HCID = m_DigiMan->GetHitsCollectionID("BesMdcTruthCollection");
280 if(HCID>0)
281 {
282 BesMdcHitsCollection* HC = 0;
283 HC = (BesMdcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
284 G4int n_hit = HC->entries();
285 if(n_hit>0)
286 {
287 //arrange hits in hits collection in order of trackIndex
288 BesMdcHit* hit;
289 vector<BesMdcHit*>* vecHC = HC->GetVector();
290 for(int i=0;i<n_hit-1;i++)
291 for(int j=i+1;j<n_hit;j++)
292 if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
293 {
294 hit = (*vecHC)[i];
295 (*vecHC)[i] = (*vecHC)[j];
296 (*vecHC)[j] = hit;
297 }
298
299 for(G4int i=0;i<n_hit;i++)
300 {
301 hit = (*HC)[i];
302 MdcTruthType mdcTruth;
303 mdcTruth.trackIndex = hit->GetTrackID();
304 mdcTruth.layerNo = hit->GetLayerNo();
305 mdcTruth.cellNo = hit->GetCellNo();
306 mdcTruth.edep = hit->GetEdep();
307 mdcTruth.driftD = hit->GetDriftD();
308 mdcTruth.posX = hit->GetPos().x();
309 mdcTruth.posY = hit->GetPos().y();
310 mdcTruth.posZ = hit->GetPos().z();
311 mdcTruth.posFlag = hit->GetPosFlag();
312
313 asciiEvt.mdcTruth.truthCol.push_back(mdcTruth);
314 }
315 }
316 }
317 asciiEvt.mdcTruth.nTruth = asciiEvt.mdcTruth.truthCol.size();
318}
319
320void BesAsciiIO::SaveCgemTruth(EVENT& asciiEvt)
321{
322 G4int HCID = -1;
323 HCID = m_DigiMan->GetHitsCollectionID("BesCgemTruthCollection");
324 if(HCID>0)
325 {
326 BesCgemHitsCollection* HC = 0;
327 HC = (BesCgemHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
328 G4int n_hit = HC->entries();
329 if(n_hit>0)
330 {
331 //arrange hits in hits collection in order of trackIndex
332 BesCgemHit* hit;
333 vector<BesCgemHit*>* vecHC = HC->GetVector();
334 for(int i=0;i<n_hit-1;i++)
335 for(int j=i+1;j<n_hit;j++)
336 if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
337 {
338 hit = (*vecHC)[i];
339 (*vecHC)[i] = (*vecHC)[j];
340 (*vecHC)[j] = hit;
341 }
342
343 for(G4int i=0;i<n_hit;i++)
344 {
345 hit = (*HC)[i];
346 CgemTruthType cgemTruth;
347
348 cgemTruth.m_ID_track = hit->GetTrackID ();
349 cgemTruth.m_ID_layer = hit->GetLayerID ();
350 cgemTruth.m_pdg_code = hit->GetPDGCode ();
351 cgemTruth.m_global_time = hit->GetGlobalTime ();
352 cgemTruth.m_E_deposit = hit->GetTotalEnergyDeposit ();
353 cgemTruth.m_L_step = hit->GetStepLength ();
354 cgemTruth.m_XYZ_pre_x = hit->GetPositionOfPrePoint ().x();
355 cgemTruth.m_XYZ_pre_y = hit->GetPositionOfPrePoint ().y();
356 cgemTruth.m_XYZ_pre_z = hit->GetPositionOfPrePoint ().z();
357 cgemTruth.m_XYZ_post_x = hit->GetPositionOfPostPoint ().x();
358 cgemTruth.m_XYZ_post_y = hit->GetPositionOfPostPoint ().y();
359 cgemTruth.m_XYZ_post_z = hit->GetPositionOfPostPoint ().z();
360 cgemTruth.m_P_pre_x = hit->GetMomentumOfPrePoint ().x();
361 cgemTruth.m_P_pre_y = hit->GetMomentumOfPrePoint ().y();
362 cgemTruth.m_P_pre_z = hit->GetMomentumOfPrePoint ().z();
363 cgemTruth.m_P_post_x = hit->GetMomentumOfPostPoint ().x();
364 cgemTruth.m_P_post_y = hit->GetMomentumOfPostPoint ().y();
365 cgemTruth.m_P_post_z = hit->GetMomentumOfPostPoint ().z();
366
367 asciiEvt.cgemTruth.truthCol.push_back(cgemTruth);
368 }
369 }
370 }
371 asciiEvt.cgemTruth.nTruth = asciiEvt.cgemTruth.truthCol.size();
372}
373
374void BesAsciiIO::SaveTofTruth(EVENT& asciiEvt)
375{
376 G4int HCID = -1;
377 HCID = m_DigiMan->GetHitsCollectionID("BesTofHitsList");
378 if(HCID>0)
379 {
380 BesTofHitsCollection* HC = 0;
381 HC = (BesTofHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
382 G4int n_hit = HC->entries();
383 if(n_hit>0)
384 {
385 //arrange hits in hits collection in order of trackIndex
386 BesTofHit* hit;
387 vector<BesTofHit*>* vecHC = HC->GetVector();
388 for(int i=0;i<n_hit-1;i++)
389 for(int j=i+1;j<n_hit;j++)
390 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
391 {
392 hit = (*vecHC)[i];
393 (*vecHC)[i] = (*vecHC)[j];
394 (*vecHC)[j] = hit;
395 }
396
397 for(G4int i=0;i<n_hit;i++)
398 {
399 hit = (*HC)[i];
400 TofTruthType tofTruth;
401 tofTruth.trackIndex = hit->GetTrackIndex();
402 tofTruth.partId = hit->GetPartId();
403 tofTruth.scinNb = hit->GetScinNb();
404 tofTruth.posX = hit->GetPos().x();
405 tofTruth.posY = hit->GetPos().y();
406 tofTruth.posZ = hit->GetPos().z();
407 tofTruth.px = hit->GetMomentum().x();
408 tofTruth.py = hit->GetMomentum().y();
409 tofTruth.pz = hit->GetMomentum().z();
410 tofTruth.trackL = hit->GetTrackL();
411 tofTruth.time = hit->GetTime();
412
413 asciiEvt.tofTruth.truthCol.push_back(tofTruth);
414 }
415 }
416 }
417 asciiEvt.tofTruth.nTruth=asciiEvt.tofTruth.truthCol.size();
418}
419
420void BesAsciiIO::SaveEmcTruth(EVENT& asciiEvt)
421{
422 G4int HCID = -1;
423 HCID = m_DigiMan->GetHitsCollectionID("BesEmcHitsList");
424 if(HCID>0)
425 {
426 BesEmcHitsCollection* HC = 0;
427 HC = (BesEmcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
428 G4int n_hit = HC->entries();
429 if(n_hit>0)
430 {
431 //arrange hits in hits collection in order of trackIndex
432 BesEmcHit* hit;
433 vector<BesEmcHit*>* vecHC = HC->GetVector();
434 for(int i=0;i<n_hit-1;i++)
435 for(int j=i+1;j<n_hit;j++)
436 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
437 {
438 hit = (*vecHC)[i];
439 (*vecHC)[i] = (*vecHC)[j];
440 (*vecHC)[j] = hit;
441 }
442
443 for(G4int i=0;i<n_hit;i++)
444 {
445 hit = (*HC)[i];
446 EmcTruthType emcTruth;
447 emcTruth.trackIndex = hit->GetTrackIndex();
448 emcTruth.partId = hit->GetPartId();
449 emcTruth.numTheta = hit->GetNumThetaCrystal();
450 emcTruth.numPhi = hit->GetNumPhiCrystal();
451 emcTruth.posX = hit->GetPosCrystal().x();
452 emcTruth.posY = hit->GetPosCrystal().y();
453 emcTruth.posZ = hit->GetPosCrystal().z();
454 emcTruth.px = hit->GetMomentum().x();
455 emcTruth.py = hit->GetMomentum().y();
456 emcTruth.pz = hit->GetMomentum().z();
457 emcTruth.totalEdep = hit->GetEdepCrystal();
458
459 asciiEvt.emcTruth.truthCol.push_back(emcTruth);
460 }
461 }
462 }
463 asciiEvt.emcTruth.nTruth=asciiEvt.emcTruth.truthCol.size();
464}
465
466void BesAsciiIO::SaveMucTruth(EVENT& asciiEvt)
467{
468 G4int HCID = -1;
469 HCID = m_DigiMan->GetHitsCollectionID("BesMucHitsList");
470 if(HCID>0)
471 {
472 BesMucHitsCollection* HC = 0;
473 HC = (BesMucHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
474 G4int n_hit = HC->entries();
475 if(n_hit>0)
476 {
477 //arrange hits in hits collection in order of trackIndex
478 BesMucHit* hit;
479 vector<BesMucHit*>* vecHC = HC->GetVector();
480 for(int i=0;i<n_hit-1;i++)
481 for(int j=i+1;j<n_hit;j++)
482 if((*vecHC)[i]->GetTrackIndex()>(*vecHC)[j]->GetTrackIndex())
483 {
484 hit = (*vecHC)[i];
485 (*vecHC)[i] = (*vecHC)[j];
486 (*vecHC)[j] = hit;
487 }
488
489 for(G4int i=0;i<n_hit;i++)
490 {
491 hit = (*HC)[i];
492 MucTruthType mucTruth;
493 mucTruth.trackIndex = hit->GetTrackIndex();
494 mucTruth.partId = hit->GetPart();
495 mucTruth.segId = hit->GetSeg();
496 mucTruth.gapId = hit->GetGap();
497 mucTruth.stripId = hit->GetStrip();
498 mucTruth.posX = hit->GetPos().x();
499 mucTruth.posY = hit->GetPos().y();
500 mucTruth.posZ = hit->GetPos().z();
501 mucTruth.px = hit->GetMomentum().x();
502 mucTruth.py = hit->GetMomentum().y();
503 mucTruth.pz = hit->GetMomentum().z();
504 asciiEvt.mucTruth.truthCol.push_back(mucTruth);
505 }
506 }
507 }
508 asciiEvt.mucTruth.nTruth=asciiEvt.mucTruth.truthCol.size();
509}
510
511void BesAsciiIO::SaveMdcDigits(EVENT& asciiEvt)
512{
513 G4int mdcDigiCollID = -1;
514 mdcDigiCollID = m_DigiMan->GetDigiCollectionID("BesMdcDigisCollection");
515 if(mdcDigiCollID>=0){
516 BesMdcDigisCollection* mdcDC = (BesMdcDigisCollection*)m_DigiMan->GetDigiCollection(mdcDigiCollID);
517 G4int nDigi = mdcDC->entries();
518 if(nDigi>0)
519 {
520 BesMdcDigi* mdcDigi;
521 for(int i=0;i<nDigi;i++)
522 {
523 MdcDigiType mdcData;
524 mdcDigi=(*mdcDC)[i];
525 mdcData.trackIndex = mdcDigi->GetTrackID();
526 mdcData.layerNo = mdcDigi->GetLayerNo();
527 mdcData.cellNo = mdcDigi->GetCellNo();
528 mdcData.energyDeposit = mdcDigi->GetEdep();
529 mdcData.driftTime = mdcDigi->GetDriftT();
530 asciiEvt.mdcDigi.digiCol.push_back(mdcData);
531 }
532 }
533 }
534 asciiEvt.mdcDigi.nDigi=asciiEvt.mdcDigi.digiCol.size();
535}
536
537void BesAsciiIO::SaveCgemDigits(EVENT& asciiEvt)
538{
539 G4int cgemDigiCollID = -1;
540 cgemDigiCollID = m_DigiMan->GetDigiCollectionID("BesCgemDigisCollection");
541 if(cgemDigiCollID>=0)
542 {
543 BesCgemDigisCollection* cgemDC = (BesCgemDigisCollection*)m_DigiMan->GetDigiCollection(cgemDigiCollID);
544 G4int nDigi = cgemDC->entries();
545 if(nDigi>0)
546 {
547 BesCgemDigi* cgemDigi;
548 for(int i=0;i<nDigi;i++)
549 {
550 CgemDigiType cgemData;
551 cgemDigi=(*cgemDC)[i];
552
553 cgemData.m_ID_track = cgemDigi->GetTrackID ();
554 cgemData.m_ID_layer = cgemDigi->GetLayerID ();
555 cgemData.m_ID_sheet = cgemDigi->GetSheetID ();
556 cgemData.m_F_XV = cgemDigi->GetStripType ();
557 cgemData.m_ID_strip = cgemDigi->GetStripID ();
558 cgemData.m_E_deposit = cgemDigi->GetEnergyDeposit ();
559 cgemData.m_global_time = cgemDigi->GetGlobalTime ();
560
561 asciiEvt.cgemDigi.digiCol.push_back(cgemData);
562 }
563 }
564 }
565 asciiEvt.cgemDigi.nDigi=asciiEvt.cgemDigi.digiCol.size();
566}
567
568
569void BesAsciiIO::SaveTofDigits(EVENT& asciiEvt)
570{
571
572 G4int tofDigiCollID = -1;
573 tofDigiCollID = m_DigiMan->GetDigiCollectionID("BesTofDigitsCollection");
574
575 if(tofDigiCollID>=0)
576 {
577 BesTofDigitsCollection* tofDC = (BesTofDigitsCollection*)m_DigiMan->GetDigiCollection(tofDigiCollID);
578 G4int nDigi = tofDC->entries();
579 if(nDigi>0)
580 {
581 //arrange digis in digitsCollection in order of trackIndex
582 BesTofDigi* digi;
583 vector<BesTofDigi*>* vecDC = tofDC->GetVector();
584 for(int i=0;i<nDigi-1;i++)
585 for(int j=i+1;j<nDigi;j++)
586 if((*vecDC)[i]->GetTrackIndex()>(*vecDC)[j]->GetTrackIndex())
587 {
588 digi = (*vecDC)[i];
589 (*vecDC)[i] = (*vecDC)[j];
590 (*vecDC)[j] = digi;
591 }
592
593 for(int i=0;i<nDigi;i++)
594 {
595 TofDigiType tofData;
596 digi = (*tofDC)[i];
597 tofData.trackIndex = digi->GetTrackIndex();
598 tofData.partId = digi->GetPartId();
599 tofData.scinNb = digi->GetScinNb();
600 tofData.forwADC = digi->GetForwADC();
601 tofData.forwTDC = digi->GetForwTDC();
602 tofData.backADC = digi->GetBackADC();
603 tofData.backTDC = digi->GetBackTDC();
604
605 asciiEvt.tofDigi.digiCol.push_back(tofData);
606 }
607 }
608 }
609 asciiEvt.tofDigi.nDigi=asciiEvt.tofDigi.digiCol.size();
610}
611
612void BesAsciiIO::SaveEmcDigits(EVENT& asciiEvt)
613{
614 G4int emcDigiCollID = -1;
615 emcDigiCollID = m_DigiMan->GetDigiCollectionID("BesEmcDigitsCollection");
616 if(emcDigiCollID>=0)
617 {
618 BesEmcDigitsCollection* emcDC = (BesEmcDigitsCollection*)m_DigiMan->GetDigiCollection(emcDigiCollID);
619 G4int nDigi = emcDC->entries();
620 if(nDigi>0)
621 {
622 //arrange digis in digitsCollection in order of trackIndex
623 BesEmcDigi* digi;
624 vector<BesEmcDigi*>* vecDC = emcDC->GetVector();
625 for(int i=0;i<nDigi-1;i++)
626 for(int j=i+1;j<nDigi;j++)
627 if((*vecDC)[i]->GetTrackIndex()>(*vecDC)[j]->GetTrackIndex())
628 {
629 digi = (*vecDC)[i];
630 (*vecDC)[i] = (*vecDC)[j];
631 (*vecDC)[j] = digi;
632 }
633 for(int i=0;i<nDigi;i++)
634 {
635 EmcDigiType emcData;
636 digi = (*emcDC)[i];
637 emcData.trackIndex = digi->GetTrackIndex();
638 emcData.partId = digi->GetPartId();
639 emcData.numTheta = digi->GetThetaNb();
640 emcData.numPhi = digi->GetPhiNb();
641 emcData.energyDeposit = digi->GetEnergy();
642 emcData.hitTime = (G4double)digi->GetTime();
643 asciiEvt.emcDigi.digiCol.push_back(emcData);
644 }
645 }
646 }
647 asciiEvt.emcDigi.nDigi=asciiEvt.emcDigi.digiCol.size();
648}
649
650void BesAsciiIO::SaveMucDigits(EVENT& asciiEvt)
651{
652 G4int mucDigiCollID =-1;
653 mucDigiCollID = m_DigiMan->GetDigiCollectionID("BesMucDigisCollection");
654 if(mucDigiCollID>=0)
655 {
656 BesMucDigisCollection* mucDC = (BesMucDigisCollection*)m_DigiMan->GetDigiCollection(mucDigiCollID);
657 G4int nDigi = mucDC->entries();
658 if(nDigi > 0) {
659 BesMucDigi* mucDigi;
660 for(int i = 0; i < nDigi; i++)
661 {
662 MucDigiType mucData;
663 mucDigi = (*mucDC)[i];
664 mucData.trackIndex = mucDigi->GetTrackIndex();
665 mucData.partNo = mucDigi->GetPartId();
666 mucData.segNo = mucDigi->GetSegId();
667 mucData.gapNo = mucDigi->GetGapId();
668 mucData.stripNo = mucDigi->GetStripId();
669
670 asciiEvt.mucDigi.digiCol.push_back(mucData);
671 }
672 }
673 }
674 asciiEvt.mucDigi.nDigi=asciiEvt.mucDigi.digiCol.size();
675}
676
677//Below used when output hits not digis
678void BesAsciiIO::SaveHitAsciiEvents(G4int runId, const G4Event* evt){
679 HitEVENT asciiEvt;
680
681 asciiEvt.set_initialized();
682 asciiEvt.header.set_initialized();
683 asciiEvt.header.eventNo=evt->GetEventID();
684 asciiEvt.header.runNo= runId;
685 asciiEvt.decayMode.set_initialized();
686 SaveDecayMode(asciiEvt);
687
688 asciiEvt.trackTruth.set_initialized();
689 asciiEvt.vertexTruth.set_initialized();
690 SaveTrackTruth(asciiEvt);
691 SaveVertexTruth(asciiEvt);
692
693 if(m_mdcTruFlag)
694 {
695 asciiEvt.mdcTruth.set_initialized();
696 SaveMdcTruth(asciiEvt);
697 }
698
699 if(m_mdcDigiFlag)
700 {
701 asciiEvt.mdcHit.set_initialized();
702 SaveMdcHits(asciiEvt);
703 }
704
705 if(m_cgemTruFlag)
706 {
707 asciiEvt.cgemTruth.set_initialized();
708 SaveCgemTruth(asciiEvt);
709 }
710
711 if(m_cgemDigiFlag)
712 {
713 asciiEvt.cgemHit.set_initialized();
714 SaveCgemHits(asciiEvt);
715 }
716 /*
717 if(m_tofTruFlag)
718 {
719 asciiEvt.tofTruth.set_initialized();
720 SaveTofTruth(asciiEvt);
721 }
722
723 if(m_tofDigiFlag)
724 {
725 asciiEvt.tofHit.set_initialized();
726 SaveTofHits(asciiEvt);
727 }
728
729 if(m_emcTruFlag)
730 {
731 asciiEvt.emcTruth.set_initialized();
732 SaveEmcTruth(asciiEvt);
733 }
734
735 if(m_emcDigiFlag)
736 {
737 asciiEvt.emcHit.set_initialized();
738 SaveEmcHits(asciiEvt);
739 }
740
741 if(m_mucTruFlag)
742 {
743 asciiEvt.mucTruth.set_initialized();
744 SaveMucTruth(asciiEvt);
745 }
746
747 if(m_mucDigiFlag)
748 {
749 asciiEvt.mucHit.set_initialized();
750 SaveMucHits(asciiEvt);
751 }
752 */
753 ofstream os;
754 if(evt->GetEventID()==0){
755 os.open(m_asciiFile);
756 FRMTVERSION version;
757 version.set_initialized();
758 version.major = 1;
759 version.minor = 0;
760 os << version;
761 }
762 else
763 os.open(m_asciiFile,ios::out|ios::app);
764
765 try {
766 os << asciiEvt;
767 }catch (AsciiWrongTag& ex) {
768 std::cerr << "wrong tag, got " << ex.got()
769 << " expected: " << ex.expected()
770 << std::endl;
771 } catch (AsciiDumpException& ) {
772 std::cerr << "AsciiDumpException was caught!" << std::endl;
773 }
774
775 os.close();
776}
777
778void BesAsciiIO::SaveDecayMode(HitEVENT& asciiEvt)
779{
780 asciiEvt.decayMode.size=10;
781 //interface to event data service
782 ISvcLocator* svcLocator = Gaudi::svcLocator();
783 IDataProviderSvc* evtSvc;
784 StatusCode sc=svcLocator->service("EventDataSvc", evtSvc);
785 if (sc.isFailure())
786 G4cout<<"Could not accesss EventDataSvc!"<<G4endl;
787
788 SmartDataPtr<DecayMode> decayMode(evtSvc,"/Event/MC/DecayMode");
789 if(!decayMode)
790 {
791 for(int i=0;i<10;i++)
792 asciiEvt.decayMode.data[i]=0;
793 }
794 else
795 {
796 int dm[10]={0,0,0,0,0,0,0,0,0,0};
797 decayMode->getData(dm,10);
798 for(int i=0;i<10;i++)
799 asciiEvt.decayMode.data[i]=dm[i];
800 }
801}
802
803void BesAsciiIO::SaveTrackTruth(HitEVENT& asciiEvt)
804{
805
806 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
807
808 vector<BesTruthTrack*>* trackList = sensitiveManager->GetTrackList();
809
810 //arrange TruthTrack in trackList in order of trackIndex
811 G4int nTrack = trackList->size();
812 BesTruthTrack* track;
813 for(int i=0;i<nTrack-1;i++)
814 for(int j=i+1;j<nTrack;j++)
815 if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
816 {
817 track=(*trackList)[i];
818 (*trackList)[i]=(*trackList)[j];
819 (*trackList)[j]=track;
820 }
821
822 for(int i=0;i<nTrack;i++)
823 {
824 TrackTruthType trackTruth;
825 track = (*(sensitiveManager->GetTrackList())) [i];
826 trackTruth.trackIndex = track->GetIndex();
827 trackTruth.PDGCode = track->GetPDGCode();
828 trackTruth.PDGCharge = track->GetPDGCharge();
829 trackTruth.v0Index = track->GetVertex()->GetIndex();
830 if( track->GetTerminalVertex() )
831 trackTruth.v1Index = track->GetTerminalVertex()->GetIndex();
832 else
833 trackTruth.v1Index = -99;
834 trackTruth.px = track->GetP4().x();
835 trackTruth.py = track->GetP4().y();
836 trackTruth.pz = track->GetP4().z();
837 trackTruth.E = track->GetP4().t();
838
839 G4int size = track->GetDaughterIndexes().size();
840 if(size>0)
841 {
842 trackTruth.minDaughterIndex = (track->GetDaughterIndexes())[0];
843 trackTruth.maxDaughterIndex = (track->GetDaughterIndexes())[size-1];
844 }
845 else
846 {
847 trackTruth.minDaughterIndex = -99;
848 trackTruth.maxDaughterIndex = -99;
849 }
850
851
852 asciiEvt.trackTruth.truthCol.push_back(trackTruth);
853 }
854
855 asciiEvt.trackTruth.nTruth=asciiEvt.trackTruth.truthCol.size();
856}
857
858void BesAsciiIO::SaveVertexTruth(HitEVENT& asciiEvt)
859{
860 BesSensitiveManager* sensitiveManager = BesSensitiveManager::GetSensitiveManager();
861
862 vector<BesTruthVertex*>* vertexList = sensitiveManager->GetVertexList();
863
864 G4int nVertex = vertexList->size();
865 BesTruthVertex* vertex;
866 for(int i=0;i<nVertex;i++)
867 {
868 VertexTruthType vertexTruth;
869 vertex = (*vertexList) [i];
870 vertexTruth.vertexIndex = vertex->GetIndex();
871 if(vertex->GetParentTrack())
872 vertexTruth.parentTrackIndex = vertex->GetParentTrack()->GetIndex();
873 else
874 vertexTruth.parentTrackIndex = -99;
875
876 vertexTruth.posX = vertex->GetPosition().x();
877 vertexTruth.posY = vertex->GetPosition().y();
878 vertexTruth.posZ = vertex->GetPosition().z();
879 vertexTruth.time = vertex->GetTime();
880 asciiEvt.vertexTruth.truthCol.push_back(vertexTruth);
881 }
882 asciiEvt.vertexTruth.nTruth = asciiEvt.vertexTruth.truthCol.size();
883}
884
885void BesAsciiIO::SaveMdcTruth(HitEVENT& asciiEvt)
886{
887 G4int HCID = -1;
888 HCID = m_DigiMan->GetHitsCollectionID("BesMdcTruthCollection");
889 if(HCID>0)
890 {
891 BesMdcHitsCollection* HC = 0;
892 HC = (BesMdcHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
893 G4int n_hit = HC->entries();
894 if(n_hit>0)
895 {
896 //arrange hits in hits collection in order of trackIndex
897 BesMdcHit* hit;
898 vector<BesMdcHit*>* vecHC = HC->GetVector();
899 for(int i=0;i<n_hit-1;i++)
900 for(int j=i+1;j<n_hit;j++)
901 if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
902 {
903 hit = (*vecHC)[i];
904 (*vecHC)[i] = (*vecHC)[j];
905 (*vecHC)[j] = hit;
906 }
907
908 for(G4int i=0;i<n_hit;i++)
909 {
910 hit = (*HC)[i];
911 MdcTruthType mdcTruth;
912 mdcTruth.trackIndex = hit->GetTrackID();
913 mdcTruth.layerNo = hit->GetLayerNo();
914 mdcTruth.cellNo = hit->GetCellNo();
915 mdcTruth.edep = hit->GetEdep();
916 mdcTruth.driftD = hit->GetDriftD();
917 mdcTruth.posX = hit->GetPos().x();
918 mdcTruth.posY = hit->GetPos().y();
919 mdcTruth.posZ = hit->GetPos().z();
920 mdcTruth.posFlag = hit->GetPosFlag();
921
922 asciiEvt.mdcTruth.truthCol.push_back(mdcTruth);
923 }
924 }
925 }
926 asciiEvt.mdcTruth.nTruth = asciiEvt.mdcTruth.truthCol.size();
927}
928
929void BesAsciiIO::SaveCgemTruth(HitEVENT& asciiEvt)
930{
931 G4int HCID = -1;
932 HCID = m_DigiMan->GetHitsCollectionID("BesCgemTruthCollection");
933 if(HCID>0)
934 {
935 BesCgemHitsCollection* HC = 0;
936 HC = (BesCgemHitsCollection*) (m_DigiMan->GetHitsCollection(HCID));
937 G4int n_hit = HC->entries();
938 if(n_hit>0)
939 {
940 //arrange hits in hits collection in order of trackIndex
941 BesCgemHit* hit;
942 vector<BesCgemHit*>* vecHC = HC->GetVector();
943 for(int i=0;i<n_hit-1;i++)
944 for(int j=i+1;j<n_hit;j++)
945 if((*vecHC)[i]->GetTrackID()>(*vecHC)[j]->GetTrackID())
946 {
947 hit = (*vecHC)[i];
948 (*vecHC)[i] = (*vecHC)[j];
949 (*vecHC)[j] = hit;
950 }
951
952 for(G4int i=0;i<n_hit;i++)
953 {
954 hit = (*HC)[i];
955 CgemTruthType cgemTruth;
956
957 cgemTruth.m_ID_track = hit->GetTrackID ();
958 cgemTruth.m_ID_layer = hit->GetLayerID ();
959 cgemTruth.m_pdg_code = hit->GetPDGCode ();
960 cgemTruth.m_global_time = hit->GetGlobalTime ();
961 cgemTruth.m_E_deposit = hit->GetTotalEnergyDeposit ();
962 cgemTruth.m_L_step = hit->GetStepLength ();
963 cgemTruth.m_XYZ_pre_x = hit->GetPositionOfPrePoint ().x();
964 cgemTruth.m_XYZ_pre_y = hit->GetPositionOfPrePoint ().y();
965 cgemTruth.m_XYZ_pre_z = hit->GetPositionOfPrePoint ().z();
966 cgemTruth.m_XYZ_post_x = hit->GetPositionOfPostPoint ().x();
967 cgemTruth.m_XYZ_post_y = hit->GetPositionOfPostPoint ().y();
968 cgemTruth.m_XYZ_post_z = hit->GetPositionOfPostPoint ().z();
969 cgemTruth.m_P_pre_x = hit->GetMomentumOfPrePoint ().x();
970 cgemTruth.m_P_pre_y = hit->GetMomentumOfPrePoint ().y();
971 cgemTruth.m_P_pre_z = hit->GetMomentumOfPrePoint ().z();
972 cgemTruth.m_P_post_x = hit->GetMomentumOfPostPoint ().x();
973 cgemTruth.m_P_post_y = hit->GetMomentumOfPostPoint ().y();
974 cgemTruth.m_P_post_z = hit->GetMomentumOfPostPoint ().z();
975
976 asciiEvt.cgemTruth.truthCol.push_back(cgemTruth);
977 }
978 }
979 }
980 asciiEvt.cgemTruth.nTruth = asciiEvt.cgemTruth.truthCol.size();
981}
982
983void BesAsciiIO::SaveMdcHits(HitEVENT& asciiEvt)
984{
985 G4int mdcHitCollID = -1;
986 mdcHitCollID = m_DigiMan->GetHitsCollectionID("BesMdcHitsCollection");
987 if(mdcHitCollID>=0)
988 {
989 BesMdcHitsCollection* mdcDC = (BesMdcHitsCollection*)m_DigiMan->GetHitsCollection(mdcHitCollID);
990 G4int nHit = mdcDC->entries();
991 if(nHit>0)
992 {
993 BesMdcHit* mdcHit;
994 for(int i=0;i<nHit;i++)
995 {
996 MdcHitType mdcData;
997 mdcHit=(*mdcDC)[i];
998 mdcData.trackIndex = mdcHit->GetTrackID();
999 mdcData.layerNo = mdcHit->GetLayerNo();
1000 mdcData.cellNo = mdcHit->GetCellNo();
1001 mdcData.posX = mdcHit->GetPos().x();
1002 mdcData.posY = mdcHit->GetPos().y();
1003 mdcData.posZ = mdcHit->GetPos().z();
1004 mdcData.energyDeposit = mdcHit->GetEdep();
1005 mdcData.driftDistance = mdcHit->GetDriftD();
1006 mdcData.globalT = mdcHit->GetGlobalT();
1007 mdcData.theta = mdcHit->GetTheta();
1008 mdcData.enterAngle = mdcHit->GetEnterAngle();
1009 mdcData.posFlag = mdcHit->GetPosFlag();
1010
1011 asciiEvt.mdcHit.hitCol.push_back(mdcData);
1012 }
1013 }
1014 }
1015 asciiEvt.mdcHit.nHit=asciiEvt.mdcHit.hitCol.size();
1016}
1017
1018void BesAsciiIO::SaveCgemHits(HitEVENT& asciiEvt)
1019{
1020 G4int cgemHitCollID = -1;
1021 cgemHitCollID = m_DigiMan->GetHitsCollectionID("BesCgemHitsCollection");
1022 if(cgemHitCollID>=0)
1023 {
1024 BesCgemHitsCollection* cgemDC = (BesCgemHitsCollection*)m_DigiMan->GetHitsCollection(cgemHitCollID);
1025 G4int nHit = cgemDC->entries();
1026 if(nHit>0)
1027 {
1028 BesCgemHit* cgemHit;
1029 for(int i=0;i<nHit;i++)
1030 {
1031 CgemHitType cgemData;
1032 cgemHit=(*cgemDC)[i];
1033
1034 cgemData.m_ID_track = cgemHit->GetTrackID ();
1035 cgemData.m_ID_layer = cgemHit->GetLayerID ();
1036 cgemData.m_pdg_code = cgemHit->GetPDGCode ();
1037 cgemData.m_global_time = cgemHit->GetGlobalTime ();
1038 cgemData.m_E_deposit = cgemHit->GetTotalEnergyDeposit ();
1039 cgemData.m_L_step = cgemHit->GetStepLength ();
1040 cgemData.m_XYZ_pre_x = cgemHit->GetPositionOfPrePoint ().x();
1041 cgemData.m_XYZ_pre_y = cgemHit->GetPositionOfPrePoint ().y();
1042 cgemData.m_XYZ_pre_z = cgemHit->GetPositionOfPrePoint ().z();
1043 cgemData.m_XYZ_post_x = cgemHit->GetPositionOfPostPoint ().x();
1044 cgemData.m_XYZ_post_y = cgemHit->GetPositionOfPostPoint ().y();
1045 cgemData.m_XYZ_post_z = cgemHit->GetPositionOfPostPoint ().z();
1046 cgemData.m_P_pre_x = cgemHit->GetMomentumOfPrePoint ().x();
1047 cgemData.m_P_pre_y = cgemHit->GetMomentumOfPrePoint ().y();
1048 cgemData.m_P_pre_z = cgemHit->GetMomentumOfPrePoint ().z();
1049 cgemData.m_P_post_x = cgemHit->GetMomentumOfPostPoint ().x();
1050 cgemData.m_P_post_y = cgemHit->GetMomentumOfPostPoint ().y();
1051 cgemData.m_P_post_z = cgemHit->GetMomentumOfPostPoint ().z();
1052
1053 asciiEvt.cgemHit.hitCol.push_back(cgemData);
1054 }
1055 }
1056 }
1057 asciiEvt.cgemHit.nHit=asciiEvt.cgemHit.hitCol.size();
1058}
BesCgemDigisCollection
G4TDigiCollection< BesCgemDigi > BesCgemDigisCollection
Definition: InstallArea/include/CgemSim/CgemSim/BesCgemDigi.hh:91
BesCgemHitsCollection
G4THitsCollection< BesCgemHit > BesCgemHitsCollection
Definition: InstallArea/include/CgemSim/CgemSim/BesCgemHit.hh:149
BesEmcDigitsCollection
G4TDigiCollection< BesEmcDigi > BesEmcDigitsCollection
Definition: InstallArea/include/EmcSim/EmcSim/BesEmcDigi.hh:69
BesEmcHitsCollection
G4THitsCollection< BesEmcHit > BesEmcHitsCollection
Definition: InstallArea/include/EmcSim/EmcSim/BesEmcHit.hh:83
BesMdcDigisCollection
G4TDigiCollection< BesMdcDigi > BesMdcDigisCollection
Definition: InstallArea/include/MdcSim/MdcSim/BesMdcDigi.hh:59
BesMdcHitsCollection
G4THitsCollection< BesMdcHit > BesMdcHitsCollection
Definition: InstallArea/include/MdcSim/MdcSim/BesMdcHit.hh:78
BesMucDigisCollection
G4TDigiCollection< BesMucDigi > BesMucDigisCollection
Definition: InstallArea/include/MucSim/MucSim/BesMucDigi.hh:56
BesMucHitsCollection
G4THitsCollection< BesMucHit > BesMucHitsCollection
Definition: InstallArea/include/MucSim/MucSim/BesMucHit.hh:104
BesTofDigitsCollection
G4TDigiCollection< BesTofDigi > BesTofDigitsCollection
Definition: InstallArea/include/TofSim/TofSim/BesTofDigi.hh:83
BesTofHitsCollection
G4THitsCollection< BesTofHit > BesTofHitsCollection
Definition: InstallArea/include/TofSim/TofSim/BesTofHit.hh:108
BesAsciiIO::BesAsciiIO
BesAsciiIO(G4int, G4int, G4int, G4int, G4int, G4int, G4int, G4int, G4int, G4int, G4String)
Definition: bak_BesSim-00-04-14/src/BesAsciiIO.cc:44
BesAsciiIO::SaveHitAsciiEvents
void SaveHitAsciiEvents(G4int, const G4Event *)
Definition: bak_BesSim-00-04-14/src/BesAsciiIO.cc:678
BesAsciiIO::SaveAsciiEvents
void SaveAsciiEvents(G4int, const G4Event *)
Definition: bak_BesSim-00-04-14/src/BesAsciiIO.cc:65