#include "KM3SimIO.hh" #include #include #include #include #include #include #include "G4ParticleTable.hh" //#include "G4SystemOfUnits.hh" #include using namespace CLHEP; int KM3SimIO::GetNoEventsEvt(std::istream& is){ int n_events = 0; string w; cout<<"Reading number of events..."<>w; if(w == mc_keys::start_event_t && is.good()) n_events += 1; } cout<<"Number of events to process: "<print(); } else{ inputfile = new TFile(infilechar.c_str(), "READ"); input_tree = (TTree*) inputfile->Get("E"); input_tree->SetBranchAddress("Evt", &evt); header = (Head*) inputfile->Get("Head"); if(!header) header = (Head*) inputfile->Get("Header"); //header->print(); } if(EVT_output) outevt.open(outfilechar.c_str(), ofstream::out); else{ outputfile = new TFile(outfilechar.c_str(), "RECREATE"); output_tree = new TTree("E", "Evt Tree"); output_tree->Branch("Evt", &evt, 32000, 4); } //evt = outfile->pevt; if(EVT_input) nevents = GetNoEventsEvt(inevt); else{ nevents = input_tree->GetEntries(); cout<<"Number of events to process: "<size(); evt_id = 0; RunHeaderIsWrite = false; cout<<"KM3SimIO initialized!"<Write("", TObject::kOverwrite); outputfile->Close(); } if(EVT_input) inevt.close(); else inputfile->Close(); } void KM3SimIO::AddHeaderInfo(string tagname, string tagval) { /* old way to write duplicated tags! string key = tagname; for (int i = 1; header->find(key) != header->end() ; i++) { key = tagname + "_" + to_string(i); } (*header)[key] = tagval; */ header->set_line(tagname, tagval); } void KM3SimIO::WriteRunHeader() { RunHeaderIsWrite = true; if(EVT_output) write(*header, outevt); else { outputfile->WriteObject( header , "Head"); } } void KM3SimIO::ReadEvent() { cout<<"Event "<GetEntry(evt_id); G4int index = 0; G4int idbeam = 0; light_tracks_indices.clear(); for(auto &trk: evt->mc_trks){ idbeam = trk.type; if(trk.is_finalstate() && idbeam != 0 && idbeam < 1000000){ //the logic for type is kept, because it's still possible that conversion errors occur in the external software light_tracks_indices.push_back(index); //cout<<"Track id: "<Fill(); evt_id+=1; } void KM3SimIO::AddHit(int id, int PMTid, double pe, double t, int trackid, int npepure, double ttpure, int creatorProcess) { Hit hit; hit.id = id; hit.pmt_id = PMTid; hit.a = pe; hit.t = t; //hit.pure_t = ttpure; //new format does not have this //hit.pure_a = npepure; //new format does not have this hit.pattern_flags = creatorProcess; //G4cout<<"MC tracks size: "<mc_trks.size()<mc_trks[light_tracks_indices[trackid - 1]].type; track_id = evt->mc_trks[light_tracks_indices[trackid - 1]].id; // -1 comes from the fact that in Geant4 TrackID = 0 is reserved for the abstract parent track, the origin of everything } else { Gid = 6; trackid = 999999; } hit.type = Gid; hit.origin = track_id; evt->mc_trks[light_tracks_indices[trackid - 1]].hit_ids.push_back(id); // -1 comes from the fact that in Geant4 TrackID = 0 is reserved for the abstract parent track, the origin of everything evt->mc_hits.push_back(hit); } void KM3SimIO::AddTrackLength(int trackid, double track_L){ evt->mc_trks[light_tracks_indices[trackid - 1]].len = track_L / m ; // -1 comes from the fact that in Geant4 TrackID = 0 is reserved for the abstract parent track, the origin of everything } ////////////////////// // // // READER FUNCTIONS // // // ////////////////////// int KM3SimIO::GetNumberOfParticles(void) { return light_tracks_indices.size(); //this corresponds to the number of particles in the primary vertex. Incoming neutrino, the target and recoil nuclei are not taken into account in the simulation } void KM3SimIO::GetParticleInfo(int &idbeam, double &xx0, double &yy0, double &zz0, double &pxx0, double &pyy0, double &pzz0, double &t0) { G4cout<<"Do nothing for now, but maybe will be useful in the km3sim Multi Threading implementation"<mc_trks[0]; if(!neutrino.is_neutrino()) return false; idneu = neutrino.type; //if (argnumber == 15) // idtarget = (int)args[14]; //else idtarget = 1445; xneu = neutrino.pos.x * CLHEP::m / CLHEP::cm; // convert to cm yneu = neutrino.pos.y * CLHEP::m / CLHEP::cm; zneu = neutrino.pos.z * CLHEP::m / CLHEP::cm; double pmom = neutrino.E; pxneu = neutrino.dir.x * pmom; pyneu = neutrino.dir.y * pmom; pzneu = neutrino.dir.z * pmom; t0 = neutrino.t; if(verbosity>2) G4cout<<"Neutrino info: x:"<mc_trks.size()<mc_trks){ idbeam = trk.type; if( !trk.is_finalstate() || idbeam == 0 || idbeam > 1000000) continue; //if (fabs(idbeam) == 411 (D+) || fabs(idbeam) == 421 (D0) || fabs(idbeam) == 431 (D+s) || // fabs(idbeam) == 4122 (LAMBDA+c) && fabs(idbeam) == 4212 (SIGMA+c) && // fabs(idbeam) == 4222 (SIGMA++c)) // these particles are not defined or have no decay modes in GEANT4 // they are kept in the simulation, because they don't do any harm and their occurace in theorey can be counted //cout<<"Track type: "<FindParticle(idbeam); double pmass = definition->GetPDGMass() / GeV; //cout<<"PMass :"< 0.0) pmom = sqrt(totenergy * totenergy - pmass * pmass); else pmom = 0.0; //G4cout<<"Pmom is: "<SetMomentum(pxx0 * GeV, pyy0 * GeV, pzz0 * GeV); vertex->SetPrimary(particle); // Put the vertex to G4Event object anEvent->AddPrimaryVertex(vertex); } }