// InHepEvt.cc // Contact person: Matthew Strait // See InHepEvt.hh for more details //———————————————————————// #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; namespace RAT { // The HEPEVT standard doesn't say anything about whether you can have // blank lines or not, nor does it give any method of having comments. // Nevertheless, we will tolerate empty lines, lines with only spaces // and tabs, and lines that look like they are trying to be comments. // Return true if this line is one of those cases and false otherwise. static bool IsCommentLine(const string& line) { if(line.find_first_not_of(" \t") == string::npos) return true; switch(line[0]){ case '#': case 'C': case '!': case '\'': case '"': case '/': return true; default: return false; } } void InHepEvt::ReadHepEvtFile(std::ifstream & hep) { // The global absolute time, incremented each time we read in a // particle. If the whole HEPEVT file is supposed to start at a // particular global time, that should be given as the delta-t // of the first particle, relative to RAT's default time set in // data/DATE.ratdb double gtime = 0; // These must persist across iterations because the HEPEVT standard // specifies that if a line doesn't have a entry for mass, delta-t, // position (x,y,z) or polarization (x,y,z), you use the value from // the most recent line that does. (Is this really supposed to apply // across events? It's hard to tell.) double deltat = 0; HepEvt::Particle part; memset(&part, 0, sizeof(part)); while(true){ string line; HepEvt::Event event; if(!std::getline(hep, line)) break; if(IsCommentLine(line)) continue; int nhep = -1; { stringstream ss(line); if(!(ss >> nhep)) Log::Die(Form("Bad HEPEVT line:\n%s", line.c_str())); } if(nhep < 1) Log::Die(Form("Bad HEPEVT entry with %d particles", nhep)); part.time = 0; for(int i = 0; i < nhep; i++){ if(!std::getline(hep, line)) Log::Die(Form("HEPEVT file ended unexpectedly on particle %d " "of %d in event %u", i, nhep, (unsigned)fHepEvts.size()+1)); if(IsCommentLine(line)) continue; stringstream ss(line); if(!(ss >> part.stat >> part.pdg >> part.daughter[0] >> part.daughter[1] >> part.px >> part.py >> part.pz)) Log::Die(Form("HEPEVT info missing in \"%s\"\n", line.c_str())); // Apparently until C++11, >> will leave the value alone if // conversion fails, but starting with C++11, it writes zero into // the value. Holy compatibility breaking, Batman! So be paranoid. double tmpflt; if(ss >> tmpflt) part.mass = tmpflt; if(ss >> tmpflt) deltat = tmpflt; // note special case if(ss >> tmpflt) part.x = tmpflt; if(ss >> tmpflt) part.y = tmpflt; if(ss >> tmpflt) part.z = tmpflt; if(ss >> tmpflt) part.plx = tmpflt; if(ss >> tmpflt) part.ply = tmpflt; if(ss >> tmpflt) part.plz = tmpflt; gtime += deltat; // Always set the first particle's time in an event to zero and // use the delta-t of the first particle to set the event time. if(i == 0) event.time = gtime; else part.time += deltat; event.part.push_back(part); } fHepEvts.push_back(event); info << fHepEvts.size() << " HEPEVT events read so far" << newline; } if(fHepEvts.empty()) Log::Die("No events in HEPEVT file!"); } void InHepEvt::GenerateEvent(G4Event *event) { const HepEvt::Event & evt = fHepEvts[fCurrentEvent]; std::vector g4parts; // Make all the G4PrimaryParticles out of fHepEvts for(unsigned int i = 0; i < evt.part.size(); i++){ const HepEvt::Particle & part = evt.part[i]; // HEPEVT files use GeV, but RAT uses MeV. G4PrimaryParticle * const particle = new G4PrimaryParticle( part.pdg, part.px*1000, part.py*1000, part.pz*1000); particle->SetMass(part.mass*1000); // To the best of my knowledge, polarization has no effect anywhere // in RAT, but we may as well support it here just in case. particle->SetPolarization(part.plx, part.ply, part.plz); g4parts.push_back(particle); } // Set the particle daughters, position and time, and push them into // the G4Event for(unsigned int i = 0; i < evt.part.size(); i++){ // I'm not sure if anything in RAT cares about particle daughters, // but let's set them just in case. const HepEvt::Particle & part = evt.part[i]; for(unsigned int d = 0; d < 2; d++){ const unsigned int dindex = part.daughter[d]; if(dindex == 0) continue; if(dindex > evt.part.size()){ warn << "HepEvt: Ignoring out-of-bounds daughter in event " << fCurrentEvent << newline; continue; } if(part.stat != 2) warn << "HepEvt: Non-decayed particle (status " << part.stat << ") with a daughter in event " << fCurrentEvent << "\nContinuing, but please check that." << newline; g4parts[i]->SetDaughter(g4parts[part.daughter[d]-1]); } // If the status is not 1, this particle exists only so it can have // daughters (or it is some other special case) and so we don't make // a G4PrimaryVertex with it. if(part.stat != 1) continue; G4PrimaryVertex * const vertex = new G4PrimaryVertex( G4ThreeVector(part.x, part.y, part.z), part.time); vertex->SetPrimary(g4parts[i]); event->AddPrimaryVertex(vertex); } fCurrentEvent++; if(fCurrentEvent >= fHepEvts.size()) return BecomeDone(); fNextTime = fHepEvts[fCurrentEvent].time; } void InHepEvt::BecomeDone() { G4RunManager::GetRunManager()->AbortRun(/* soft abort */ true); fDone = true; } void InHepEvt::ResetTime(double offset) { if(fCurrentEvent < fHepEvts.size()){ if(fTimeGen){ fNextTime = fTimeGen->GenerateEventTime() + offset; } else{ fNextTime = fHepEvts[fCurrentEvent].time + offset; if(fCurrentEvent > 0) fNextTime -= fHepEvts[fCurrentEvent-1].time; } } else{ fNextTime = 1e300; // Should never be consulted } } void InHepEvt::SetState(G4String filename) { ifstream infile(filename.data()); if(!infile.is_open()) Log::Die(Form("Could not open %s\n", filename.data())); ReadHepEvtFile(infile); infile.close(); } void InHepEvt::SetTimeState(G4String state) { if (fTimeGen) fTimeGen->SetState(state); else warn << "InHepEvt error: Cannot set time state, no time generator " "selected\n"; } G4String InHepEvt::GetTimeState() const { if (fTimeGen) return fTimeGen->GetState(); else return G4String("InHepEvt error: no time generator selected"); } void InHepEvt::SetPosState(G4String state) { if (fPosGen) fPosGen->SetState(state); else warn << "InHepEvt error: Cannot set position state, no position " "generator selected\n"; } G4String InHepEvt::GetPosState() const { if (fPosGen) return fPosGen->GetState(); else return G4String("InHepEvt error: no pos generator selected"); } } // namespace RAT