Geant4  10.01
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 
32 #include "G4RunManager.hh"
33 
34 #include "CCalAnalysis.hh"
35 #include "CCalutils.hh"
36 
38 
40  energy(0), hcalE(0), ecalE(0), timeHist(0), lateralProfile(0),
41  timeProfile(0)
42 {
43 
44 #ifdef debug
45  fVerbosity = 1;
46 #else
47  fVerbosity = 0;
48 #endif
49 
50  numberOfTimeSlices = 200;
51 
52  // Create analysis manager
53  G4AnalysisManager* analysisManager = G4AnalysisManager::Instance();
54  analysisManager->SetVerboseLevel(1);
55  analysisManager->SetFirstHistoId(1);
56  analysisManager->SetFirstNtupleId(1);
57 
58  // Open an output file
59  analysisManager->OpenFile("ccal");
60  G4cout << "********************************************" << G4endl
61  << "* o/p file ccal" << G4endl
62  << "********************************************" << G4endl
63  << G4endl;
64 
65 
66  // Create a tuple :
67  // Create ntuple
68  analysisManager->CreateNtuple("ntuple1", "Event info");
69  char tupleid[6];
70  for (int i=0;i<28;i++)
71  {
72  sprintf(tupleid,"hcal%d",i);
73  analysisManager->CreateNtupleFColumn(tupleid);
74  }
75  for (int i=0; i<49; i++)
76  {
77  sprintf(tupleid,"ecal%d",i);
78  analysisManager->CreateNtupleFColumn(tupleid);
79  }
80  analysisManager->CreateNtupleFColumn("ELAB");
81  analysisManager->CreateNtupleFColumn("XPOS");
82  analysisManager->CreateNtupleFColumn("YPOS");
83  analysisManager->CreateNtupleFColumn("ZPOS");
84  analysisManager->CreateNtupleFColumn("EDEP");
85  analysisManager->CreateNtupleFColumn("EDEC");
86  analysisManager->CreateNtupleFColumn("EHDC");
87  analysisManager->FinishNtuple();
88 
89  //
90  //Create histograms. Save the ID of the first histogram in each block
91  //
92  char id[5], ntupletag[50];
93  //Energy deposit in Hcal layers
94  for (int i = 0; i<28; i++) {
95  sprintf(id,"h%d",i+100);
96  sprintf(ntupletag, "Energy Deposit in Hcal Layer%d in GeV",i);
97  G4int histoID = analysisManager->CreateH1(id,ntupletag, 100, 0., 1.0);
98  if (!i)
99  hcalE = histoID;
100  }
101  // Energy deposits in Ecal towers
102  for (int i = 0; i<49; i++)
103  {
104  sprintf(id, "h%d",i+200);
105  sprintf(ntupletag, "Energy Deposit in Ecal Tower%d in GeV",i);
106  G4int histoID = analysisManager->CreateH1(id,ntupletag, 100, 0., 1.0);
107  if (!i)
108  ecalE = histoID;
109  }
110  // Total energy deposit
111  energy = analysisManager->CreateH1("h4000", "Total energy deposited in GeV",
112  100, 0., 100.0);
113 
114  // Time slices
115  for (int i=0; i<numberOfTimeSlices; i++){
116  sprintf(id, "h%d",i+300);
117  sprintf(ntupletag, "Time slice %d nsec energy profile in GeV",i);
118  G4int histoID = analysisManager->CreateH1(id,ntupletag,100, 0., 100.0);
119  if (!i)
120  timeHist = histoID;
121  }
122 
123  // Profile of lateral energy deposit in Hcal
124  for (int i = 0; i<70; i++) {
125  sprintf(id, "h%d",i+500);
126  sprintf(ntupletag, "Lateral energy profile at %d cm in GeV",i);
127  G4int histoID = analysisManager->CreateH1(id, ntupletag, 100, 0., 10.0);
128  if (!i)
129  lateralProfile = histoID;
130  }
131 
132  // Time profile
133  timeProfile = analysisManager->CreateH1("h901", "Time Profile in Sensitive Detector",
134  200, 0., 200.);
135  analysisManager->CreateH1("h902", "Time Profile in Sensitive+Passive",
136  200, 0., 200.);
137 }
138 
139 
141  Finish();
142  delete G4AnalysisManager::Instance();
143 }
144 
145 
147 }
148 
150 {;}
151 
153  if (instance == 0) instance = new CCalAnalysis();
154  return instance;
155 }
156 
157 
158 // This function fill the 1d histogram of the energies in HCal layers
160 {
161  G4AnalysisManager* man = G4AnalysisManager::Instance();
162  G4double totalFilledEnergyHcal = 0.0;
163  for (int i=0; i<28; i++) {
164  G4double x = v[i];
165  man->FillH1(hcalE+i,x);
166  if (fVerbosity)
167  {
168  G4cout << "Fill Hcal histo " << i << " with " << x << G4endl;
169  totalFilledEnergyHcal += x;
170  }
171  }
172 
173  if (fVerbosity)
174  G4cout <<
175  "CCalAnalysis::InsertEnergyHcal: Total filled Energy Hcal histo "
176  << totalFilledEnergyHcal << G4endl;
177 }
178 
179 
180 // This function fill the 1d histogram of the energies in ECal layers
182 {
183  G4AnalysisManager* man = G4AnalysisManager::Instance();
184 
185  G4double totalFilledEnergyEcal = 0.0;
186 
187  for (G4int i=0; i<49; i++) {
188  G4double x = v[i];
189  man->FillH1(ecalE+i,x);
190  if (fVerbosity)
191  {
192  G4cout << "Fill Ecal histo " << i << " with " << x << G4endl;
193  totalFilledEnergyEcal += x;
194  }
195  }
196  if (fVerbosity)
197  G4cout <<
198  "CCalAnalysis::InsertEnergyEcal: Total filled Energy Ecal histo "
199  << totalFilledEnergyEcal << G4endl;
200 }
201 
202 
203 // This function fill the 1d histogram of the lateral profile
205 {
206  G4AnalysisManager* man = G4AnalysisManager::Instance();
207  G4double totalFilledProfileHcal = 0.0;
208 
209  for (G4int i=0; i<70; i++) {
210  G4double x = v[i];
211  man->FillH1(lateralProfile+1,x);
212  if (fVerbosity)
213  {
214  G4cout << "Fill Profile Hcal histo " << i << " with " << x << G4endl;
215  totalFilledProfileHcal += x;
216  }
217  }
218  if (fVerbosity)
219  G4cout << "CCalAnalysis::InsertLateralProfile: Total filled Profile Hcal"
220  << " histo " << totalFilledProfileHcal << G4endl;
221 }
222 
223 
224 // This function fill the 1d histogram of the energy
226 {
227 
228  G4double x = v;
229  G4AnalysisManager::Instance()->FillH1(energy,x);
230  if (fVerbosity)
231  G4cout << "CCalAnalysis::InsertEnergy: Fill Total energy Hcal histo with "
232  << x << G4endl;
233 }
234 
235 
236 // This function fill the 1d histograms of time profiles at stepping action
238 {
239  G4AnalysisManager* man = G4AnalysisManager::Instance();
240 
241  G4double totalFilledTimeProfile = 0.0;
242  for (G4int j=0; j<numberOfTimeSlices; j++)
243  {
244  G4double x = v[j];
245  man->FillH1(timeHist+j,x);
246  if (fVerbosity)
247  {
248  G4cout << "Fill Time slice histo " << j << " with " << x << G4endl;
249  totalFilledTimeProfile += x;
250  }
251 
252  G4double t = j + 0.5;
253  man->FillH1(timeProfile+1,t,x);
254  if (fVerbosity)
255  G4cout << "Fill Time profile histo 1 with " << t << " " << x << G4endl;
256  }
257  if (fVerbosity)
258  G4cout << "CCalAnalysis::InsertTime: Total filled Time profile histo "
259  << totalFilledTimeProfile << G4endl;
260 }
261 
262 
263 // This function fill the 1d histograms of time profiles in SD
264 void CCalAnalysis::InsertTimeProfile(int hit, double time, double edep)
265 {
266  G4AnalysisManager* man = G4AnalysisManager::Instance();
267  man->FillH1(timeProfile,time,edep);
268 
269  if (fVerbosity)
270  G4cout << "CCalAnalysis:: Fill Time Profile with Hit " << hit
271  << " Edeposit " << edep << " Gev at " << time << " ns" << G4endl;
272 }
273 
274 
275 
276 void CCalAnalysis::setNtuple(float* HCalE, float* ECalE, float elab,
277  float x, float y, float z, float edep,
278  float edec, float edhc)
279 {
280  G4AnalysisManager* man = G4AnalysisManager::Instance();
281  G4int counter=0;
282  for (int i=0; i<28; i++)
283  {
284  man->FillNtupleFColumn(counter,HCalE[i]);
285  counter++;
286  }
287  for (int i=0; i<49; i++)
288  {
289  man->FillNtupleFColumn(counter,ECalE[i]);
290  counter++;
291  }
292  man->FillNtupleFColumn(counter,elab);
293  man->FillNtupleFColumn(counter+1,x);
294  man->FillNtupleFColumn(counter+2,y);
295  man->FillNtupleFColumn(counter+3,z);
296  man->FillNtupleFColumn(counter+4,edep);
297  man->FillNtupleFColumn(counter+5,edec);
298  man->FillNtupleFColumn(counter+6,edhc);
299 
300  man->AddNtupleRow();
301 
302  if (fVerbosity)
303  G4cout << "CCalAnalysis:: Fill Ntuple " << G4endl;
304 }
305 
306 
307 /*
308  This member reset the histograms and it is called at the begin
309  of each run; here we put the inizialization so that the histograms have
310  always the right dimensions depending from the detector geometry
311 */
313 {
314  G4AnalysisManager* man = G4AnalysisManager::Instance();
315  for (int i=0; i<numberOfTimeSlices; i++){
316  man->GetH1(timeHist+1)->reset();
317  }
318 }
319 
320 //============================================================================
321 
322 // This member is called at the end of each run
324 {
325  // Save histograms
326  G4AnalysisManager* man = G4AnalysisManager::Instance();
327  man->Write();
328  man->CloseFile();
329 }
330 
331 
332 // This member is called at the end of every event
334  // The plotter is updated only if there is some
335  // hits in the event
336  if (!flag) return;
337 }
338 
339 
virtual ~CCalAnalysis()
G4bool SetFirstHistoId(G4int firstId)
void setNtuple(float *hcalE, float *ecalE, float elab, float x, float y, float z, float edep, float edec, float edhc)
G4int CreateH1(const G4String &name, const G4String &title, G4int nbins, G4double xmin, G4double xmax, const G4String &unitName="none", const G4String &fcnName="none", const G4String &binSchemeName="linear")
G4int timeProfile
Definition: CCalAnalysis.hh:91
G4double z
Definition: TRTMaterials.hh:39
void SetVerboseLevel(G4int verboseLevel)
void InsertEnergyHcal(float *)
G4int CreateNtuple(const G4String &name, const G4String &title)
G4int fVerbosity
Definition: CCalAnalysis.hh:76
int G4int
Definition: G4Types.hh:78
G4bool OpenFile(const G4String &fileName="")
void InsertEnergyEcal(float *)
G4bool FillNtupleFColumn(G4int id, G4float value)
G4GLOB_DLL std::ostream G4cout
G4int numberOfTimeSlices
Definition: CCalAnalysis.hh:78
void InsertEnergy(float v)
void InsertLateralProfile(float *)
void EndOfEvent(G4int flag)
G4bool SetFirstNtupleId(G4int firstId)
Defines the format which is used for the output file default is a ROOT file. Comment the g4root and u...
Definition: CCalAnalysis.hh:44
G4double energy(const ThreeVector &p, const G4double m)
G4bool FillH1(G4int id, G4double value, G4double weight=1.0)
void BeginOfRun(G4int n)
static CCalAnalysis * getInstance()
G4int CreateNtupleFColumn(const G4String &name)
void EndOfRun(G4int n)
static CCalAnalysis * instance
Definition: CCalAnalysis.hh:74
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4int lateralProfile
Definition: CCalAnalysis.hh:89
tools::histo::h1d * GetH1(G4int id, G4bool warn=true, G4bool onlyIfActive=true) const
void InsertTimeProfile(int, double, double)
void InsertTime(float *)