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

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