#include <LCDG4EventAction.hh>
Inheritance diagram for LCDG4EventAction:
Public Member Functions | |
| LCDG4EventAction (SIO_stream *, LCDG4McPartManager *) | |
| ~LCDG4EventAction () | |
| void | SetRandSeed (G4int r) |
| G4int | GetRandSeed () |
| void | BeginOfEventAction (const G4Event *) |
| void | EndOfEventAction (const G4Event *) |
| void | SetupOutputLcio (const char *lcio_file) |
Public Attributes | |
| SIO_stream * | stream |
Private Member Functions | |
| void | SaveCalorimeterHits (const LCDG4CalHitsCollection *pHitColl, std::vector< CalorimeterHit > &calhits, IMPL::LCCollectionVec *lcioColl) const |
Private Attributes | |
| G4int | m_vertexCollID |
| G4int | m_trackerCollID |
| G4int | m_lumCollID |
| G4int | m_muonCollID |
| G4int | m_ecalBCollID |
| G4int | m_ecalECCollID |
| G4int | m_hcalBCollID |
| G4int | m_hcalECCollID |
| G4bool | m_event |
| G4int | rndSeed |
| G4int | fixedSeed |
| LCDG4McPartManager * | m_mcmgr |
| std::map< int, MCParticleSIO * > | mcpartmap |
| LCDG4LcioHelper * | _lcioHelper |
|
||||||||||||
|
Definition at line 58 of file LCDG4EventAction.cc. References m_ecalBCollID, m_ecalECCollID, m_event, m_hcalBCollID, m_hcalECCollID, m_lumCollID, m_muonCollID, m_trackerCollID, m_vertexCollID, and stream.
00059 : fixedSeed(0) 00060 , m_mcmgr(mcmgr) 00061 , _lcioHelper(NULL) 00062 #ifdef TXT_MODE 00063 , outf_txt(NULL) 00064 #endif 00065 { 00066 m_vertexCollID = -9; 00067 m_trackerCollID = -9; 00068 m_ecalBCollID = -9; 00069 m_ecalECCollID = -9; 00070 m_hcalBCollID = -9; 00071 m_hcalECCollID = -9; 00072 m_muonCollID = -9; 00073 m_lumCollID = -9; 00074 m_event = 1; 00075 stream = stream_out; 00076 // create a new dump file 00077 ofstream dumpfile("debug-helper"); 00078 dumpfile.close(); 00079 } |
|
|
Definition at line 117 of file LCDG4EventAction.cc. References _lcioHelper.
00117 {
00118 if(_lcioHelper) delete _lcioHelper;
00119 #ifdef TXT_MODE
00120 if(outf_txt) {
00121 outf_txt->close();
00122 delete outf_txt;
00123 }
00124 #endif
00125 }
|
|
|
Definition at line 127 of file LCDG4EventAction.cc. References fixedSeed, m_ecalBCollID, m_ecalECCollID, m_hcalBCollID, m_hcalECCollID, m_lumCollID, m_mcmgr, m_muonCollID, m_trackerCollID, m_vertexCollID, LCDG4McPartManager::SetDebugFlag(), and SetRandSeed().
00127 {
00128
00129 if(evt==0) return;
00130
00131 // For debugging a specific event ==================
00132 // create file debevt.dat with just an event number inside it (start from 1).
00133 // This number will be read into ievtdgb.
00134 // debugging flag in McPartMgr will be set just for event chosen
00135 static int ievt=0;
00136 static int ievtdbg;
00137 if(ievt==0) {
00138 ifstream debevtfile("debevt.dat");
00139 debevtfile>> ievtdbg >> fixedSeed;
00140 debevtfile.close();
00141 }
00142
00143 // increment event counter: 1st event is 1, not zero
00144 ievt++;
00145 if(ievt==ievtdbg) {
00146 m_mcmgr->SetDebugFlag(0x1f0);
00147 G4UImanager* UI = G4UImanager::GetUIpointer();
00148 UI->ApplyCommand("/tracking/verbose 1");
00149 }
00150
00151 // At first event only
00152 if(ievt==1) {
00153 // Initialize the random number generator
00154 // HepRandom::createInstance();
00155 // HepRandomEngine* anEngine = new RanluxEngine((unsigned)time(NULL),3);
00156 // HepRandom::setTheEngine(anEngine);
00157
00158 // Set the random seed for best randomness
00159 G4int seed = (unsigned)time(NULL);
00160 if(fixedSeed<0) {
00161 seed=-fixedSeed;
00162 G4cout<<"***** Using user provided seed at init only: seed="
00163 << seed << G4endl;
00164 }
00165 else {
00166 G4cout<<"***** Using time for best randomness: seed=" << seed << G4endl;
00167 }
00168 SetRandSeed((G4int)seed);
00169 }
00170
00171 if(fixedSeed>0) {
00172 // Use a fixed random seed for all events: for easy reproducibility
00173 G4cout<< "***** Using fixed random seed = " << fixedSeed << G4endl;
00174 SetRandSeed(fixedSeed);
00175 }
00176
00177 G4SDManager * SDman = G4SDManager::GetSDMpointer();
00178 G4String colNam;
00179 if(m_vertexCollID < -1) {
00180 m_vertexCollID = SDman->GetCollectionID(colNam="vxdCollection");
00181 }
00182 if(m_trackerCollID < -1) {
00183 m_trackerCollID = SDman->GetCollectionID(colNam="trackerCollection");
00184 }
00185 if(m_ecalBCollID < -1) {
00186 m_ecalBCollID = SDman->GetCollectionID(colNam="ecalBCollection");
00187 }
00188 if(m_ecalECCollID < -1) {
00189 m_ecalECCollID = SDman->GetCollectionID(colNam="ecalECCollection");
00190 }
00191 if(m_hcalBCollID < -1) {
00192 m_hcalBCollID = SDman->GetCollectionID(colNam="hcalBCollection");
00193 }
00194 if(m_hcalECCollID < -1) {
00195 m_hcalECCollID = SDman->GetCollectionID(colNam="hcalECCollection");
00196 }
00197 if(m_muonCollID < -1) {
00198 m_muonCollID = SDman->GetCollectionID(colNam="muonCollection");
00199 }
00200 if(m_lumCollID < -1) {
00201 m_lumCollID = SDman->GetCollectionID(colNam="lumCollection");
00202 }
00203
00204 G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
00205 if(pVVisManager) {
00206 G4UImanager::GetUIpointer()->ApplyCommand("/vis/drawTree/current");
00207 }
00208 }
|
|
|
Definition at line 210 of file LCDG4EventAction.cc. References _lcioHelper, SIO_blockManager::add(), SIO_recordManager::add(), MCParticleSIO::addChild(), MuonHit::addMcHit(), CalorimeterHit::addMcHit(), MCParticleSIO::ALIVE, LCDG4LcioHelper::buildMCParticleColl(), LCDG4McPartManager::Clear(), SIO_record::connect(), MCParticleSIO::DECAYED, LCDMcPart::Get4Momentum(), LCDMcPart::Get4MomentumPtr(), LCDG4TrackerHit::GetBarEnd(), LCDG4LcioHelper::getCalHitLcio(), LCDMcPart::GetCalorimeterEntranceMomentum(), LCDMcPart::GetCalorimeterEntrancePoint(), LCDMcPart::GetCharge(), LCDG4CalHit::GetEdep(), LCDG4TrackerHit::GetEdep(), LCDG4CalHit::GetEdepAbs(), LCDG4CalHit::GetEMC(), MCParticleSIO::getfsflag(), MCParticleSIO::GetG4Status(), LCDMcPart::GetG4Status(), LCDMcPart::GetGenStatus(), LCDG4CalHit::GetLayer(), LCDG4TrackerHit::GetLayer(), MCParticleSIO::getMcPart(), LCDG4CalHit::GetNMC(), LCDG4TrackerHit::GetNorthSouth(), LCDMcPart::GetParentTag(), LCDMcPart::GetParticleID(), LCDG4CalHit::GetPhi(), LCDG4TrackerHit::GetPos(), LCDG4McPartManager::GetShowerFlag(), LCDG4CalHit::GetSystem(), LCDG4TrackerHit::GetSystem(), LCDtowerID::GetTag(), LCDG4TrackerHit::GetTdep(), LCDMcPart::GetTermPosition(), LCDG4CalHit::GetTheta(), LCDG4TrackerHit::GetTrackID(), LCDG4LcioHelper::getTrkHitLcio(), LCDG4CalHit::GetTrkID(), MCParticleSIO::INTERACTED, LCDG4CalHitsCollection, LCDG4TrackerHitsCollection, MCParticleSIO::LEFT, MCParticleSIO::LOOPING, MCParticleSIO::LOST, m_ecalBCollID, m_ecalECCollID, m_hcalBCollID, m_hcalECCollID, m_lumCollID, m_mcmgr, m_muonCollID, m_trackerCollID, m_vertexCollID, MCParticleSIO::MAXSTEPS, LCDG4McPartManager::MCparticles, mcpartmap, myFormat(), MCParticleSIO::PRIMARY, SIO_recordManager::remove(), SaveCalorimeterHits(), LCDG4LcioHelper::saveCollection(), MCPrintSIO::set(), MCParticleSIO::setCalMomentum(), MCParticleSIO::setCalPosition(), MCParticleSIO::setfsflag(), LCDMcPart::SetG4Status(), MCParticleSIO::SetG4Status(), LCDtowerID::SetLayer(), MCParticleSIO::setMcPart(), MCParticleSIO::setMomentum(), MCParticleSIO::setParent(), LCDtowerID::SetPhi(), MCParticleSIO::setPosition(), MCParticleSIO::setStatus(), LCDtowerID::SetSystem(), LCDtowerID::SetTheta(), MCParticleSIO::SHOWER, MCParticleSIO::status(), MCParticleSIO::Status, MCParticleSIO::STOPPED, stream, MCParticleSIO::STUCK, LCDG4McPartManager::trkidToMctag(), SIO_stream::write(), and LCDG4LcioHelper::writeEvent().
00210 {
00211
00212 ofstream dumpfile("debug-helper",ios::app);
00213
00214 G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
00215 LCDG4TrackerHitsCollection* VHC = NULL;
00216 LCDG4TrackerHitsCollection* THC = NULL;
00217 LCDG4CalHitsCollection* EBHC = NULL;
00218 LCDG4CalHitsCollection* EECHC = NULL;
00219 LCDG4CalHitsCollection* HBHC = NULL;
00220 LCDG4CalHitsCollection* HECHC = NULL;
00221 LCDG4CalHitsCollection* MHC = NULL;
00222 LCDG4CalHitsCollection* LHC = NULL;
00223
00224 if(HCE) {
00225 if (m_vertexCollID >= 0) {
00226 VHC = (LCDG4TrackerHitsCollection*)(HCE->GetHC(m_vertexCollID));
00227 }
00228 if (m_trackerCollID >= 0) {
00229 THC = (LCDG4TrackerHitsCollection*)(HCE->GetHC(m_trackerCollID));
00230 }
00231 if(m_ecalBCollID >= 0) {
00232 EBHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_ecalBCollID));
00233 }
00234 if(m_ecalECCollID >= 0) {
00235 EECHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_ecalECCollID));
00236 }
00237 if(m_hcalBCollID >= 0) {
00238 HBHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_hcalBCollID));
00239 }
00240 if(m_hcalECCollID >= 0) {
00241 HECHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_hcalECCollID));
00242 }
00243 if (m_muonCollID >= 0) {
00244 MHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_muonCollID));
00245 }
00246 if (m_lumCollID >= 0) {
00247 LHC = (LCDG4CalHitsCollection*)(HCE->GetHC(m_lumCollID));
00248 }
00249 }
00250
00251 #ifdef TXT_MODE
00252 if(outf_txt) {
00253 *outf_txt << "EVENT " << evt->GetEventID() << G4endl;
00254 }
00255 #endif
00256
00257 // SIO
00258 SIO_record *record = SIO_recordManager::add("LCD_LCDG4_Event");
00259 SIO_iodev *event = new SIO_iodev(evt->GetEventID(),"EventMarker");
00260 SIO_blockManager::add(event);
00261 record->connect(event);
00262
00263 VxdSIO *vxd = NULL;
00264 TrackerSIO *tracker = NULL;
00265 CalorimeterSIO *ecalsio = NULL;
00266 CalorimeterSIO *hcalsio = NULL;
00267 CalorimeterSIO *lum = NULL;
00268 MuonSIO *muon = NULL;
00269 MCPrintSIO* mcprint = new MCPrintSIO("MCPrint");
00270
00271 //************** MC Particle Information *************
00272 vector<LCDMcPart*>& mclist = m_mcmgr->MCparticles;
00273 vector<LCDMcPart*>::const_iterator mciter = mclist.begin();
00274
00275 int mapindex = 0; // tags start at one, not zero
00276
00277 // sio format requires that first element is parent of all other elements.
00278 // Its ParticleID is set to 99999999 (HePEVT id).
00279 // Sometimes the input stdhep file already has this HepEVT element, then
00280 // keep it. Otherwise, create a new one for the .sio output
00281 // if( (*mciter)->GetParticleID()!=99999999 ) { // no, it did not work...
00282 MCParticleSIO* mcparthead = new MCParticleSIO(99999999,0,0);
00283 mcparthead->setfsflag(300);
00284 mcpartmap[mapindex] = mcparthead;
00285 mapindex++;
00286 // }
00287
00288 dumpfile <<" "<< evt->GetEventID();
00289 dumpfile <<" "<< mclist.size();
00290 for(mciter = mclist.begin(); mciter != mclist.end(); mciter++, mapindex++) {
00291 LCDMcPart* pMcPart = *mciter;
00292 if(!pMcPart) continue; // some McParts have been erased!
00293
00294 MCParticleSIO* mcp;
00295 mcp = new MCParticleSIO(pMcPart->GetParticleID(),(int)pMcPart->GetCharge(),pMcPart->Get4MomentumPtr()->e());
00296 mcp->setMcPart(pMcPart);
00297
00298 MCParticleSIO::Status Status = MCParticleSIO::ALIVE;
00299
00300 switch(pMcPart->GetG4Status()) {
00301 case 0 :
00302 Status = MCParticleSIO::ALIVE; break;
00303 case 1 :
00304 Status = MCParticleSIO::DECAYED; break;
00305 case 2 :
00306 Status = MCParticleSIO::INTERACTED; break;
00307 case 3 :
00308 Status = MCParticleSIO::LEFT; break;
00309 case 4 :
00310 Status = MCParticleSIO::STOPPED; break;
00311 case 5 :
00312 Status = MCParticleSIO::LOOPING; break;
00313 case 6 :
00314 Status = MCParticleSIO::LOST; break;
00315 case 7 :
00316 Status = MCParticleSIO::STUCK; break;
00317 case 8 :
00318 Status = MCParticleSIO::PRIMARY; break;
00319 case 9 :
00320 Status = MCParticleSIO::SHOWER; break;
00321 case 10 :
00322 Status = MCParticleSIO::MAXSTEPS; break;
00323 }
00324 mcp->setfsflag( 100*pMcPart->GetGenStatus() );
00325
00326
00327 mcp->setStatus(Status);
00328 // McPart position is set to where it dies
00329 // starting position is taken from where its parent dies
00330 // secondaries from middle of track lose info on starting point
00331 mcp->setPosition(pMcPart->GetTermPosition().x(),
00332 pMcPart->GetTermPosition().y(),
00333 pMcPart->GetTermPosition().z());
00334 mcp->setMomentum(pMcPart->Get4Momentum().px(),
00335 pMcPart->Get4Momentum().py(),
00336 pMcPart->Get4Momentum().pz());
00337 mcp->setCalPosition(pMcPart->GetCalorimeterEntrancePoint().x(),
00338 pMcPart->GetCalorimeterEntrancePoint().y(),
00339 pMcPart->GetCalorimeterEntrancePoint().z() );
00340 mcp->setCalMomentum(pMcPart->GetCalorimeterEntranceMomentum().x(),
00341 pMcPart->GetCalorimeterEntranceMomentum().y(),
00342 pMcPart->GetCalorimeterEntranceMomentum().z() );
00343
00344 // if(mapindex==1) G4cout<< "\n **** Array of LCDMcPart's ****\n";
00345 // G4cout << "*** tag="<< mapindex
00346 // <<", id="<< pMcPart->GetParticleID()
00347 // <<", status="<< pMcPart->GetGenStatus()
00348 // <<' '<< pMcPart->GetG4Status()
00349 // <<", q="<< (int)pMcPart->GetCharge()
00350 // <<", E="<< pMcPart->Get4Momentum().e()
00351 // <<", pvec=("<< pMcPart->Get4Momentum().px()
00352 // <<"; "<< pMcPart->Get4Momentum().py()
00353 // <<"; "<< pMcPart->Get4Momentum().pz()
00354 // <<") par="<< pMcPart->GetParentTag()
00355 // <<", startPos:x,y,z=("<< pMcPart->GetPositionPtr()->x()/cm
00356 // <<' '<< pMcPart->GetPositionPtr()->y()/cm
00357 // <<' '<< pMcPart->GetPositionPtr()->z()/cm <<')'
00358 // <<", endPos:x,y,z=("<< pMcPart->GetTermPositionPtr()->x()/cm
00359 // <<' '<< pMcPart->GetTermPositionPtr()->y()/cm
00360 // <<' '<< pMcPart->GetTermPositionPtr()->z()/cm <<')'
00361 // <<", x,y,z=("<< pMcPart->GetCalorimeterEntrancePoint().x()
00362 // <<' '<< pMcPart->GetCalorimeterEntrancePoint().y()
00363 // <<' '<< pMcPart->GetCalorimeterEntrancePoint().z() <<')'
00364 // <<", mom=("<< pMcPart->GetCalorimeterEntranceMomentum().x()
00365 // <<' '<< pMcPart->GetCalorimeterEntranceMomentum().y()
00366 // <<' '<< pMcPart->GetCalorimeterEntranceMomentum().z() <<')'
00367 // << G4endl;
00368
00369 mcpartmap[mapindex] = mcp;
00370 }
00371
00372 // Setup hierarchy links (parents <--> children) based on McParts
00373 for(unsigned int k = 1; k < mcpartmap.size(); k++) {
00374 MCParticleSIO* pmcp = mcpartmap[k];
00375 if(!pmcp) continue;
00376
00377 const LCDMcPart* pMcPart = pmcp->getMcPart();
00378
00379 if( pMcPart->GetParentTag()==-1) {
00380 mcpartmap[0]->addChild(pmcp);
00381 pmcp->setParent(mcpartmap[0]);
00382 }
00383 else {
00384 int parentTag = pMcPart->GetParentTag();
00385 MCParticleSIO* parent = mcpartmap[parentTag];
00386 parent->addChild(pmcp);
00387 pmcp->setParent(parent);
00388
00389 // set parent status to DECAYED
00390 if(parent->GetG4Status()==MCParticleSIO::ALIVE) {
00391 // for LCDMcPart
00392 if(parentTag>0) mclist[parentTag-1]->SetG4Status(1);
00393 // and for MCParticleSIO
00394 parent->SetG4Status(MCParticleSIO::DECAYED);
00395 }
00396 }
00397 }
00398
00399 // set other status codes
00400 mcpartmap[0]->setStatus( MCParticleSIO::PRIMARY );
00401 for(unsigned int k = 1; k < mcpartmap.size(); k++) {
00402 MCParticleSIO* pmcp = mcpartmap[k];
00403 if(!pmcp) continue;
00404
00405 LCDMcPart* pMcPart = mclist[k-1];
00406
00407 // Documentation particles
00408 if( pmcp->getfsflag() == 300 ) {
00409 // HEPEVT status DOCUMENTATION
00410 pmcp->setStatus( MCParticleSIO::PRIMARY );
00411 pMcPart->SetG4Status(8);
00412 // continue;
00413 }
00414 else if( pmcp->getfsflag() == 200 ) {
00415 // HEPEVT status INTERMEDIATE
00416 if( pmcp->status() == MCParticleSIO::ALIVE ) {
00417 // several gluons and quarks show up here
00418 pmcp->setStatus( MCParticleSIO::PRIMARY );
00419 pMcPart->SetG4Status(8);
00420 }
00421 }
00422
00423 // LCDMcPart* lcdmcpart = const_cast<LCDMcPart*>(pmcp->getMcPart());
00424 // Hep3Vector& hep3v = lcdmcpart->GetPosition();
00425 // if(k==1) G4cout<< "\n **** Array of MCParticleSIO's ****\n";
00426 // G4cout << pmcp
00427 // <<", tag="<< k <<", id="<< pmcp->id()
00428 // <<", stat(Gen,Sim)="<< pmcp->GetGenStatus() <<','<< pmcp->GetG4Status()
00429 // <<", q="<< pmcp->charge()
00430 // <<", E="<< pmcp->energy()
00431 // <<", xyzProd=("<< hep3v.x()/mm
00432 // <<' '<< hep3v.y()/mm
00433 // <<' '<< hep3v.z()/mm
00434 // <<", xyzEnd=("<< pmcp->x()/mm
00435 // <<' '<< pmcp->y()/mm
00436 // <<' '<< pmcp->z()/mm <<')'
00437 // <<", pvec=("<< pmcp->px()/GeV
00438 // <<' '<< pmcp->py()/GeV
00439 // <<' '<< pmcp->pz()/GeV <<')'
00440 // <<", Calxyz=("<< pmcp->calX()/mm
00441 // <<' '<< pmcp->calY()/mm
00442 // <<' '<< pmcp->calZ()/mm <<')'
00443 // <<", Calmom=("<< pmcp->calPX()/GeV
00444 // <<' '<< pmcp->calPY()/GeV
00445 // <<' '<< pmcp->calPZ()/GeV <<')'
00446 // <<", parID="<< pmcp->parent()->id()
00447 // <<", parE="<< pmcp->parent()->energy()/GeV
00448 // <<' '<< pmcp->parent()
00449 // << G4endl << G4endl;
00450
00451 // G4cout << "daulist ("<< pmcp->numChildren() <<" kids)";
00452 // for(int ii=0; ii<pmcp->numChildren(); ii++) {
00453 // const MCParticleSIO* kid = pmcp->child(ii);
00454 // G4cout <<" : "<< kid->id() <<' '<< kid;
00455 // }
00456 // G4cout << G4endl;
00457
00458 }
00459
00460 // Setup MCParticle table for output to SIO file
00461 mcprint->set(mcpartmap[0]);
00462
00463 // Build the particle table for LCIO
00464 if(_lcioHelper) _lcioHelper->buildMCParticleColl(mclist);
00465
00466
00467 //************** VXD hits *************
00468 if(VHC) {
00469 int n_hit = VHC->entries();
00470 dumpfile << " "<< n_hit;
00471 vector<TkHit> vxdhits;
00472
00473 LCCollectionVec* lcioVxdColl = new LCCollectionVec(LCIO::SIMTRACKERHIT);
00474
00475 for(int ihit = 0; ihit < n_hit; ihit++) {
00476 LCDG4TrackerHit* g4vxd=(*VHC)[ihit];
00477
00478 // LCIO hit handling
00479 SimTrackerHitImpl* lcioVxdHit=_lcioHelper->getTrkHitLcio(*g4vxd);
00480 lcioVxdColl->push_back(lcioVxdHit);
00481
00482 double edep = g4vxd->GetEdep();
00483 float tdep = g4vxd->GetTdep();
00484 int trkid = g4vxd->GetTrackID();
00485 assert(trkid>0);
00486 G4int mctag = m_mcmgr->trkidToMctag(trkid);
00487 unsigned int layn = g4vxd->GetLayer();
00488 unsigned int sysn = g4vxd->GetSystem();
00489 unsigned int beflg = g4vxd->GetBarEnd();
00490 unsigned int nsflg = g4vxd->GetNorthSouth();
00491 G4ThreeVector xpos = g4vxd->GetPos();
00492
00493 double pos[3];
00494 pos[0] = xpos.x();
00495 pos[1] = xpos.y();
00496 pos[2] = xpos.z();
00497
00498 // pack hit info into tag
00499 unsigned int tag = layn;
00500 tag |= (sysn << 28);
00501 tag |= (nsflg << 30);
00502 tag |= (beflg << 31);
00503
00504 // Get MCParticle associated to this hit
00505 // make sure this is not a backscattering from calorimeter
00506 MCParticleSIO *mcp = NULL;
00507 int showerFlag = m_mcmgr->GetShowerFlag(trkid);
00508 if(showerFlag<=1) {
00509 mcp = mcpartmap[mctag];
00510 }
00511
00512 // save hit into hits array
00513 TkHit tk(tag,pos[0],pos[1],pos[2],edep,mcp,tdep);
00514 vxdhits.push_back(tk);
00515 } // for ihit
00516
00517 vxd = new VxdSIO("VXD",vxdhits);
00518 SIO_blockManager::add(vxd);
00519
00520 // LCIO format
00521 // set flag for long format (including position)
00522 LCFlagImpl chFlag(0) ;
00523 lcioVxdColl->setFlag( chFlag.getFlag() );
00524 _lcioHelper->saveCollection("VXDtCollection",lcioVxdColl);
00525 } // if VHC
00526
00527
00528 //************** Tracker hits *************
00529 if(THC) {
00530 int n_hit = THC->entries();
00531 dumpfile << " "<< n_hit;
00532 int ihit = 0;
00533 vector<TkHit> trkhits;
00534
00535 LCCollectionVec* lcioTrkColl = new LCCollectionVec(LCIO::SIMTRACKERHIT);
00536
00537 for(ihit = 0; ihit < n_hit; ihit++) {
00538 LCDG4TrackerHit* g4trk=(*THC)[ihit];
00539
00540 // LCIO hit handling
00541 SimTrackerHitImpl* lcioTrkHit = _lcioHelper->getTrkHitLcio(*g4trk);
00542 lcioTrkColl->push_back(lcioTrkHit);
00543
00544 double edep = g4trk->GetEdep();
00545 double tdep = g4trk->GetTdep();
00546 int trkid = g4trk->GetTrackID();
00547 assert(trkid>0);
00548 G4int mctag = m_mcmgr->trkidToMctag(trkid);
00549 int layn = g4trk->GetLayer();
00550 unsigned int sysn = g4trk->GetSystem();
00551 unsigned int beflg = g4trk->GetBarEnd();
00552 unsigned int nsflg = g4trk->GetNorthSouth();
00553 G4ThreeVector xpos = g4trk->GetPos();
00554
00555 // pack hit info into tag
00556 unsigned int tag = layn;
00557 tag |= (sysn << 28);
00558 tag |= (nsflg << 30);
00559 tag |= (beflg << 31);
00560
00561 // Get MCParticle associated to this hit
00562 // make sure this is not a backscattering from calorimeter
00563 MCParticleSIO *mcp = NULL;
00564 int showerFlag = m_mcmgr->GetShowerFlag(trkid);
00565 if(showerFlag<=1) {
00566 mcp = mcpartmap[mctag];
00567 }
00568 else {
00569 // this is a backscattering hit!
00570 // G4cout<< "backscattering hit: mcp="<< mcp << ", mctag="<< mctag
00571 // << ", trkid="<< trkid <<", pos="<< xpos << G4endl;
00572 }
00573
00574 double pos[3];
00575 pos[0] = xpos.x();
00576 pos[1] = xpos.y();
00577 pos[2] = xpos.z();
00578
00579 TkHit tkhit(tag,pos[0],pos[1],pos[2],edep,mcp,tdep);
00580 trkhits.push_back(tkhit);
00581 } // for
00582 tracker = new TrackerSIO("Tracker",trkhits);
00583 SIO_blockManager::add(tracker);
00584
00585 // LCIO format
00586 // set flag for long format (including position)
00587 LCFlagImpl chFlag(0);
00588 lcioTrkColl->setFlag( chFlag.getFlag() );
00589 _lcioHelper->saveCollection("CENtCollection",lcioTrkColl);
00590
00591 } // if THC
00592
00593 // ************* EMCAL **********************
00594 // EM Calorimeter Hits
00595 #ifdef TXT_MODE
00596 if(outf_txt) *outf_txt << "EMCAL ";
00597 #endif
00598 if(EBHC || EECHC) {
00599 // number of hits, Barrel + EndCap
00600 int n_hit =0;
00601 if(EBHC) {
00602 n_hit += EBHC->entries();
00603 }
00604 if(EECHC) {
00605 n_hit += EECHC->entries();
00606 }
00607 // write on dumpfile
00608 dumpfile << " "<< n_hit;
00609 #ifdef TXT_MODE
00610 // Write out number of hits for this event (ASCII format)
00611 if(outf_txt) *outf_txt << n_hit << endl;
00612 #endif
00613
00614 // collect barrel and endcap hits into a single collection
00615 vector<CalorimeterHit> calhits;
00616 LCCollectionVec* lcioEMColl = new LCCollectionVec(LCIO::SIMCALORIMETERHIT);
00617 if(EBHC) SaveCalorimeterHits( EBHC, calhits, lcioEMColl );
00618 if(EECHC) SaveCalorimeterHits( EECHC, calhits, lcioEMColl );
00619
00620 // save hits into sio format
00621 ecalsio = new CalorimeterSIO("EMCalorimeter",calhits);
00622 SIO_blockManager::add(ecalsio);
00623
00624 // LCIO format
00625 // set flag for long format (including position)
00626 // and PDG and cellid1
00627 LCFlagImpl chFlag(0) ;
00628 chFlag.setBit( LCIO::CHBIT_LONG );
00629 lcioEMColl->setFlag( chFlag.getFlag() );
00630 // save hits into lcio format
00631 _lcioHelper->saveCollection("EMcalCollection",lcioEMColl);
00632
00633 } // if
00634 #ifdef TXT_MODE
00635 else {
00636 if(outf_txt) *outf_txt << "0" << G4endl;
00637 }
00638 #endif
00639
00640 //************** HDCal hits *************
00641 // HAD Calorimeter Hits
00642 #ifdef TXT_MODE
00643 if(outf_txt) *outf_txt << "HCAL ";
00644 #endif
00645 if(HBHC || HECHC) {
00646 // number of hits, Barrel + EndCap
00647 int n_hit =0;
00648 if(HBHC) n_hit += HBHC->entries();
00649 if(HECHC) n_hit += HECHC->entries();
00650
00651 // write on dumpfile
00652 dumpfile << " "<< n_hit;
00653 #ifdef TXT_MODE
00654 // Write out number of hits for this event (ASCII format)
00655 if(outf_txt) *outf_txt << n_hit << endl;
00656 #endif
00657
00658 // collect barrel and endcap hits into a single collection
00659 vector<CalorimeterHit> calhits;
00660 LCCollectionVec* lcioHADColl = new LCCollectionVec(LCIO::SIMCALORIMETERHIT);
00661 if(HBHC) SaveCalorimeterHits( HBHC, calhits, lcioHADColl );
00662 if(HECHC) SaveCalorimeterHits( HECHC, calhits, lcioHADColl );
00663
00664 // save hits into sio format
00665 hcalsio = new CalorimeterSIO("HADCalorimeter",calhits);
00666 SIO_blockManager::add(hcalsio);
00667
00668 // LCIO format
00669 // set flag for long format (including position)
00670 // and PDG and cellid1
00671 LCFlagImpl chFlag(0) ;
00672 chFlag.setBit( LCIO::CHBIT_LONG );
00673 lcioHADColl->setFlag( chFlag.getFlag() );
00674 _lcioHelper->saveCollection("HADcalCollection",lcioHADColl);
00675
00676 } // if(HBHC || HECHC)
00677 #ifdef TXT_MODE
00678 else {
00679 if(outf_txt) *outf_txt << "0" << endl;
00680 }
00681 #endif
00682
00683 // ******************** LUM Cal Hits ********************
00684 if(LHC) {
00685 int n_hit = LHC->entries();
00686 dumpfile << " "<< n_hit;
00687 vector<CalorimeterHit> lumhits;
00688
00689 LCCollectionVec* lcioLumColl = new LCCollectionVec(LCIO::SIMCALORIMETERHIT);
00690
00691 int ihit = 0;
00692 for(ihit = 0; ihit < n_hit; ihit++) {
00693 LCDG4CalHit* g4lum=(*LHC)[ihit];
00694 double edeps = g4lum->GetEdep();
00695 double edepa = g4lum->GetEdepAbs();
00696
00697 int layn = g4lum->GetLayer();
00698 int sysn = g4lum->GetSystem();
00699 // int beflg = g4lum->GetBarEnd();
00700 int ithe = g4lum->GetTheta();
00701 int iphi = g4lum->GetPhi();
00702 int nmc = g4lum->GetNMC();
00703
00704 if(nmc != 0) {
00705 // LCIO hit handling
00706 SimCalorimeterHitImpl* lcioHit = _lcioHelper->getCalHitLcio(*g4lum);
00707 lcioLumColl->push_back(lcioHit);
00708
00709 // SIO hit handling
00710 LCDtowerID *towerID = new LCDtowerID(0);
00711 towerID->SetTheta(ithe);
00712 towerID->SetPhi(iphi);
00713 towerID->SetLayer(layn);
00714 towerID->SetSystem(sysn);
00715
00716 CalorimeterHit lumhit(towerID->GetTag());
00717 int imc = 0;
00718 for(imc = 0; imc < nmc; imc++) {
00719 int trkid = g4lum->GetTrkID(imc);
00720 double emc = g4lum->GetEMC(imc);
00721 G4int mctag = m_mcmgr->trkidToMctag(trkid);
00722 if(mctag != -1) {
00723 MCParticleSIO *mcp = mcpartmap[mctag];
00724 lumhit.addMcHit(mcp,emc);
00725 }
00726 else {
00727 G4cout <<" EndOfEvent LUM problem: mctag=-1 for imc,trkid="
00728 << imc <<' '<< trkid << G4endl;
00729 G4cout <<" ihit,ir,ith,iphi,nmc,edeps,abs = "
00730 << ihit <<' '<< layn <<' '<< ithe <<' '<< iphi <<' '
00731 << nmc <<' '<< edeps/MeV <<' '<< edepa/MeV << endl;
00732 }
00733 } // for nmc
00734 lumhits.push_back(lumhit);
00735 delete towerID;
00736 } // if nmc
00737 } // for hit
00738 lum = new CalorimeterSIO("LUMCalorimeter",lumhits);
00739 SIO_blockManager::add(lum);
00740
00741 // LCIO format
00742 // set flag for long format (including position)
00743 LCFlagImpl chFlag(0) ;
00744 chFlag.setBit( LCIO::CHBIT_LONG );
00745 lcioLumColl->setFlag( chFlag.getFlag() );
00746 _lcioHelper->saveCollection("LUMcalCollection",lcioLumColl);
00747
00748 } // if LUM Cal Hits
00749
00750
00751 // ******************** Muon Cal Hits ********************
00752 if(MHC) {
00753 int n_hit = MHC->entries();
00754 dumpfile << " "<< n_hit;
00755 vector<MuonHit> calhits;
00756
00757 LCCollectionVec* lcioMuonColl = new LCCollectionVec(LCIO::SIMCALORIMETERHIT);
00758 int ihit = 0;
00759 for(ihit = 0; ihit < n_hit; ihit++) {
00760 LCDG4CalHit *g4muon = (*MHC)[ihit];
00761 double edeps = g4muon->GetEdep();
00762 double edepa = g4muon->GetEdepAbs();
00763
00764 int layn = g4muon->GetLayer();
00765 // int sysn = g4muon->GetSystem();
00766 // int beflg = g4muon->GetBarEnd();
00767 int ithe = g4muon->GetTheta();
00768 int iphi = g4muon->GetPhi();
00769 int nmc = g4muon->GetNMC();
00770 if (nmc != 0) {
00771 // LCIO hit handling
00772 SimCalorimeterHitImpl* lcioMuonHit=_lcioHelper->getCalHitLcio(*g4muon);
00773 lcioMuonColl->push_back(lcioMuonHit);
00774
00775 // SIO hit handling
00776 LCDtowerID *towerID = new LCDtowerID(0);
00777 towerID->SetTheta(g4muon->GetTheta());
00778 towerID->SetPhi(g4muon->GetPhi());
00779 towerID->SetLayer(g4muon->GetLayer());
00780 towerID->SetSystem(g4muon->GetSystem());
00781
00782 MuonHit calhit(towerID->GetTag(),towerID->GetTag());
00783
00784 for(int imc = 0; imc < nmc; imc++) {
00785 int trkid = g4muon->GetTrkID(imc);
00786 double emc = g4muon->GetEMC(imc);
00787 G4int mctag = m_mcmgr->trkidToMctag(trkid);
00788 if(mctag != -1) {
00789 MCParticleSIO *mcp = mcpartmap[mctag];
00790 calhit.addMcHit(mcp,emc);
00791 }
00792 else {
00793 G4cout <<" EndOfEvent MUON problem: mctag=-1 for imc,trkid="
00794 << imc <<' '<< trkid << G4endl;
00795 G4cout <<" ihit,ir,ith,iphi,nmc,edeps,abs = "
00796 << ihit <<' '<< layn <<' '<< ithe <<' '<< iphi <<' '
00797 << nmc <<' '<< edeps/MeV <<' '<< edepa/MeV << endl;
00798 }
00799 } // for imc<nmc
00800 calhits.push_back(calhit);
00801 delete towerID;
00802 } // if nmc!=0
00803 } // for geant4 hits
00804 muon = new MuonSIO("MUCalorimeter",calhits);
00805 SIO_blockManager::add(muon);
00806
00807 // LCIO format
00808 // set flag for long format (including position)
00809 LCFlagImpl chFlag(0) ;
00810 chFlag.setBit( LCIO::CHBIT_LONG );
00811 lcioMuonColl->setFlag( chFlag.getFlag() );
00812 _lcioHelper->saveCollection("MUONcalCollection",lcioMuonColl);
00813 } // if MHC
00814
00815 dumpfile << endl;
00816 dumpfile.close();
00817
00818 // SIO
00819 stream->write(record);
00820
00821 delete vxd;
00822 delete tracker;
00823 delete ecalsio;
00824 delete hcalsio;
00825 delete lum;
00826 delete muon;
00827 delete mcprint;
00828 delete event;
00829
00830 SIO_recordManager::remove(record);
00831
00832 int evnum = evt->GetEventID();
00833 int necal=0, nhcal=0;
00834 // ECal
00835 if(EBHC||EECHC) {
00836 if(EBHC) necal += EBHC->entries();
00837 if(EECHC) necal += EECHC->entries();
00838 }
00839 // HCal
00840 if(HBHC||HECHC) {
00841 if(HBHC) nhcal += HBHC->entries();
00842 if(HECHC) nhcal += HECHC->entries();
00843 }
00844
00845 if( evnum%10 == 0 ) {
00846 G4cout<<" Event #Hits: VXD Tracker ECal HCal Muon Lum"
00847 << G4endl;
00848 }
00849 G4cout << myFormat( 5, evnum )
00850 << myFormat(13, ( VHC ? VHC->entries() : 0 ) )
00851 << myFormat( 9, ( THC ? THC->entries() : 0 ) )
00852 << myFormat( 8, necal )
00853 << myFormat( 8, nhcal )
00854 << myFormat( 6, ( MHC ? MHC->entries() : 0 ) )
00855 << myFormat( 6, ( LHC ? LHC->entries() : 0 ) )
00856 << G4endl;
00857
00858 // Write out event to LCIO output
00859 _lcioHelper->writeEvent(evt);
00860
00861 // Clear LCDG4McPartManager arrays
00862 m_mcmgr->Clear();
00863
00864 // Clear mcpartmap collection
00865 for(unsigned int i=0; i<mcpartmap.size(); ++i) {
00866 if(mcpartmap[i]) delete mcpartmap[i];
00867 }
00868 mcpartmap.clear();
00869 }
|
|
|
Definition at line 26 of file LCDG4EventAction.hh. References rndSeed.
00026 { return rndSeed;}
|
|
||||||||||||||||
|
Referenced by EndOfEventAction(). |
|
|
Definition at line 25 of file LCDG4EventAction.hh. References rndSeed. Referenced by BeginOfEventAction().
00025 { rndSeed = r; HepRandom::setTheSeed(rndSeed);}
|
|
|
Definition at line 84 of file LCDG4EventAction.cc. References _lcioHelper, LCDG4LcioHelper::getInstance(), and LCDG4LcioHelper::openLcioFile(). Referenced by main().
00085 {
00086 G4cout << "Setting up LCIO output file " << filename << G4endl;
00087 if (filename == "") {
00088 G4cerr<<"Argument for LCIO output file has not been set. Aborting.";
00089 exit(1);
00090 }
00091 _lcioHelper = LCDG4LcioHelper::getInstance();
00092 _lcioHelper->openLcioFile(string(filename));
00093 }
|
|
|
Definition at line 59 of file LCDG4EventAction.hh. Referenced by EndOfEventAction(), SetupOutputLcio(), and ~LCDG4EventAction(). |
|
|
Definition at line 54 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(). |
|
|
Definition at line 48 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), EndOfEventAction(), and LCDG4EventAction(). |
|
|
Definition at line 49 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), EndOfEventAction(), and LCDG4EventAction(). |
|
|
Definition at line 52 of file LCDG4EventAction.hh. Referenced by LCDG4EventAction(). |
|
|
Definition at line 50 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), EndOfEventAction(), and LCDG4EventAction(). |
|
|
Definition at line 51 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), EndOfEventAction(), and LCDG4EventAction(). |
|
|
Definition at line 46 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), EndOfEventAction(), and LCDG4EventAction(). |
|
|
Definition at line 56 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), and EndOfEventAction(). |
|
|
Definition at line 47 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), EndOfEventAction(), and LCDG4EventAction(). |
|
|
Definition at line 45 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), EndOfEventAction(), and LCDG4EventAction(). |
|
|
Definition at line 44 of file LCDG4EventAction.hh. Referenced by BeginOfEventAction(), EndOfEventAction(), and LCDG4EventAction(). |
|
|
Definition at line 57 of file LCDG4EventAction.hh. Referenced by EndOfEventAction(). |
|
|
Definition at line 53 of file LCDG4EventAction.hh. Referenced by GetRandSeed(), and SetRandSeed(). |
|
|
Definition at line 31 of file LCDG4EventAction.hh. Referenced by EndOfEventAction(), and LCDG4EventAction(). |
1.3.4