Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // 27 /// \file Par02Output.cc 28 /// \brief Implementation of the Par02Output class 29 30 #include "Par02Output.hh" 31 32 #include "Par02EventInformation.hh" 33 34 #include "G4AnalysisManager.hh" 35 #include "G4Event.hh" 36 #include "G4RunManager.hh" 37 #include "G4SystemOfUnits.hh" 38 #include "G4UnitsTable.hh" 39 40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 41 42 Par02Output* Par02Output::fPar02Output = nullptr; 43 G4ThreadLocal G4int Par02Output::fCurrentNtupleId = 0; 44 G4ThreadLocal G4int Par02Output::fCurrentID = 0; 45 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 47 48 Par02Output::Par02Output() : fFileNameWithRunNo(false) 49 { 50 fFileName = "DefaultOutput.root"; 51 } 52 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 54 55 Par02Output::~Par02Output() = default; 56 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 58 59 Par02Output* Par02Output::Instance() 60 { 61 if (!fPar02Output) { 62 fPar02Output = new Par02Output(); 63 } 64 return fPar02Output; 65 } 66 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 68 69 void Par02Output::SetFileName(G4String aName) 70 { 71 fFileName = aName; 72 } 73 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 75 76 void Par02Output::AppendName(G4bool aApp) 77 { 78 fFileNameWithRunNo = aApp; 79 } 80 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 82 83 G4String Par02Output::GetFileName() 84 { 85 return fFileName; 86 } 87 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 89 90 void Par02Output::StartAnalysis(G4int aRunID) 91 { 92 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 93 if (fFileNameWithRunNo) { 94 fFileName += "_run"; 95 fFileName += G4UIcommand::ConvertToString(aRunID); 96 } 97 analysisManager->SetDefaultFileType("root"); 98 analysisManager->SetVerboseLevel(1); 99 analysisManager->SetFileName(fFileName); 100 analysisManager->OpenFile(fFileName); 101 } 102 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 104 105 void Par02Output::EndAnalysis() 106 { 107 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 108 analysisManager->Write(); 109 analysisManager->CloseFile(); 110 } 111 112 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 113 114 void Par02Output::CreateNtuples() 115 { 116 const G4Event* event = G4RunManager::GetRunManager()->GetCurrentEvent(); 117 G4String evName = "Event_"; 118 evName += G4UIcommand::ConvertToString(event->GetEventID()); 119 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 120 fCurrentNtupleId = analysisManager->CreateNtuple(evName, evName); 121 122 analysisManager->CreateNtupleIColumn("particleID"); // column Id = 0 123 analysisManager->CreateNtupleIColumn("PID"); // column Id = 1 124 analysisManager->CreateNtupleDColumn("MC_pX"); // column Id = 2 125 analysisManager->CreateNtupleDColumn("MC_pY"); // column Id = 3 126 analysisManager->CreateNtupleDColumn("MC_pZ"); // column Id = 4 127 128 analysisManager->CreateNtupleDColumn("tracker_res"); // column Id = 5 129 analysisManager->CreateNtupleDColumn("tracker_eff"); // column Id = 6 130 analysisManager->CreateNtupleDColumn("tracker_pX"); // column Id = 7 131 analysisManager->CreateNtupleDColumn("tracker_pY"); // column Id = 8 132 analysisManager->CreateNtupleDColumn("tracker_pZ"); // column Id = 9 133 134 analysisManager->CreateNtupleDColumn("emcal_res"); // column Id = 10 135 analysisManager->CreateNtupleDColumn("emcal_eff"); // column Id = 11 136 analysisManager->CreateNtupleDColumn("emcal_X"); // column Id = 12 137 analysisManager->CreateNtupleDColumn("emcal_Y"); // column Id = 13 138 analysisManager->CreateNtupleDColumn("emcal_Z"); // column Id = 14 139 analysisManager->CreateNtupleDColumn("emcal_E"); // column Id = 15 140 141 analysisManager->CreateNtupleDColumn("hcal_res"); // column Id = 16 142 analysisManager->CreateNtupleDColumn("hcal_eff"); // column Id = 17 143 analysisManager->CreateNtupleDColumn("hcal_X"); // column Id = 18 144 analysisManager->CreateNtupleDColumn("hcal_Y"); // column Id = 19 145 analysisManager->CreateNtupleDColumn("hcal_Z"); // column Id = 20 146 analysisManager->CreateNtupleDColumn("hcal_E"); // column Id = 21 147 148 analysisManager->FinishNtuple(fCurrentNtupleId); 149 } 150 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 152 153 void Par02Output::CreateHistograms() 154 { 155 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 156 analysisManager->CreateH1("Pdiff", "momentum smeared in tracker", 100, 0.8, 1.2); 157 analysisManager->SetH1XAxisTitle(0, "p_{smeared}/p_{true}"); 158 analysisManager->SetH1YAxisTitle(0, "Entries"); 159 analysisManager->CreateH1("EMCalEdiff", "energy smeared in EMCal", 100, 0.8, 1.2); 160 analysisManager->SetH1XAxisTitle(1, "E_{smeared}/E_{true}"); 161 analysisManager->SetH1YAxisTitle(1, "Entries"); 162 analysisManager->CreateH1("HCalEdiff", "energy smeared in HCal", 100, 0.0, 2.0); 163 analysisManager->SetH1XAxisTitle(2, "E_{smeared}/E_{true}"); 164 analysisManager->SetH1YAxisTitle(2, "Entries"); 165 } 166 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 168 169 void Par02Output::SaveTrack(SaveType aWhatToSave, G4int aPartID, G4int aPDG, G4ThreeVector aVector, 170 G4double aResolution, G4double aEfficiency, G4double aEnergy) 171 { 172 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 173 switch (aWhatToSave) { 174 case Par02Output::eNoSave: 175 break; 176 case Par02Output::eSaveMC: { 177 analysisManager->FillNtupleIColumn(fCurrentNtupleId, 0, aPartID); 178 analysisManager->FillNtupleIColumn(fCurrentNtupleId, 1, aPDG); 179 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 2, aVector.x()); 180 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 3, aVector.y()); 181 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 4, aVector.z()); 182 fCurrentID = aPartID; 183 break; 184 } 185 case Par02Output::eSaveTracker: { 186 if (aPartID != fCurrentID) 187 G4cout << " Wrong particle - trying to save Tracker information of different particle" 188 << G4endl; 189 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 5, aResolution); 190 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 6, aEfficiency); 191 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 7, aVector.x()); 192 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 8, aVector.y()); 193 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 9, aVector.z()); 194 break; 195 } 196 case Par02Output::eSaveEMCal: { 197 if (aPartID != fCurrentID) 198 G4cout << " Wrong particle - trying to save EMCal information of different particle" 199 << G4endl; 200 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 10, aResolution); 201 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 11, aEfficiency); 202 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 12, aVector.x()); 203 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 13, aVector.y()); 204 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 14, aVector.z()); 205 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 15, aEnergy); 206 break; 207 } 208 case Par02Output::eSaveHCal: { 209 if (aPartID != fCurrentID) 210 G4cout << " Wrong particle - trying to save HCal information of different particle" 211 << G4endl; 212 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 16, aResolution); 213 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 17, aEfficiency); 214 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 18, aVector.x()); 215 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 19, aVector.y()); 216 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 20, aVector.z()); 217 analysisManager->FillNtupleDColumn(fCurrentNtupleId, 21, aEnergy); 218 analysisManager->AddNtupleRow(fCurrentNtupleId); 219 break; 220 } 221 } 222 } 223 224 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 225 226 void Par02Output::FillHistogram(G4int aHistNo, G4double aValue) const 227 { 228 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 229 analysisManager->FillH1(aHistNo, aValue); 230 } 231 232 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 233