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

Last change on this file was 1337, checked in by garnier, 14 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.