Geant4  10.00.p01
CCalAnalysis.cc
Go to the documentation of this file.
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 //
27 // File: CCalAnalysis.cc
28 // Description: CCalAnalysis interfaces all user analysis code
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 
44 CCalAnalysis* CCalAnalysis::instance = 0;
45 
46 CCalAnalysis::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 
147 CCalAnalysis::~CCalAnalysis() {
148  Finish();
149 }
150 
151 
152 void CCalAnalysis::Init() {
153 }
154 
155 void 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 
166 CCalAnalysis* 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
173 void 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
195 void 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
217 void 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
239 void 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
252 void 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
282 void 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 
296 void 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 */
330 void 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
352 void 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
363 void 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
void Init()
Definition: G4IonTable.cc:89
G4double z
Definition: TRTMaterials.hh:39
#define assert(x)
Definition: mymalloc.cc:1309
Definition: xmlparse.cc:179
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4double energy(const ThreeVector &p, const G4double m)
#define G4endl
Definition: G4ios.hh:61
static MCTruthManager * instance