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