source: trunk/examples/advanced/underground_physics/src/DMXAnalysisManager.cc@ 1343

Last change on this file since 1343 was 1337, checked in by garnier, 15 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 15.9 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// $Id: DMXAnalysisManager.cc
28// GEANT4 tag $Name:
29//
30// Author: Alex Howard (alexander.howard@cern.ch)
31//
32// History:
33// -----------
34// 16 Jan 2002 Alex Howard Created
35// 17 June 2002 Alex Howard Successfully Modified to AIDA 2.2
36// 17 November 2002 Alex Howard Migrated to AIDA 3.0 and added fitting
37// 22 October 2009 Luciano Pandola Changed variables names for ntuple2
38// (avoid problem with columns) and some
39// clean-up (e.g. remove ntuple4)
40//
41// -------------------------------------------------------------------
42#ifdef G4ANALYSIS_USE
43
44#include "DMXAnalysisManager.hh"
45
46DMXAnalysisManager* DMXAnalysisManager::instance = 0;
47
48//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50DMXAnalysisManager::DMXAnalysisManager() :
51 af(0), tf(0),tree(0), hf(0), tpf(0),
52 ntuple1(0), ntuple2(0), ntuple3(0),
53 hEsourcep(0), hEdepp(0), hEdepRecoil(0), hNumPhLow(0), hNumPhHigh(0),
54 hAvPhArrival(0), hPhArrival(0), hPMTHits(0), h1stPMTHit(0),hGammaEdep(0),
55 hNeutronEdep(0), hElectronEdep(0), hPositronEdep(0), hOtherEdep(0)
56 //,funFact(0),fitFact(0),exponFun(0),gaussFun(0),e_fitter(0),fitResult(0)
57{;}
58
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60
61DMXAnalysisManager::~DMXAnalysisManager()
62{;}
63
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65
66DMXAnalysisManager* DMXAnalysisManager::getInstance()
67{
68 if (instance == 0) instance = new DMXAnalysisManager();
69 return instance;
70}
71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72
73
74void DMXAnalysisManager::book(G4String histogramfile, G4bool)
75
76{
77 // histoManager->selectStore("DMX.his");
78 G4cout << " Histogramfile: " << histogramfile << G4endl;
79
80
81 //build up the factories
82 af = AIDA_createAnalysisFactory();
83
84
85 //parameters for the TreeFactory
86 G4bool fileExists = true;
87 G4bool readOnly = false;
88
89 tf = af->createTreeFactory();
90
91 tree = tf->create(histogramfile, "hbook", readOnly, fileExists);
92
93 G4cout << "Tree store : " << tree->storeName() << G4endl;
94
95 G4cout << " Booked Hbook File " << G4endl;
96
97 //HistoFactory and TupleFactory depend on theTree
98 hf = af->createHistogramFactory( *tree );
99 tpf = af->createTupleFactory(*tree );
100
101 // ---- primary ntuple ------
102
103 ntuple1 = tpf->create( "1", "Particle Source Energy",
104 "double energy" );
105
106 // ---- secondary ntuple ------
107
108 ntuple2 = tpf->create( "2", "Scintillation Hits Info",
109 "float Event,e_prim,tot_e,s_hits,xe_time,num_ph,avphtime,firstpart,firstparte,gamma,neutron,posi,elec,other,seed1,seed2" );
110
111 // ---- tertiary ntuple ------
112
113 ntuple3 = tpf->create( "3", "PMT Hits Info",
114 "float event, hits, xpos, ypos, zpos" );
115
116 // Creating an 1-dimensional histogram in the root directory of the tree
117
118 hEsourcep = hf->createHistogram1D("10","Source Energy /keV", 1000,0.,10000.);
119
120 hEdepp = hf->createHistogram1D("20","Energy Deposit /keV", 1000,0.,1000.);
121
122 hEdepRecoil = hf->createHistogram1D("30","Nuclear Recoil Edep /keV", 100,0.,100.);
123
124 hNumPhLow = hf->createHistogram1D("40","Number of Photons - LowE", 200,0.,200.);
125
126 hNumPhHigh = hf->createHistogram1D("50","Number of Photons - HighE", 100,0.,10000.);
127
128 hAvPhArrival = hf->createHistogram1D("60","Average Photon Arrival/ns", 200,0.,200.);
129
130 hPhArrival = hf->createHistogram1D("61","1st event Photon Arrival", 200,0.,200.);
131
132 //2D histograms:
133 hPMTHits = hf->createHistogram2D("70","PMT Hit Pattern",
134 300 ,-30.,30.,300,-30.,30.);
135
136 h1stPMTHit = hf->createHistogram2D("71","1st event PMT Hit Pattern",
137 300 ,-30.,30.,300,-30.,30.);
138
139 // extra 1D Histos:
140 hGammaEdep = hf->createHistogram1D("91","Gamma Energy Deposit/keV", 1000,0.,1000.);
141
142 hNeutronEdep = hf->createHistogram1D("92","Neutron Ener Deposit/keV", 1000,0.,1000.);
143
144 hElectronEdep = hf->createHistogram1D("93","Electron Ener Deposit/keV",1000,0.,1000.);
145
146 hPositronEdep = hf->createHistogram1D("94","Positron Ener Deposit/keV",1000,0.,1000.);
147
148 hOtherEdep = hf->createHistogram1D("95","Other Ener Deposit/keV", 1000,0.,1000.);
149
150 // // Creating the plotter factory
151// pf = af->createPlotterFactory();
152
153// if(pf && plotevent) {
154// plotter = pf->create();
155// if(plotter) {
156// plotter->show();
157// plotter->setParameter("pageTitle","DMX Output Summary");
158// }
159// delete pf;
160// }
161
162 // create fitting factories ..
163
164 // G4cout << " creating fit factories " << G4endl;
165
166 // funFact = af->createFunctionFactory(*tree);
167// fitFact = af->createFitFactory();
168
169 // G4cout << " creating Exponential Function " << G4endl;
170// exponFun = funFact->createFunctionByName("Exponential","E");
171// G4cout << " created " << G4endl;
172
173// G4cout << " creating Gaussian Function " << G4endl;
174// gaussFun = funFact->createFunctionByName("Gaussian","G");
175// G4cout << " created " << G4endl;
176
177 // assert(exponFun);
178
179 // // e_fitter = fitFact->createFitter("Chi2","","printLevel=-1 errorUP=1.0");
180// G4cout << " creating fitter " << G4endl;
181// e_fitter = fitFact->createFitter();
182
183// // assert(e_fitter);
184
185// G4cout << " Created e_fitter " << G4endl;
186 return;
187
188}
189
190//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
191
192void DMXAnalysisManager::Init()
193{;}
194
195//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196
197void DMXAnalysisManager::Finish()
198{
199
200 // Committing the transaction with the tree
201 G4cout << "Committing..." << G4endl;
202 // write all histograms to file
203 tree->commit();
204
205 G4cout << "Closing the tree..." << G4endl;
206
207 // close (will again commit)
208 tree->close();
209
210
211 delete tpf;
212 delete hf;
213 delete tree;
214 delete tf;
215 delete af;
216
217 // extra delete as objects are created in book() method rather than during
218 // initialisation of class
219 //delete funFact;
220 //delete fitFact;
221
222}
223
224//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
225void DMXAnalysisManager::PulseTimeFit()
226{
227
228 // created fitter factory then fitter
229
230// if(!plotter) plotter = pf->create();
231// if(plotter) {
232// plotter->show();
233// plotter->setParameter("pageTitle","DMX Fitting");
234// }
235
236// G4cout << " Started Pulse Time Fit " << G4endl;
237// int i;
238// for (i = 0; i < 2 ; ++i)
239// G4cout << "\n Fit Parameters: " << i << ":\t"
240// << exponFun->parameterNames()[i] << G4endl;
241// // first fit to time histogram:
242// exponFun->setParameter("amp", 3.);
243// exponFun->setParameter("slope", -0.03);
244// G4cout << " parameters set " << G4endl;
245
246// fitResult = e_fitter->fit(*hPhArrival,*exponFun);
247
248// G4cout << " Fit Completed " << G4endl;
249// G4cout << " chi2/ndf is: " << fitResult->quality() << " / "
250// << fitResult->ndf() << G4endl;
251
252// // retrieve fitted parameters and errors and print them
253
254// for (i = 0; i < 2 ; ++i)
255// G4cout << fitResult->fittedParameterNames()[i] << " = "
256// << fitResult->fittedParameters()[i] << " +/- "
257// << fitResult->errors()[i];
258
259// // plot function with its annotation and then histogram
260
261// plotter->currentRegion().plot(*exponFun,"[0,100] annotation");
262// plotter->currentRegion().plot(*hPhArrival, "overlay");
263// plotter->refresh();
264// plotter->writeToFile("ExponentialFit.ps","ps");
265
266
267// for (i = 0; i < 3 ; ++i)
268// G4cout << "\n Fit Parameters: " << i << ":\t"
269// << gaussFun->parameterNames()[i] << G4endl;
270// // first fit to time histogram:
271// gaussFun->setParameter("mean", 50.);
272// gaussFun->setParameter("sigma", 10.);
273// gaussFun->setParameter("amp", 2.);
274// G4cout << " parameters set " << G4endl;
275
276// fitResult = e_fitter->fit(*hAvPhArrival,*gaussFun);
277
278// G4cout << " Fit Completed " << G4endl;
279// G4cout << " chi2/ndf is: " << fitResult->quality() << " / "
280// << fitResult->ndf() << G4endl;
281
282// // retrieve fitted parameters and errors and print them
283
284// for (i = 0; i < 3 ; ++i)
285// G4cout << fitResult->fittedParameterNames()[i] << " = "
286// << fitResult->fittedParameters()[i] << " +/- "
287// << fitResult->errors()[i];
288
289// // plot function with its annotation and then histogram
290
291// plotter->createRegions(1,1);
292// plotter->currentRegion().plot(*gaussFun,"[0,100] annotation");
293// plotter->currentRegion().plot(*hAvPhArrival, "overlay");
294// plotter->refresh();
295// plotter->writeToFile("TimeConstantFit.ps","ps");
296
297// G4cout << " Fit finished " << G4endl;
298
299}
300
301
302//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
303void DMXAnalysisManager::analyseScintHits(G4int event_id, G4double energy_pri, G4double totEnergy, G4int S_hits, G4double firstLXeHitTime, G4int P_hits, G4double aveTimePmtHits, G4String firstparticleName, G4double firstParticleE, G4bool gamma_ev, G4bool neutron_ev, G4bool positron_ev,G4bool electron_ev,G4bool other_ev, long seed1, long seed2)
304
305{
306
307 G4int firstparticleIndex = 0;
308 if(firstparticleName == "gamma") firstparticleIndex = 1;
309 if(firstparticleName == "neutron") firstparticleIndex = 2;
310 if(firstparticleName == "electron") firstparticleIndex = 3;
311 if(firstparticleName == "positron") firstparticleIndex = 4;
312 if(firstparticleName == "other") {
313 firstparticleIndex = 5;
314 hEdepRecoil->fill(totEnergy);
315 }
316
317 hNumPhLow->fill(P_hits,10.); // fill(x,weight)
318
319 hNumPhHigh->fill(P_hits); // fill(x,y,weight)
320
321
322 hEsourcep->fill( energy_pri/keV );
323
324 hEdepp->fill( totEnergy/keV );
325
326 hAvPhArrival->fill(aveTimePmtHits/ns);
327
328 // Fill the ntuple
329 ntuple2->fill( ntuple2->findColumn( "Event" ), (G4float) event_id );
330 ntuple2->fill( ntuple2->findColumn( "e_prim" ), (G4float) (energy_pri/keV) );
331 ntuple2->fill( ntuple2->findColumn( "tot_e" ), (G4float) totEnergy );
332 ntuple2->fill( ntuple2->findColumn( "s_hits" ), (G4float) S_hits );
333 ntuple2->fill( ntuple2->findColumn( "xe_time" ), (G4float) firstLXeHitTime );
334 ntuple2->fill( ntuple2->findColumn( "num_ph" ), (G4float) P_hits );
335 ntuple2->fill( ntuple2->findColumn( "avphtime"), (G4float) aveTimePmtHits );
336 ntuple2->fill( ntuple2->findColumn( "firstpart" ), (G4float) firstparticleIndex);
337 ntuple2->fill( ntuple2->findColumn( "firstparte"), (G4float) firstParticleE );
338 ntuple2->fill( ntuple2->findColumn( "gamma" ), (G4float) gamma_ev );
339 ntuple2->fill( ntuple2->findColumn( "neutron" ), (G4float) neutron_ev );
340 ntuple2->fill( ntuple2->findColumn( "posi" ), (G4float) positron_ev );
341 ntuple2->fill( ntuple2->findColumn( "elec" ), (G4float) electron_ev );
342 ntuple2->fill( ntuple2->findColumn( "other" ), (G4float) other_ev );
343 ntuple2->fill( ntuple2->findColumn( "seed1" ), (G4float) seed1 );
344 ntuple2->fill( ntuple2->findColumn( "seed2" ), (G4float) seed2 );
345
346 //Values of attributes are prepared; store them to the nTuple:
347 ntuple2->addRow();
348}
349
350//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
351
352void DMXAnalysisManager::analysePMTHits(G4int event, G4int i, G4double x, G4double y, G4double z)
353{
354
355 hPMTHits->fill(x/mm, y/mm); // fill(x,y,weight)
356
357 if (event == 0 ) {
358 h1stPMTHit->fill(x,y);
359 }
360
361 // Fill the ntuple
362 ntuple3->fill( ntuple3->findColumn( "event" ), (G4float) event );
363 ntuple3->fill( ntuple3->findColumn( "hits" ), (G4float) i );
364 ntuple3->fill( ntuple3->findColumn( "xpos" ), (G4float) x );
365 ntuple3->fill( ntuple3->findColumn( "ypos" ), (G4float) y );
366 ntuple3->fill( ntuple3->findColumn( "zpos" ), (G4float) z );
367
368 // NEW: Values of attributes are prepared; store them to the nTuple:
369 ntuple3->addRow();
370
371}
372
373void DMXAnalysisManager::analysePrimaryGenerator(G4double energy)
374{
375
376 // Fill energy ntple:
377 ntuple1->fill( ntuple1->findColumn( "energy" ), energy );
378
379 // NEW: Values of attributes are prepared; store them to the nTuple:
380 ntuple1->addRow();
381
382}
383
384void DMXAnalysisManager::analyseParticleSource(G4double energy, G4String name)
385{
386
387 if(name == "gamma") {
388 hGammaEdep->fill(energy/keV);
389 }
390 if(name == "neutron") {
391 hNeutronEdep->fill(energy/keV); // fill(x,weight)
392 }
393 if(name == "electron") {
394 hElectronEdep->fill(energy/keV); // fill(x,weight)
395 }
396 if(name == "positron") {
397 hPositronEdep->fill(energy/keV); // fill(x,weight)
398 }
399 if(name == "other") {
400 hOtherEdep->fill(energy/keV); // fill(x,weight)
401 }
402
403}
404
405//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406
407void DMXAnalysisManager::HistTime(G4double time)
408{
409 hPhArrival->fill(time/ns); // fill(x,y,weight)
410}
411
412//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
413
414void DMXAnalysisManager::PlotHistosInter(G4int)
415{
416 // // The plotter is updated only if there are some hits in the event
417// if(!flag) return;
418// // Set the plotter ; set the number of regions and attach histograms
419// // to plot for each region.
420// // It is done here, since then EndOfRun set regions
421// // for paper output.
422// if(plotter) {
423// plotter->currentRegion().plot(*hNumPhLow);
424// plotter->refresh();
425// }
426}
427
428void DMXAnalysisManager::PlotHistos(G4bool)
429{
430 /*
431 if(!plotter) plotter = pf->create();
432 plotter->show();
433 plotter->setParameter("pageTitle","DMX Output Summary");
434
435 // G4int n = 1;
436 if(plotter) {
437 // We set one single region for the plotter
438 // We now print the histograms, each one in a separate file
439 plotter->createRegions(1,1);
440 plotter->currentRegion().plot(*hNumPhHigh);
441 plotter->refresh();
442 plotter->writeToFile("NumberPhotons.ps","ps");
443
444 plotter->createRegions(1,1);
445 plotter->currentRegion().plot(*hEdepp);
446 plotter->refresh();
447 plotter->writeToFile("EnergyDeposit.ps","ps");
448
449 plotter->createRegions(1,1);
450 plotter->currentRegion().plot(*hEsourcep);
451 plotter->refresh();
452 plotter->writeToFile("SourceEnergy.ps","ps");
453
454 plotter->createRegions(1,1);
455 plotter->currentRegion().plot(*hNumPhLow);
456 plotter->refresh();
457 plotter->writeToFile("NumberPhotonsLow.ps","ps");
458 }
459 // // Creating two regions
460 // plotter->clearRegions();
461 // // plotter->createRegions(2, 2, 0); // set the current working region to the first one
462 // plotter->show();
463
464 if (interactive) {
465 // Wait for the keyboard return to avoid destroying the plotter window too quickly.
466 G4cout << "Press <ENTER> to exit" << G4endl;
467 G4cin.get();
468 }
469 */
470}
471
472#endif
Note: See TracBrowser for help on using the repository browser.