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

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

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