| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 #ifdef G4ANALYSIS_USE
26 // 27 //
>> 28 // $Id: GammaRayTelAnalysis.cc,v 1.22 2010-11-10 00:15:42 asaim Exp $
>> 29 // GEANT4 tag $Name: not supported by cvs2svn $
27 // ------------------------------------------- 30 // ------------------------------------------------------------
28 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file
29 // CERN Geneva Switzerland 32 // CERN Geneva Switzerland
30 // 33 //
31 // 34 //
32 // ------------ GammaRayAnalysisManager 35 // ------------ GammaRayAnalysisManager ------
33 // by R.Giannitrapani, F.Longo & G.S 36 // by R.Giannitrapani, F.Longo & G.Santin (03 dic 2000)
34 // 37 //
35 // 03.04.2013 F.Longo/L.Pandola <<
36 // - migrated to G4tools <<
37 // <<
38 // 29.05.2003 F.Longo 38 // 29.05.2003 F.Longo
39 // - Anaphe 5.0.5 compliant << 39 // - anaphe 5.0.5 compliant
40 // 40 //
41 // 18.06.2002 R.Giannitrapani, F.Longo & G.San 41 // 18.06.2002 R.Giannitrapani, F.Longo & G.Santin
42 // - new release for Anaphe 4.0.3 42 // - new release for Anaphe 4.0.3
43 // 43 //
44 // 07.12.2001 A.Pfeiffer 44 // 07.12.2001 A.Pfeiffer
45 // - integrated Guy's addition of the ntuple 45 // - integrated Guy's addition of the ntuple
46 // 46 //
47 // 06.12.2001 A.Pfeiffer 47 // 06.12.2001 A.Pfeiffer
48 // - updating to new design (singleton) 48 // - updating to new design (singleton)
49 // 49 //
50 // 22.11.2001 G.Barrand 50 // 22.11.2001 G.Barrand
51 // - Adaptation to AIDA 51 // - Adaptation to AIDA
52 // 52 //
53 // ******************************************* 53 // ************************************************************
54 <<
55 #include <fstream> 54 #include <fstream>
56 #include <iomanip> << 55
>> 56 #include "G4RunManager.hh"
57 57
58 #include "GammaRayTelAnalysis.hh" 58 #include "GammaRayTelAnalysis.hh"
59 #include "GammaRayTelAnalysisMessenger.hh" <<
60 #include "GammaRayTelDetectorConstruction.hh" 59 #include "GammaRayTelDetectorConstruction.hh"
>> 60 #include "GammaRayTelAnalysisMessenger.hh"
61 61
62 #include "G4RunManager.hh" << 62 GammaRayTelAnalysis* GammaRayTelAnalysis::instance = 0;
63 63
>> 64 //--------------------------------------------------------------------------------
64 //....oooOO0OOooo........oooOO0OOooo........oo 65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 66
66 GammaRayTelAnalysis *GammaRayTelAnalysis::inst << 67 GammaRayTelAnalysis::GammaRayTelAnalysis()
>> 68 :GammaRayTelDetector(0),analysisFactory(0), tree(0)//, plotter(0),
>> 69 ,tuple(0)
>> 70 ,energy(0), hits(0), posXZ(0), posYZ(0)
>> 71 ,histo1DDraw("enable"),histo1DSave("enable"),histo2DDraw("enable")
>> 72 ,histo2DSave("enable"),histo2DMode("strip")
>> 73 {
>> 74 G4RunManager* runManager = G4RunManager::GetRunManager();
>> 75 GammaRayTelDetector =
>> 76 (GammaRayTelDetectorConstruction*)(runManager->GetUserDetectorConstruction());
>> 77
>> 78 #ifdef G4ANALYSIS_USE
>> 79 // Define the messenger and the analysis system
>> 80
>> 81 analysisMessenger = new GammaRayTelAnalysisMessenger(this);
>> 82
>> 83 analysisFactory = AIDA_createAnalysisFactory(); // create the Analysis Factory
>> 84 if(analysisFactory) {
>> 85 AIDA::ITreeFactory* treeFactory = analysisFactory->createTreeFactory();
>> 86 // create Tree Factory
>> 87
>> 88 if(treeFactory) {
>> 89 // Tree in memory :
>> 90 // Create a "tree" associated to an xml file
>> 91
>> 92 // tree = treeFactory->create("gammaraytel.hbook", "hbook", false, false);
>> 93 // (hbook implementation)
>> 94 tree = treeFactory->create("gammaraytel.aida","xml",false,true,"");
>> 95 if(tree) {
>> 96 // Get a tuple factory :
>> 97 ITupleFactory* tupleFactory = analysisFactory->createTupleFactory(*tree);
>> 98 if(tupleFactory) {
>> 99 // Create a tuple :
>> 100 tuple = tupleFactory->create("1","1", "float energy, plane, x, y, z");
>> 101 assert(tuple);
>> 102
>> 103 delete tupleFactory;
>> 104 }
>> 105
>> 106 IHistogramFactory* histoFactory = analysisFactory->createHistogramFactory(*tree);
67 107
68 //....oooOO0OOooo........oooOO0OOooo........oo << 108 if(histoFactory) {
>> 109 // Create histos :
69 110
70 GammaRayTelAnalysis::GammaRayTelAnalysis() : d << 111 int Nplane = GammaRayTelDetector->GetNbOfTKRLayers();
71 detector = static_cast<const GammaRayTelDe << 112 int Nstrip = GammaRayTelDetector->GetNbOfTKRStrips();
>> 113 int Ntile = GammaRayTelDetector->GetNbOfTKRTiles();
>> 114 double sizexy = GammaRayTelDetector->GetTKRSizeXY();
>> 115 double sizez = GammaRayTelDetector->GetTKRSizeZ();
>> 116 int N = Nstrip*Ntile;
>> 117
>> 118 // 1D histogram that store the energy deposition of the
>> 119 // particle in the last (number 0) TKR X-plane
>> 120
>> 121 energy = histoFactory->createHistogram1D("10","Edep in the last X plane (keV)", 100, 50, 200);
>> 122
>> 123 // 1D histogram that store the hits distribution along the TKR X-planes
>> 124
>> 125 hits = histoFactory->createHistogram1D("20","Hits dist in TKR X planes",Nplane, 0, Nplane-1);
>> 126 // 2D histogram that store the position (mm) of the hits (XZ projection)
>> 127
>> 128 if (histo2DMode == "strip"){
>> 129 posXZ = histoFactory->createHistogram2D("30","Tracker Hits XZ (strip,plane)",
>> 130 N, 0, N-1,
>> 131 2*Nplane, 0, Nplane-1);
>> 132 }
>> 133 else
>> 134 {
>> 135 posXZ = histoFactory->createHistogram2D("30","Tracker Hits XZ (x,z) in mm",
>> 136 int(sizexy/5), -sizexy/2, sizexy/2,
>> 137 int(sizez/5), -sizez/2, sizez/2);
>> 138 }
>> 139 // 2D histogram that store the position (mm) of the hits (YZ projection)
>> 140
>> 141 if(histo2DMode=="strip")
>> 142 posYZ = histoFactory->createHistogram2D("40","Tracker Hits YZ (strip,plane)",
>> 143 N, 0, N-1,
>> 144 2*Nplane, 0, Nplane-1);
>> 145 else
>> 146 posYZ = histoFactory->createHistogram2D("40","Tracker Hits YZ (y,z) in mm",
>> 147 int(sizexy/5), -sizexy/2, sizexy/2,
>> 148 int(sizez/5), -sizez/2, sizez/2);
>> 149
>> 150 delete histoFactory;
>> 151 }
>> 152
>> 153 }
>> 154 }
>> 155 delete treeFactory; // Will not delete the ITree.
>> 156
>> 157 }
>> 158
>> 159 // IPlotterFactory* plotterFactory = analysisFactory->createPlotterFactory(0,0);
>> 160 // if(plotterFactory) {
>> 161 // plotter = plotterFactory->create();
>> 162 // if(plotter) {
>> 163 // plotter->show();
>> 164 // plotter->setParameter("pageTitle","Gamma Ray Tel");
>> 165
>> 166 // }
>> 167 // delete plotterFactory;
>> 168 // }
72 169
73 // Define the messenger and the analysis s << 170 #endif
74 analysisMessenger = new GammaRayTelAnalysi <<
75 histogramFileName = "gammaraytel"; <<
76 } 171 }
77 172
78 //....oooOO0OOooo........oooOO0OOooo........oo <<
79 173
80 GammaRayTelAnalysis::~GammaRayTelAnalysis() { 174 GammaRayTelAnalysis::~GammaRayTelAnalysis() {
81 Finish(); << 175 Finish();
82 } 176 }
83 177
84 //....oooOO0OOooo........oooOO0OOooo........oo <<
85 <<
86 void GammaRayTelAnalysis::Init() { <<
87 } <<
88 <<
89 //....oooOO0OOooo........oooOO0OOooo........oo <<
90 178
91 void GammaRayTelAnalysis::Finish() { << 179 void GammaRayTelAnalysis::Init()
92 delete analysisMessenger; << 180 {
93 analysisMessenger = nullptr; << 181 }
94 } << 182
95 << 183 void GammaRayTelAnalysis::Finish()
96 //....oooOO0OOooo........oooOO0OOooo........oo << 184 {
97 << 185 #ifdef G4ANALYSIS_USE
98 auto GammaRayTelAnalysis::getInstance() -> Gam << 186 delete tree;
99 if (instance == nullptr) { << 187 //delete plotter;
100 instance = new GammaRayTelAnalysis(); << 188 // delete analysisFactory; // Will delete tree and histos.
101 } << 189 delete analysisMessenger;
102 return instance; << 190
>> 191 analysisMessenger = 0;
>> 192 #endif
>> 193 }
>> 194
>> 195 GammaRayTelAnalysis* GammaRayTelAnalysis::getInstance()
>> 196 {
>> 197 if (instance == 0) instance = new GammaRayTelAnalysis();
>> 198 return instance;
103 } 199 }
104 200
105 //....oooOO0OOooo........oooOO0OOooo........oo <<
106 <<
107 // This function fill the 2d histogram of the 201 // This function fill the 2d histogram of the XZ positions
108 void GammaRayTelAnalysis::InsertPositionXZ(G4d << 202 void GammaRayTelAnalysis::InsertPositionXZ(double x, double z)
109 auto *manager = G4AnalysisManager::Instance( << 203 {
110 manager->FillH2(1, x, z); << 204 #ifdef G4ANALYSIS_USE
>> 205 if(posXZ) posXZ->fill(x, z);
>> 206 #endif
111 } 207 }
112 208
113 //....oooOO0OOooo........oooOO0OOooo........oo <<
114 <<
115 // This function fill the 2d histogram of the 209 // This function fill the 2d histogram of the YZ positions
116 void GammaRayTelAnalysis::InsertPositionYZ(G4d << 210 void GammaRayTelAnalysis::InsertPositionYZ(double y, double z)
117 auto *manager = G4AnalysisManager::Instance( << 211 {
118 manager->FillH2(2, y, z); << 212 #ifdef G4ANALYSIS_USE
>> 213 if(posYZ) posYZ->fill(y, z);
>> 214 #endif
119 } 215 }
120 216
121 //....oooOO0OOooo........oooOO0OOooo........oo <<
122 <<
123 // This function fill the 1d histogram of the 217 // This function fill the 1d histogram of the energy released in the last Si plane
124 void GammaRayTelAnalysis::InsertEnergy(G4doubl << 218 void GammaRayTelAnalysis::InsertEnergy(double en)
125 auto *manager = G4AnalysisManager::Instance( << 219 {
126 manager->FillH1(1, energy); << 220 #ifdef G4ANALYSIS_USE
>> 221 if(energy) energy->fill(en);
>> 222 #endif
127 } 223 }
128 224
129 //....oooOO0OOooo........oooOO0OOooo........oo <<
130 <<
131 // This function fill the 1d histogram of the 225 // This function fill the 1d histogram of the hits distribution along the TKR planes
132 void GammaRayTelAnalysis::InsertHits(G4int pla << 226 void GammaRayTelAnalysis::InsertHits(int nplane)
133 auto *manager = G4AnalysisManager::Instance( << 227 {
134 manager->FillH1(2, planeNumber); << 228 #ifdef G4ANALYSIS_USE
>> 229 if(hits) hits->fill(nplane);
>> 230 #endif
>> 231 }
>> 232
>> 233 void GammaRayTelAnalysis::setNtuple(float E, float p, float x, float y, float z)
>> 234 {
>> 235 #ifdef G4ANALYSIS_USE
>> 236 tuple->fill(tuple->findColumn("energy"),E);
>> 237 tuple->fill(tuple->findColumn("plane"),p);
>> 238 tuple->fill(tuple->findColumn("x"),x);
>> 239 tuple->fill(tuple->findColumn("y"),y);
>> 240 tuple->fill(tuple->findColumn("z"),z);
>> 241 tuple->addRow();
>> 242 #endif
135 } 243 }
136 244
137 //....oooOO0OOooo........oooOO0OOooo........oo 245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
>> 246 /*
>> 247 This member reset the histograms and it is called at the begin
>> 248 of each run; here we put the inizialization so that the histograms have
>> 249 always the right dimensions depending from the detector geometry
>> 250 */
>> 251
>> 252 //void GammaRayTelAnalysis::BeginOfRun(G4int n) // to be reintroduced
>> 253 void GammaRayTelAnalysis::BeginOfRun()
>> 254 {
>> 255 #ifdef G4ANALYSIS_USE
>> 256
>> 257 if(energy) energy->reset();
>> 258 if(hits) hits->reset();
>> 259 if(posXZ) posXZ->reset();
>> 260 if(posYZ) posYZ->reset();
138 261
139 void GammaRayTelAnalysis::setNtuple(G4double e << 262 #endif
140 auto *manager = G4AnalysisManager::Instance( <<
141 manager->FillNtupleDColumn(0, energy); <<
142 manager->FillNtupleDColumn(1, planeNumber); <<
143 manager->FillNtupleDColumn(2, x); <<
144 manager->FillNtupleDColumn(3, y); <<
145 manager->FillNtupleDColumn(4, z); <<
146 manager->AddNtupleRow(); <<
147 } 263 }
148 264
149 //....oooOO0OOooo........oooOO0OOooo........oo 265 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
150 266
151 /* 267 /*
152 This member reset the histograms and it is ca << 268 This member is called at the end of each run
153 here we put the inizialization so that the hi << 269 */
154 always the right dimensions depending from th << 270 void GammaRayTelAnalysis::EndOfRun(G4int)
155 */ << 271 {
156 void GammaRayTelAnalysis::BeginOfRun() { << 272 #ifdef G4ANALYSIS_USE
157 auto *manager = G4AnalysisManager::Instance( << 273 if(tree) {
158 manager->SetDefaultFileType("root"); << 274 tree->commit();
159 << 275 tree->close();
160 // Open an output file << 276 }
161 << 277
162 G4cout << "Opening output file " << histogra << 278 // if(plotter) {
163 manager->OpenFile(histogramFileName); << 279
164 manager->SetFirstHistoId(1); << 280 // // We set one single region for the plotter
165 G4cout << " done" << G4endl; << 281 // // We now print the histograms, each one in a separate file
166 << 282
167 auto Nplane = detector->GetNbOfTKRLayers(); << 283 // if(histo2DSave == "enable") {
168 auto numberOfStrips = detector->GetNbOfTKRSt << 284 // char name[15];
169 auto numberOfTiles = detector->GetNbOfTKRTil << 285 // plotter->createRegions(1,1);
170 auto sizeXY = detector->GetTKRSizeXY(); << 286 // sprintf(name,"posxz_%d.ps", n);
171 auto sizeZ = detector->GetTKRSizeZ(); << 287 // plotter->currentRegion().plot(*posXZ);
172 auto N = numberOfStrips * numberOfTiles; << 288 // plotter->refresh();
173 << 289 // // plotter->write(name,"ps"); // temporary unavailable
174 // Book 1D histograms << 290
175 //------------------- << 291 // plotter->createRegions(1,1);
176 << 292 // sprintf(name,"posyz_%d.ps", n);
177 constexpr auto NUMBER_OF_BINS{100}; << 293 // plotter->currentRegion().plot(*posYZ);
178 constexpr auto LOWER_BOUND{50}; << 294 // plotter->next().plot(*posYZ);
179 constexpr auto UPPER_BOUND{200}; << 295 // plotter->refresh();
180 << 296 // // plotter->write(name,"ps"); // temporary unavailable
181 // 1D histogram that store the energy deposi << 297 // }
182 manager->CreateH1("1", "Deposited energy in << 298
183 << 299 // if(histo1DSave == "enable") {
184 // 1D histogram that store the hits distribu << 300 // plotter->createRegions(1,1);
185 manager->CreateH1("2", "Hits distribution in << 301 // char name[15];
186 << 302 // sprintf(name,"energy_%d.ps", n);
187 // Book 2D histograms << 303 // plotter->currentRegion().plot(*energy);
188 //------------------- << 304 // plotter->refresh();
189 << 305 // // plotter->write(name,"ps"); // temporary unavailable
190 // 2D histogram that store the position (mm) << 306 // plotter->createRegions(1,1);
191 << 307 // sprintf(name,"hits_%d.ps", n);
192 if (histo2DMode == "strip") { << 308 // plotter->currentRegion().plot(*hits);
193 manager->CreateH2("1", "Tracker hits X << 309 // plotter->refresh();
194 } else { << 310 // // plotter->write(name,"ps"); // temporary unavailable
195 manager->CreateH2("1", "Tracker hits X << 311 // plotter->createRegions(1,2);
196 } << 312 // plotter->currentRegion().plot(*energy);
>> 313 // plotter->next().plot(*hits);
>> 314 // plotter->refresh();
>> 315 // }
>> 316 #endif
>> 317 }
>> 318
>> 319 /* This member is called at the end of every event */
>> 320 void GammaRayTelAnalysis::EndOfEvent(G4int flag)
>> 321 {
>> 322 // The plotter is updated only if there is some
>> 323 // hits in the event
>> 324 if(!flag) return;
>> 325 #ifdef G4ANALYSIS_USE
>> 326 // Set the plotter ; set the number of regions and attach histograms
>> 327 // to plot for each region.
>> 328 // It is done here, since then EndOfRun set regions
>> 329 // for paper output.
>> 330 // if(plotter) {
>> 331 // if((histo2DDraw == "enable") && (histo1DDraw == "enable")) {
>> 332 // plotter->createRegions(1,2);
>> 333 // //plotter->currentRegion().plot(*posXZ); //temporary unavailable
>> 334 // plotter->currentRegion().plot(*hits);
>> 335 // // plotter->next().plot(*posYZ); //temporary unavailable
>> 336 // plotter->next().plot(*energy);
>> 337 // //plotter->next().plot(*energy);
>> 338 // // plotter->currentRegion().plot(*hits);
>> 339 // //plotter->next().plot(*hits);
>> 340 // } else if((histo1DDraw == "enable") && (histo2DDraw != "enable")) {
>> 341 // plotter->createRegions(1,2);
>> 342 // plotter->currentRegion().plot(*energy);
>> 343 // plotter->next().plot(*hits);
>> 344 // } else if((histo1DDraw != "enable") && (histo2DDraw == "enable")) {
>> 345 // /* plotter->createRegions(1,2);
>> 346 // plotter->currentRegion().plot(*posXZ);
>> 347 // plotter->next().plot(*posYZ);*/
>> 348 // G4cout << "Temporary Unavailable " << G4endl;
>> 349 // } else { // Nothing to plot.
>> 350 // plotter->createRegions(1,1);
>> 351 // }
>> 352 // plotter->refresh();
>> 353 // }
197 354
198 // 2D histogram that store the position (m << 355 #endif
199 if (histo2DMode == "strip") { <<
200 manager->CreateH2("2", "Tracker hits Y <<
201 } else { <<
202 manager->CreateH2("2", "Tracker hits Y <<
203 } <<
204 <<
205 // Book n-tuple (energy, plane, x, y, z) <<
206 //------------------------------------------ <<
207 manager->CreateNtuple("1", "Track n-tuple"); <<
208 manager->CreateNtupleDColumn("energy"); <<
209 manager->CreateNtupleDColumn("plane"); <<
210 manager->CreateNtupleDColumn("x"); <<
211 manager->CreateNtupleDColumn("y"); <<
212 manager->CreateNtupleDColumn("z"); <<
213 manager->FinishNtuple(); <<
214 } 356 }
>> 357 #endif
>> 358
>> 359
>> 360
>> 361
215 362
216 /* <<
217 This member is called at the end of each run <<
218 */ <<
219 void GammaRayTelAnalysis::EndOfRun() { <<
220 // Save histograms <<
221 auto *manager = G4AnalysisManager::Instance( <<
222 manager->Write(); <<
223 manager->CloseFile(); <<
224 } <<
225 363
226 //....oooOO0OOooo........oooOO0OOooo........oo <<
227 364
228 // This member is called at the end of every e <<
229 void GammaRayTelAnalysis::EndOfEvent(G4int /* <<
230 } <<
231 365