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