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

Last change on this file since 1230 was 807, checked in by garnier, 16 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.