source: MEMPHYS/HEAD/applications/MEMPHYS_analysis.cxx @ 171

Last change on this file since 171 was 171, checked in by barrand, 18 years ago
File size: 16.2 KB
Line 
1// AIDA analysis program over MEMPHYS.root file.
2
3// AIDA :
4#include <AIDA/AIDA.h>
5
6// std :
7#include <iostream>
8
9// Lib :
10#include <Lib/smanip.h>
11
12//////////////////////////////////////////////////////////////////////////////
13bool get_XYZ(
14 AIDA::ITuple& aParent
15,int aColumn
16,double& aX
17,double& aY
18,double& aZ
19) 
20//////////////////////////////////////////////////////////////////////////////
21//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
22{
23  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(aColumn);
24  if(!entry) return false;
25
26  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
27  if(!tuple) return false;
28
29  tuple->start();
30  if(!tuple->next()) return false;
31
32  aX = tuple->getDouble(0);
33  aY = tuple->getDouble(1);
34  aZ = tuple->getDouble(2);
35
36  return true;
37}
38//////////////////////////////////////////////////////////////////////////////
39//////////////////////////////////////////////////////////////////////////////
40//////////////////////////////////////////////////////////////////////////////
41bool dump_tracks(
42 AIDA::ITuple& aParent
43) 
44//////////////////////////////////////////////////////////////////////////////
45//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
46{
47  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(8);
48  if(!entry) return false;
49
50  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
51  if(!tuple) return false;
52
53  //if(nPart!=tracks->rows()) {
54  //  std::cout << "read: nPart / tracks mismatch "
55  //              << nPart << " " << tracks->rows() << std::endl;
56  //  return false;
57  //}
58
59  tuple->start();
60  int irow = 0;
61  while(tuple->next()) {
62
63    int pId = tuple->getInt(0);
64    int parent = tuple->getInt(1);
65    float timeStart = tuple->getFloat(2);
66
67    double dx,dy,dz;
68    if(!get_XYZ(*tuple,3,dx,dy,dz)) return false;
69
70    double mass = tuple->getDouble(4);
71    double pTot = tuple->getDouble(5);
72    double ETot = tuple->getDouble(6);
73
74    double px,py,pz;
75    if(!get_XYZ(*tuple,7,px,py,pz)) return false;
76
77    double start_x,start_y,start_z;
78    if(!get_XYZ(*tuple,8,start_x,start_y,start_z)) return false;
79
80    double stop_x,stop_y,stop_z;
81    if(!get_XYZ(*tuple,9,stop_x,stop_y,stop_z)) return false;
82
83    int startVol = tuple->getInt(10);
84    int stopVol = tuple->getInt(11);
85
86    std::cout << "----> Tk{"<<irow<<"}: " 
87              << " pId " << pId
88              << " parent " << parent
89              << " creation time " << timeStart
90              << " Volumes " << startVol << " " << stopVol << "\n"
91              << " Start Pos (" << start_x
92              << "," << start_y << "," << start_z << ")\n"
93              << " Stop Pos (" << stop_x
94              << "," << stop_y << "," << stop_z << ")\n"
95              << " dx,dy,dz " << dx << " " << dy << " " << dz << "\n"
96              << " m " << mass
97              << " ETot " << ETot
98              << " pTot " << pTot
99              << " px,py,pz " << px << " " << py << " " << pz << "\n"
100              << std::endl;
101
102    irow++;
103  }
104  return true;
105}
106//////////////////////////////////////////////////////////////////////////////
107bool process_hits_times(
108 AIDA::ITuple& aParent
109,AIDA::IHistogram1D& aHisto
110) 
111//////////////////////////////////////////////////////////////////////////////
112//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
113{
114  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(2);
115  if(!entry) return false;
116
117  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
118  if(!tuple) return false;
119
120  tuple->start();
121  while(tuple->next()) {
122
123    float time = tuple->getFloat(0);
124
125    aHisto.fill(time);
126  }
127  return true;
128}
129//////////////////////////////////////////////////////////////////////////////
130bool process_hits(
131 AIDA::ITuple& aParent
132,AIDA::IHistogram1D& aHisto
133) 
134//////////////////////////////////////////////////////////////////////////////
135//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
136{
137  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(10);
138  if(!entry) return false;
139
140  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
141  if(!tuple) return false;
142
143  tuple->start();
144  while(tuple->next()) {
145
146    //int tubeId = tuple->getInt(0);
147    //int totalPE = tuple->getInt(1);
148
149    if(!process_hits_times(*tuple,aHisto)) return false;
150  }
151  return true;
152}
153//////////////////////////////////////////////////////////////////////////////
154bool process_digits(
155 AIDA::ITuple& aParent
156,AIDA::IHistogram2D& aHisto
157) 
158//////////////////////////////////////////////////////////////////////////////
159//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
160{
161  AIDA::ITupleEntry* entry = (AIDA::ITupleEntry*)aParent.getObject(13);
162  if(!entry) return false;
163
164  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(entry);
165  if(!tuple) return false;
166
167  tuple->start();
168  while(tuple->next()) {
169
170    //int tubeId = tuple->getInt(0);
171    double pe = tuple->getDouble(1);
172    double time = tuple->getDouble(2);
173    //printf("debug : ++++ : %g %g\n",time,pe);
174
175    aHisto.fill(time,pe);
176  }
177  return true;
178}
179//////////////////////////////////////////////////////////////////////////////
180void set_region_style(
181 AIDA::IPlotterRegion& aRegion
182) 
183//////////////////////////////////////////////////////////////////////////////
184// Taken from G4SPL_analysis.cxx.
185//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
186{
187  //FIXME : have to do the below with AIDA styles.
188
189  // ROOT is in NDC, then we take a plotter with width = height = 1
190
191  // PAW :
192  //float XSIZ = 20.0F;
193  //float XLAB = 1.4F;
194  //float YSIZ = 20.0F;
195  //float YLAB = 1.4F;
196
197  float XLAB = 1.4F/20.0F; //0.07 //x distance of y title to data frame.
198  float YLAB = 0.8F/20.0F; //0.04 //y distance of x title to data frame.
199  YLAB = 0.05F; //FIXME : cooking.
200
201  // To have a good matching with ROOT for text size :
202  //double majic = 0.014/0.027;
203  double majic = 0.6;
204  std::string s = Lib::smanip::tostring(majic);
205  //FIXME : have methods for the below :
206  aRegion.style().setParameter("textScale",s);
207  aRegion.style().xAxisStyle().setParameter("magStyle.scale",s);
208  aRegion.style().yAxisStyle().setParameter("magStyle.scale",s);
209  aRegion.style().xAxisStyle().tickLabelStyle().setParameter("scale",s);
210  aRegion.style().xAxisStyle().labelStyle().setParameter("scale",s);
211  aRegion.style().yAxisStyle().tickLabelStyle().setParameter("scale",s);
212  aRegion.style().yAxisStyle().labelStyle().setParameter("scale",s);
213
214  // ROOT def margins 0.1. SoPlotter def 0.1.
215  aRegion.layout().setParameter("rightMargin",0.1);
216  aRegion.layout().setParameter("topMargin",0.1);
217  //gStyle->SetPadBottomMargin(0.15);
218  aRegion.layout().setParameter("leftMargin",0.15);
219  //gStyle->SetPadLeftMargin(0.15);
220  aRegion.layout().setParameter("bottomMargin",0.15);
221
222  aRegion.style().setParameter("superposeBins","TRUE");
223  aRegion.setParameter("plotter.wallStyle.visible","FALSE");
224  aRegion.setParameter("plotter.gridStyle.visible","FALSE");
225
226  //std::string font = "Hershey"; //PAW font.
227  //std::string font = "TTF/arialbd"; //Default ROOT font 62.
228  //std::string font = "TTF/couri";
229  std::string font = "TTF/times"; //= ROOT 132 ?
230  //std::string smoothing = "FALSE";
231  std::string smoothing = "TRUE";
232
233  // X axis :
234  aRegion.style().xAxisStyle().tickLabelStyle().setFont(font);
235  aRegion.style().xAxisStyle().tickLabelStyle().setParameter("smoothing",smoothing);
236  aRegion.style().xAxisStyle().labelStyle().setFont(font);
237  aRegion.style().xAxisStyle().labelStyle().setParameter("smoothing",smoothing);
238  aRegion.setParameter("plotter.xAxis.magStyle.fontName",font);
239  aRegion.setParameter("plotter.xAxis.magStyle.smoothing",smoothing);
240  aRegion.style().xAxisStyle().lineStyle().setThickness(2);
241  aRegion.setParameter("plotter.xAxis.ticksStyle.width","2");
242  aRegion.style().setParameter("topAxisVisible","TRUE");
243  // Set hplot tick modeling :
244  aRegion.style().xAxisStyle().setParameter("modeling","hplot");
245  aRegion.style().xAxisStyle().setParameter("divisions","505");
246
247  // Y axis :
248  //aRegion.style().setParameter("yAxisLogScale","TRUE");
249  aRegion.style().yAxisStyle().tickLabelStyle().setFont(font);
250  aRegion.style().yAxisStyle().tickLabelStyle().setParameter("smoothing",smoothing);
251  aRegion.style().yAxisStyle().labelStyle().setFont(font);
252  aRegion.style().yAxisStyle().labelStyle().setParameter("smoothing",smoothing);
253  aRegion.setParameter("plotter.yAxis.magStyle.fontName",font);
254  aRegion.setParameter("plotter.yAxis.magStyle.smoothing",smoothing);
255  aRegion.style().yAxisStyle().lineStyle().setThickness(2);
256  aRegion.setParameter("plotter.yAxis.ticksStyle.width","2");
257
258  //gStyle->SetTitleSize(0.06,"XYZ"); //ROOT def 0.04. SoPlotter : 0.014
259  aRegion.setParameter("plotter.xAxis.titleHeight","0.06");
260  aRegion.setParameter("plotter.yAxis.titleHeight","0.06");
261  //gStyle->SetLabelOffset(0.01,"Y"); //ROOT def 0.005. SoPlotter def 0.02
262  aRegion.setParameter("plotter.yAxis.labelToAxis","0.01");
263  //gStyle->SetTitleOffset(1.1,"Y");
264  aRegion.setParameter("plotter.yAxis.titleToAxis",
265   Lib::smanip::tostring(1.1*XLAB));
266  aRegion.setParameter("plotter.xAxis.titleToAxis",
267   Lib::smanip::tostring(1.1*YLAB));
268  //gStyle->SetLabelSize(0.05,"XYZ"); //ROOT def 0.04. SoPlotter def 0.014.
269  aRegion.setParameter("plotter.yAxis.labelHeight","0.05");
270  aRegion.setParameter("plotter.xAxis.labelHeight","0.05");
271
272  aRegion.style().setParameter("rightAxisVisible","TRUE");
273  // Set hplot tick modeling :
274  aRegion.style().yAxisStyle().setParameter("modeling","hplot");
275  aRegion.style().yAxisStyle().setParameter("divisions","505");
276
277  // title :
278  aRegion.style().setParameter("titleHeight","0.04");
279  aRegion.style().setParameter("titleToAxis","0.02");
280  aRegion.style().titleStyle().textStyle().setFont(font);
281  aRegion.style().titleStyle().textStyle().setParameter("smoothing",smoothing);
282  aRegion.style().titleStyle().textStyle().setColor("0 0 1");
283
284  // legend box :
285  aRegion.setParameter("legendRegionVisible","TRUE");
286  aRegion.setParameter("legendRegion.viewportRegion.backgroundColor","1 1 0");
287  aRegion.setParameter("legendRegionSize","0.15 0.16");
288  aRegion.setParameter("legendRegionOrigin","0.1 0.1");
289  aRegion.setParameter("legendRegion.horizontalMargin","2");
290  aRegion.setParameter("legendRegion.verticalMargin","2");
291
292  // Frame :
293  aRegion.setParameter("plotter.innerFrameStyle.lineWidth","2");
294}
295//////////////////////////////////////////////////////////////////////////////
296bool plot(
297 AIDA::IAnalysisFactory& aAIDA
298,AIDA::IHistogram1D& aHisto1D
299,AIDA::IHistogram2D& aHisto2D
300) 
301//////////////////////////////////////////////////////////////////////////////
302//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
303{
304  std::string header = "MEMPHYS_analysis : ";
305
306  AIDA::IPlotterFactory* plotterFactory = 
307    aAIDA.createPlotterFactory(0,0);
308  if(!plotterFactory) {
309    std::cout << header << "can't get a PlotterFactory." << std::endl;
310    return false;
311  }
312   
313  AIDA::IPlotter* plotter = plotterFactory->create();
314  if(!plotter) {
315    std::cout << header << "can't get a plotter." << std::endl;
316    return false;
317  }
318   
319  plotter->clearRegions();
320  plotter->createRegions(1,2,0);
321
322 {plotter->setCurrentRegionNumber(0);
323  AIDA::IPlotterRegion& region = plotter->currentRegion();
324  set_region_style(region);
325  region.style().xAxisStyle().setLabel("time");
326  region.style().yAxisStyle().setLabel("Entries");
327  region.plot(aHisto1D);}
328
329 {plotter->setCurrentRegionNumber(1);
330  AIDA::IPlotterRegion& region = plotter->currentRegion();
331  set_region_style(region);
332  region.style().xAxisStyle().setLabel("time");
333  region.style().yAxisStyle().setLabel("PE");
334  region.plot(aHisto2D);}
335
336  plotter->show();
337
338  /* FIXME : try to have multiple plotters in the tab stack :
339 {AIDA::IPlotter* plotter = plotterFactory->create("Viewer_2");
340  if(!plotter) {
341    std::cout << header << "can't get a plotter." << std::endl;
342    return false;
343  }
344  plotter->clearRegions();
345  plotter->createRegions(1,1,0);}
346  */
347
348  plotter->interact();
349
350  delete plotter;
351  delete plotterFactory;
352
353  return true;
354}
355//////////////////////////////////////////////////////////////////////////////
356//////////////////////////////////////////////////////////////////////////////
357//////////////////////////////////////////////////////////////////////////////
358int main(
359 int //aArgc
360,char** //aArgv
361) 
362//////////////////////////////////////////////////////////////////////////////
363//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
364{
365  //std::string arg1 = aArgc>=2?aArgv[1]:"";
366
367  std::string program = "MEMPHYS_analysis";
368  std::string header = program + " : ";
369
370  AIDA::IAnalysisFactory* aida = AIDA_createAnalysisFactory();
371  if(!aida) {
372    std::cout << header << "AIDA not found." << std::endl;
373    return 1;
374  }
375
376  AIDA::ITreeFactory* treeFactory = aida->createTreeFactory();
377  if(!treeFactory) {
378    std::cout << header << "can't get a TreeFactory." << std::endl;
379    return 1;
380  }
381
382  ////////////////////////////////////////////////////////
383  // Create histograms in memory tree ////////////////////
384  ////////////////////////////////////////////////////////
385  AIDA::ITree* memory = treeFactory->create();
386  if(!memory) {
387    std::cout << header << "can't get memory tree." << std::endl;
388    return 1;
389  }
390  AIDA::IHistogramFactory* histogramFactory = 
391    aida->createHistogramFactory(*memory);
392  if(!histogramFactory) {
393    std::cout << header << "can't get an histogram factory." << std::endl;
394    return 1;
395  }
396 
397  AIDA::IHistogram1D* hits_times = 
398    histogramFactory->createHistogram1D("hits_times","Hits times",100,0,3000);
399  if(!hits_times) {
400    std::cout << header
401              << "can't create histogram : time." 
402              << std::endl;
403    return 1;
404  }
405
406  AIDA::IHistogram2D* digits_time_pe = 
407    histogramFactory->createHistogram2D("digits_pe_time","Digits PE time",
408                                        100,0,3000,100,0,10);
409  if(!digits_time_pe) {
410    std::cout << header
411              << "can't create histogram : digits_time_pe." 
412              << std::endl;
413    return 1;
414  }
415
416  delete histogramFactory;
417
418  ////////////////////////////////////////////////////////
419  /// Read data //////////////////////////////////////////
420  ////////////////////////////////////////////////////////
421
422  AIDA::ITree* tree = treeFactory->create("MEMPHYS.root","root",true,false);
423  if(!tree) {
424    std::cout << header << "can't open data file." << std::endl;
425    return 1;
426  }
427
428  AIDA::IManagedObject* object = tree->find("Event");
429  if(!object) {
430    std::cout << header
431              << "object Event not found in tree." 
432              << std::endl;
433    return 1;
434  }
435  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(object);
436  if(!tuple) {
437    std::cout << header << "object not an AIDA::ITuple." << std::endl;
438    return 1;
439  }
440
441  int coln = tuple->columns();
442  for(int index=0;index<coln;index++) {
443    std::cout << header
444              << " icol = " << index
445              << ", label = " << tuple->columnName(index) 
446              << ", type = " << tuple->columnType(index) 
447              << std::endl;
448  }
449
450  std::cout << header << "rows = " << tuple->rows() << std::endl;
451
452  int nentries = tuple->rows();
453  //nentries = 100000;
454  //nentries = 40;
455  std::cout << header
456            << "traitements de " << nentries << " entrees" 
457            << std::endl;
458
459  int irow = 0;
460  tuple->start();
461  while(tuple->next() && (irow<nentries)) {
462
463    /*
464    int eventId = tuple->getInt(0);
465    //int inputEvtId = tuple->getInt(1);
466    //int interMode = tuple->getInt(2);
467    //int vtxVol = tuple->getInt(3);
468
469    int nPart = tuple->getInt(5);
470    //int leptonIndex = tuple->getInt(6);
471    //int protonIndex = tuple->getInt(7);
472
473    int nHits = tuple->getInt(9);
474    int nDigits = tuple->getInt(11);
475    double sumPE = tuple->getDouble(12);
476
477    std::cout << ">>>>>>>>>>>>> Event{" << irow << "}: "
478              << " evt Id " << eventId
479    //        << " evt Input Id " << inputEvtId
480    //        << "\n interaction mode " << interMode
481    //        << " start in volume " << vtxVol << "\n"
482              <<" #tracks: " << nPart
483              <<" #hits: " << nHits
484              <<" #digits: " << nDigits
485              <<" sumPE " << sumPE
486              << std::endl;
487    */
488
489    //if(!dump_tracks(*tuple)) return 1;
490
491    if(!process_hits(*tuple,*hits_times)) return 1;
492    if(!process_digits(*tuple,*digits_time_pe)) return 1;
493
494    irow++;
495  }
496
497  delete treeFactory;
498  delete tree;
499
500  plot(*aida,*hits_times,*digits_time_pe);
501
502  delete aida;
503
504  return 0;
505}
Note: See TracBrowser for help on using the repository browser.