#include "WCSimRunAction.hh" #include "WCSimRunActionMessenger.hh" #include "G4Run.hh" #include "G4UImanager.hh" #include "G4VVisManager.hh" #include "G4ios.hh" #include "jhfNtuple.h" #ifdef REFLEX_DICTIONARY #include "Cintex/Cintex.h" #endif #include "TROOT.h" #include "TFile.h" #include "TTree.h" #include "TBranch.h" #include "TStreamerInfo.h" #include "WCSimRootEvent.hh" #include "WCSimRootGeom.hh" #include "WCSimPmtInfo.hh" #include int pawc_[500000]; // Declare the PAWC common struct ntupleStruct jhfNtuple; WCSimRunAction::WCSimRunAction(WCSimDetectorConstruction* test, WCSimRandomParameters* rand) : wcsimrandomparameters(rand) { ntuples = 1; // Messenger to allow IO options wcsimdetector = test; messenger = new WCSimRunActionMessenger(this); wcsimrootoptions = new WCSimRootOptions(); } WCSimRunAction::~WCSimRunAction() { } void WCSimRunAction::BeginOfRunAction(const G4Run* /*aRun*/) { // G4cout << "### Run " << aRun->GetRunID() << " start." << G4endl; numberOfEventsGenerated = 0; numberOfTimesWaterTubeHit = 0; numberOfTimesCatcherHit = 0; #ifdef REFLEX_DICTIONARY ROOT::Cintex::Cintex::Enable(); #endif // Needed for Root 4.00/04 WCSimRootEvent::Class()->GetStreamerInfo()->Optimize(kFALSE); // MF, aug 2006 ... you never know... WCSimRootTrigger::Class()->GetStreamerInfo()->Optimize(kFALSE); // Create the Root file // Now controlled by the messenger G4String rootname = GetRootFileName(); TFile* hfile = new TFile(rootname.c_str(),"RECREATE","WCSim ROOT file"); hfile->SetCompressionLevel(2); // Event tree WCSimTree = new TTree("wcsimT","WCSim Tree"); wcsimrootsuperevent = new WCSimRootEvent(); //empty list // wcsimrootsuperevent->AddSubEvent(); // make at least one event wcsimrootsuperevent->Initialize(); // make at least one event Int_t branchStyle = 1; //new style by default TTree::SetBranchStyle(branchStyle); Int_t bufsize = 64000; // TBranch *branch = tree->Branch("wcsimrootsuperevent", "Jhf2kmrootsuperevent", &wcsimrootsuperevent, bufsize,0); TBranch *branch = WCSimTree->Branch("wcsimrootevent", "WCSimRootEvent", &wcsimrootsuperevent, bufsize,2); // Geometry tree geoTree = new TTree("wcsimGeoT","WCSim Geometry Tree"); wcsimrootgeom = new WCSimRootGeom(); TBranch *geoBranch = geoTree->Branch("wcsimrootgeom", "WCSimRootGeom", &wcsimrootgeom, bufsize,0); FillGeoTree(); // Options tree optionsTree = new TTree("wcsimRootOptionsT","WCSim Options Tree"); optionsTree->Branch("wcsimrootoptions", "WCSimRootOptions", &wcsimrootoptions, bufsize, 0); //set detector & random options wcsimdetector->SaveOptionsToOutput(wcsimrootoptions); wcsimrandomparameters->SaveOptionsToOutput(wcsimrootoptions); } void WCSimRunAction::EndOfRunAction(const G4Run*) { //G4cout << "Number of Events Generated: "<< numberOfEventsGenerated << G4endl; //G4cout << "Number of times MRD hit: " << numberOfTimesMRDHit << G4endl; //G4cout << "Number of times FGD hit: " << numberOfTimesFGDHit << G4endl; //G4cout << "Number of times lArD hit: " << numberOfTimeslArDHit << G4endl; //G4cout<<"Number of times waterTube hit: " << numberOfTimesWaterTubeHit<Fill(); optionsTree->Write(); // Close the Root file at the end of the run TFile* hfile = WCSimTree->GetCurrentFile(); hfile->Close(); // Clean up stuff on the heap; I think deletion of hfile and trees // is taken care of by the file close delete wcsimrootsuperevent; wcsimrootsuperevent=0; delete wcsimrootgeom; wcsimrootgeom=0; } void WCSimRunAction::FillGeoTree(){ // Fill the geometry tree G4int geo_type; G4double cylinfo[3]; G4double pmtradius; G4int numpmt; G4int orientation; Float_t offset[3]; Int_t tubeNo; Float_t pos[3]; Float_t rot[3]; Int_t cylLoc; if (wcsimdetector->GetIsEggShapedHyperK()) { geo_type = 2; } else { geo_type = 0; } // geo_type 1 is for defunct mailbox design wcsimrootgeom-> SetGeo_Type(geo_type); if (geo_type == 0) { //cylinder cylinfo[1] = wcsimdetector->GetGeo_Dm(3); cylinfo[2] = wcsimdetector->GetGeo_Dm(2); wcsimrootgeom-> SetWCCylRadius(cylinfo[1]); wcsimrootgeom-> SetWCCylLength(cylinfo[2]); } pmtradius = wcsimdetector->GetPMTSize1(); numpmt = wcsimdetector->GetTotalNumPmts(); orientation = 0; wcsimrootgeom-> SetWCPMTRadius(pmtradius); wcsimrootgeom-> SetOrientation(orientation); G4ThreeVector offset1= wcsimdetector->GetWCOffset(); offset[0] = offset1[0]; offset[1] = offset1[1]; offset[2] = offset1[2]; wcsimrootgeom-> SetWCOffset(offset[0],offset[1],offset[2]); std::vector *fpmts = wcsimdetector->Get_Pmts(); WCSimPmtInfo *pmt; for (unsigned int i=0;i!=fpmts->size();i++){ pmt = ((WCSimPmtInfo*)fpmts->at(i)); pos[0] = pmt->Get_transx(); pos[1] = pmt->Get_transy(); pos[2] = pmt->Get_transz(); rot[0] = pmt->Get_orienx(); rot[1] = pmt->Get_orieny(); rot[2] = pmt->Get_orienz(); tubeNo = pmt->Get_tubeid(); cylLoc = pmt->Get_cylocation(); wcsimrootgeom-> SetPMT(i,tubeNo,cylLoc,rot,pos); } if (fpmts->size() != (unsigned int)numpmt) { G4cout << "Mismatch between number of pmts and pmt list in geofile.txt!!"<size() <<" vs. "<< numpmt < SetWCNumPMT(numpmt); geoTree->Fill(); geoTree->Write(); }