source: trunk/examples/advanced/composite_calorimeter/src/CCalAnalysis.cc@ 1230

Last change on this file since 1230 was 807, checked in by garnier, 17 years ago

update
File size: 12.8 KB
Line 
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: CCalAnalysis.cc
28// Description: CCalAnalysis interfaces all user analysis code
29///////////////////////////////////////////////////////////////////////////////
30
31#ifdef G4ANALYSIS_USE
32
33#include "G4RunManager.hh"
34
35#include "CCalAnalysis.hh"
36#include "CCalutils.hh"
37
38#include <AIDA/AIDA.h>
39
40#include <typeinfo>
41
42//#define debug
43
44CCalAnalysis* CCalAnalysis::instance = 0;
45
46CCalAnalysis::CCalAnalysis() :analysisFactory(0), tree(0), tuple(0), energy(0) {
47
48 for (int i=0; i<28; i++) {hcalE[i] = 0;}
49 for (int i=0; i<70; i++) {lateralProfile[i] = 0;}
50 for (int i=0; i<49; i++) {ecalE[i] = 0;}
51 for (int i=0; i<numberOfTimeSlices; i++) {timeHist[i] = 0;}
52 for (int i=0; i<2; i++) {timeProfile[i] = 0;}
53
54 analysisFactory = AIDA_createAnalysisFactory();
55 if (analysisFactory) {
56
57 AIDA::ITreeFactory* treeFactory = analysisFactory->createTreeFactory();
58 if (treeFactory) {
59 // Tree in memory :
60 // Create a "tree" associated to an hbook
61 const char* opFileptr = getenv("CCAL_FILENAME");
62 G4String opFilestr = "ccal.his";
63 if (opFileptr) opFilestr = opFileptr;
64 G4cout << "********************************************" << G4endl
65 << "* o/p file on " << opFilestr << G4endl
66 << "********************************************" << G4endl
67 << G4endl;
68 bool readOnly = false; // we want to write.
69 bool createNew = true; // create file if it doesn't exist.
70 tree = treeFactory->create(opFilestr, "hbook", readOnly,createNew);
71 if (tree) {
72 // Get a tuple factory :
73 AIDA::ITupleFactory* tupleFactory = analysisFactory->createTupleFactory(*tree);
74 if (tupleFactory) {
75 // Create a tuple :
76 G4String tag2, tag = "float";
77 for (int i=0; i<28; i++) {
78 tag2 = tag + " hcal" + i + ",";
79 tag = tag2;
80 }
81 for (int i=0; i<49; i++) {
82 tag2 = tag + " ecal" + i + ",";
83 tag = tag2;
84 }
85 tag2 = tag + " ELAB, XPOS, YPOS, ZPOS";
86 tag = tag2 + ", EDEP, EDEC, EDHC";
87
88 tuple = tupleFactory->create("1","Event info", tag); // Column wise (default)
89 //tuple = tupleFactory->create("1","Event info", tag, "--preferRWN"); // Row wise
90
91 assert(tuple);
92
93 delete tupleFactory;
94 }
95
96 AIDA::IHistogramFactory* histoFactory =
97 analysisFactory->createHistogramFactory(*tree);
98 if (histoFactory) {
99 // Create histos :
100 char id[4], ntupletag[50];
101 //Energy deposit in Hcal layers
102 for (int i = 0; i<28; i++) {
103 sprintf(id, "%d",i+100);
104 sprintf(ntupletag, "Energy Deposit in Hcal Layer%d in GeV",i);
105 hcalE[i] = histoFactory->createHistogram1D(id, ntupletag, 100, 0., 1.0);
106 }
107 // Energy deposits in Ecal towers
108 for (int i = 0; i<49; i++) {
109 sprintf(id, "%d",i+200);
110 sprintf(ntupletag, "Energy Deposit in Ecal Tower%d in GeV",i);
111 ecalE[i] = histoFactory->createHistogram1D(id, ntupletag, 100, 0., 100.0);
112 }
113 // Total energy deposit
114 energy = histoFactory->createHistogram1D("4000", "Total energy deposited in GeV",
115 100, 0., 100.0);
116
117 // Time slices
118 for (int i=0; i<numberOfTimeSlices; i++){
119 sprintf(id, "%d",i+300);
120 sprintf(ntupletag, "Time slice %d nsec energy profile in GeV",i);
121 timeHist[i] = histoFactory->createHistogram1D(id, ntupletag, 100, 0., 100.0);
122 }
123
124 // Profile of lateral energy deposit in Hcal
125 for (int i = 0; i<70; i++) {
126 sprintf(id, "%d",i+500);
127 sprintf(ntupletag, "Lateral energy profile at %d cm in GeV",i);
128 lateralProfile[i] = histoFactory->createHistogram1D(id, ntupletag, 100, 0., 10.0);
129 }
130
131 // Time profile
132 timeProfile[0] = histoFactory->createHistogram1D("901", "Time Profile in Sensitive Detector", 200, 0., 200.);
133 timeProfile[1] = histoFactory->createHistogram1D("902", "Time Profile in Sensitive+Passive", 200, 0., 200.);
134
135 delete histoFactory;
136 }
137
138 }
139 delete treeFactory; // Will not delete the ITree.
140 }
141
142 }
143
144}
145
146
147CCalAnalysis::~CCalAnalysis() {
148 Finish();
149}
150
151
152void CCalAnalysis::Init() {
153}
154
155void CCalAnalysis::Finish() {
156 if (tree) {
157 delete tree;
158 tree = 0;
159 }
160 if (analysisFactory) {
161 delete analysisFactory; // Will delete tree and histos.
162 analysisFactory = 0;
163 }
164}
165
166CCalAnalysis* CCalAnalysis::getInstance() {
167 if (instance == 0) instance = new CCalAnalysis();
168 return instance;
169}
170
171
172// This function fill the 1d histogram of the energies in HCal layers
173void CCalAnalysis::InsertEnergyHcal(float* v) {
174#ifdef debug
175 double totalFilledEnergyHcal = 0.0;
176#endif
177 for (int i=0; i<28; i++) {
178 if (hcalE[i]) {
179 double x = v[i];
180 hcalE[i]->fill(x);
181#ifdef debug
182 G4cout << "Fill Hcal histo " << i << " with " << x << G4endl;
183 totalFilledEnergyHcal += x;
184#endif
185 }
186 }
187#ifdef debug
188 G4cout << "CCalAnalysis::InsertEnergyHcal: Total filled Energy Hcal histo "
189 << totalFilledEnergyHcal << G4endl;
190#endif
191}
192
193
194// This function fill the 1d histogram of the energies in ECal layers
195void CCalAnalysis::InsertEnergyEcal(float* v) {
196#ifdef debug
197 double totalFilledEnergyEcal = 0.0;
198#endif
199 for (int i=0; i<49; i++) {
200 if (ecalE[i]) {
201 double x = v[i];
202 ecalE[i]->fill(x);
203#ifdef debug
204 G4cout << "Fill Ecal histo " << i << " with " << x << G4endl;
205 totalFilledEnergyEcal += x;
206#endif
207 }
208 }
209#ifdef debug
210 G4cout << "CCalAnalysis::InsertEnergyEcal: Total filled Energy Ecal histo "
211 << totalFilledEnergyEcal << G4endl;
212#endif
213}
214
215
216// This function fill the 1d histogram of the lateral profile
217void CCalAnalysis::InsertLateralProfile(float* v) {
218#ifdef debug
219 double totalFilledProfileHcal = 0.0;
220#endif
221 for (int i=0; i<70; i++) {
222 if (lateralProfile[i]) {
223 double x = v[i];
224 lateralProfile[i]->fill(x);
225#ifdef debug
226 G4cout << "Fill Profile Hcal histo " << i << " with " << x << G4endl;
227 totalFilledProfileHcal += x;
228#endif
229 }
230 }
231#ifdef debug
232 G4cout << "CCalAnalysis::InsertLateralProfile: Total filled Profile Hcal"
233 << " histo " << totalFilledProfileHcal << G4endl;
234#endif
235}
236
237
238// This function fill the 1d histogram of the energy
239void CCalAnalysis::InsertEnergy(float v) {
240 if (energy) {
241 double x = v;
242 energy->fill(x);
243#ifdef debug
244 G4cout << "CCalAnalysis::InsertEnergy: Fill Total energy Hcal histo with "
245 << x << G4endl;
246#endif
247 }
248}
249
250
251// This function fill the 1d histograms of time profiles at stepping action
252void CCalAnalysis::InsertTime(float* v) {
253#ifdef debug
254 double totalFilledTimeProfile = 0.0;
255#endif
256 for (int j=0; j<numberOfTimeSlices; j++) {
257 if (timeHist[j]) {
258 double x = v[j];
259 timeHist[j]->fill(x);
260#ifdef debug
261 G4cout << "Fill Time slice histo " << j << " with " << x << G4endl;
262 totalFilledTimeProfile += x;
263#endif
264 }
265 if (timeProfile[1]) {
266 double x = v[j];
267 double t = j + 0.5;
268 timeProfile[1]->fill(t,x);
269#ifdef debug
270 G4cout << "Fill Time profile histo 1 with " << t << " " << x << G4endl;
271#endif
272 }
273 }
274#ifdef debug
275 G4cout << "CCalAnalysis::InsertTime: Total filled Time profile histo "
276 << totalFilledTimeProfile << G4endl;
277#endif
278}
279
280
281// This function fill the 1d histograms of time profiles in SD
282void CCalAnalysis::InsertTimeProfile(int hit, double time, double edep) {
283
284 if (timeProfile[0]) {
285 timeProfile[0]->fill(time,edep);
286#ifdef debug
287 G4cout << "CCalAnalysis:: Fill Time Profile with Hit " << hit
288 << " Edeposit " << edep << " Gev at " << time << " ns" << G4endl;
289#else
290 hit=0; // Just to avoid compiler warning!
291#endif
292 }
293}
294
295
296void CCalAnalysis::setNtuple(float* hcalE, float* ecalE, float elab,
297 float x, float y, float z, float edep,
298 float edec, float edhc) {
299
300 if (tuple) {
301 char tag[10];
302 for (int i=0; i<28; i++) {
303 sprintf (tag, "hcal%d", i);
304 tuple->fill(tuple->findColumn(tag),hcalE[i]);
305 }
306 for (int i=0; i<49; i++) {
307 sprintf (tag, "ecal%d", i);
308 tuple->fill(tuple->findColumn(tag),ecalE[i]);
309 }
310 tuple->fill(tuple->findColumn("ELAB"),elab);
311 tuple->fill(tuple->findColumn("XPOS"),x);
312 tuple->fill(tuple->findColumn("YPOS"),y);
313 tuple->fill(tuple->findColumn("ZPOS"),z);
314 tuple->fill(tuple->findColumn("EDEP"),edep);
315 tuple->fill(tuple->findColumn("EDEC"),edec);
316 tuple->fill(tuple->findColumn("EDHC"),edhc);
317 tuple->addRow();
318#ifdef debug
319 G4cout << "CCalAnalysis:: Fill Ntuple " << G4endl;
320#endif
321 }
322}
323
324
325/*
326 This member reset the histograms and it is called at the begin
327 of each run; here we put the inizialization so that the histograms have
328 always the right dimensions depending from the detector geometry
329*/
330void CCalAnalysis::BeginOfRun(G4int ) {
331
332 /*
333 if (energy) energy->reset();
334 for (int i=0; i<28; i++) {
335 if (hcalE[i]) hcalE[i]->reset();
336 }
337 for (int i=0; i<70; i++) {
338 if (lateralProfile[i]) lateralProfile[i]->reset();
339 }
340 for (int i=0; i<49; i++) {
341 if (ecalE[i]) ecalE[i]->reset();
342 }
343 */
344 for (int i=0; i<numberOfTimeSlices; i++) {
345 if (timeHist[i]) timeHist[i]->reset();
346 }
347
348}
349
350
351// This member is called at the end of each run
352void CCalAnalysis::EndOfRun(G4int ) {
353
354 if (tree) {
355 tree->commit();
356 tree->close();
357 }
358
359}
360
361
362// This member is called at the end of every event
363void CCalAnalysis::EndOfEvent(G4int flag) {
364
365#ifdef debug
366 G4cout << " Check if empty histograms " << G4endl;
367 if ( energy ) {
368 if ( energy->allEntries() == 0 ) {
369 G4cout << "EMPTY HISTO energy " << G4endl;
370 } else if ( energy->allEntries() == energy->extraEntries() ) {
371 G4cout << "EXTRA entries only HISTO energy " << G4endl;
372 }
373 } else {
374 G4cout << "UNDEFINED HISTO energy " << G4endl;
375 }
376 for (int i=0; i<28; i++) {
377 if ( hcalE[i] ) {
378 if ( hcalE[i]->allEntries() == 0 ) {
379 G4cout << "EMPTY HISTO hcal " << i << G4endl;
380 } else if ( hcalE[i]->allEntries() == hcalE[i]->extraEntries() ) {
381 G4cout << "EXTRA entries only HISTO hcal " << i << G4endl;
382 }
383 } else {
384 G4cout << "UNDEFINED HISTO hcal " << i << G4endl;
385 }
386 }
387 for (int i=0; i<70; i++) {
388 if ( lateralProfile[i] ) {
389 if ( lateralProfile[i]->allEntries() == 0 ) {
390 G4cout << "EMPTY HISTO lateralProfile " << i << G4endl;
391 } else if ( lateralProfile[i]->allEntries() ==
392 lateralProfile[i]->extraEntries() ) {
393 G4cout << "EXTRA entries only HISTO lateralProfile " << i << G4endl;
394 }
395 } else {
396 G4cout << "UNDEFINED HISTO lateralProfile " << i << G4endl;
397 }
398 }
399 for (int i=0; i<49; i++) {
400 if ( ecalE[i] ) {
401 if ( ecalE[i]->allEntries() == 0 ) {
402 G4cout << "EMPTY HISTO ecalE " << i << G4endl;
403 } else if ( ecalE[i]->allEntries() == ecalE[i]->extraEntries() ) {
404 G4cout << "EXTRA entries only HISTO ecal " << i << G4endl;
405 }
406 } else {
407 G4cout << "UNDEFINED HISTO hcal " << i << G4endl;
408 }
409 }
410 for (int i=0; i<numberOfTimeSlices; i++) {
411 if ( timeHist[i] ) {
412 if ( timeHist[i]->allEntries() == 0 ) {
413 G4cout << "EMPTY HISTO timeHist " << i << G4endl;
414 } else if ( timeHist[i]->allEntries() == timeHist[i]->extraEntries() ) {
415 G4cout << "EXTRA entries only HISTO timeHist " << i << G4endl;
416 }
417 } else {
418 G4cout << "UNDEFINED HISTO timeHist " << i << G4endl;
419 }
420 }
421#endif
422
423 // The plotter is updated only if there is some
424 // hits in the event
425 if (!flag) return;
426}
427
428#endif
Note: See TracBrowser for help on using the repository browser.