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

Last change on this file since 169 was 169, checked in by barrand, 18 years ago

G.Barrand : ajout MEMPHYS_analysis
File size: 15.9 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  plotter->interact();
338
339  delete plotter;
340  delete plotterFactory;
341
342  return true;
343}
344//////////////////////////////////////////////////////////////////////////////
345//////////////////////////////////////////////////////////////////////////////
346//////////////////////////////////////////////////////////////////////////////
347int main(
348 int //aArgc
349,char** //aArgv
350) 
351//////////////////////////////////////////////////////////////////////////////
352//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
353{
354  //std::string arg1 = aArgc>=2?aArgv[1]:"";
355
356  std::string program = "MEMPHYS_analysis";
357  std::string header = program + " : ";
358
359  AIDA::IAnalysisFactory* aida = AIDA_createAnalysisFactory();
360  if(!aida) {
361    std::cout << header << "AIDA not found." << std::endl;
362    return 1;
363  }
364
365  AIDA::ITreeFactory* treeFactory = aida->createTreeFactory();
366  if(!treeFactory) {
367    std::cout << header << "can't get a TreeFactory." << std::endl;
368    return 1;
369  }
370
371  ////////////////////////////////////////////////////////
372  // Create histograms in memory tree ////////////////////
373  ////////////////////////////////////////////////////////
374  AIDA::ITree* memory = treeFactory->create();
375  if(!memory) {
376    std::cout << header << "can't get memory tree." << std::endl;
377    return 1;
378  }
379  AIDA::IHistogramFactory* histogramFactory = 
380    aida->createHistogramFactory(*memory);
381  if(!histogramFactory) {
382    std::cout << header << "can't get an histogram factory." << std::endl;
383    return 1;
384  }
385 
386  AIDA::IHistogram1D* hits_times = 
387    histogramFactory->createHistogram1D("hits_times","Hits times",100,0,3000);
388  if(!hits_times) {
389    std::cout << header
390              << "can't create histogram : time." 
391              << std::endl;
392    return 1;
393  }
394
395  AIDA::IHistogram2D* digits_time_pe = 
396    histogramFactory->createHistogram2D("digits_pe_time","Digits PE time",
397                                        100,0,3000,100,0,10);
398  if(!digits_time_pe) {
399    std::cout << header
400              << "can't create histogram : digits_time_pe." 
401              << std::endl;
402    return 1;
403  }
404
405  delete histogramFactory;
406
407  ////////////////////////////////////////////////////////
408  /// Read data //////////////////////////////////////////
409  ////////////////////////////////////////////////////////
410
411  AIDA::ITree* tree = treeFactory->create("MEMPHYS.root","root",true,false);
412  if(!tree) {
413    std::cout << header << "can't open data file." << std::endl;
414    return 1;
415  }
416
417  AIDA::IManagedObject* object = tree->find("Event");
418  if(!object) {
419    std::cout << header
420              << "object Event not found in tree." 
421              << std::endl;
422    return 1;
423  }
424  AIDA::ITuple* tuple = dynamic_cast<AIDA::ITuple*>(object);
425  if(!tuple) {
426    std::cout << header << "object not an AIDA::ITuple." << std::endl;
427    return 1;
428  }
429
430  int coln = tuple->columns();
431  for(int index=0;index<coln;index++) {
432    std::cout << header
433              << " icol = " << index
434              << ", label = " << tuple->columnName(index) 
435              << ", type = " << tuple->columnType(index) 
436              << std::endl;
437  }
438
439  std::cout << header << "rows = " << tuple->rows() << std::endl;
440
441  int nentries = tuple->rows();
442  //nentries = 100000;
443  nentries = 40;
444  std::cout << header
445            << "traitements de " << nentries << " entrees" 
446            << std::endl;
447
448  int irow = 0;
449  tuple->start();
450  while(tuple->next() && (irow<nentries)) {
451
452    /*
453    int eventId = tuple->getInt(0);
454    //int inputEvtId = tuple->getInt(1);
455    //int interMode = tuple->getInt(2);
456    //int vtxVol = tuple->getInt(3);
457
458    int nPart = tuple->getInt(5);
459    //int leptonIndex = tuple->getInt(6);
460    //int protonIndex = tuple->getInt(7);
461
462    int nHits = tuple->getInt(9);
463    int nDigits = tuple->getInt(11);
464    double sumPE = tuple->getDouble(12);
465
466    std::cout << ">>>>>>>>>>>>> Event{" << irow << "}: "
467              << " evt Id " << eventId
468    //        << " evt Input Id " << inputEvtId
469    //        << "\n interaction mode " << interMode
470    //        << " start in volume " << vtxVol << "\n"
471              <<" #tracks: " << nPart
472              <<" #hits: " << nHits
473              <<" #digits: " << nDigits
474              <<" sumPE " << sumPE
475              << std::endl;
476    */
477
478    //if(!dump_tracks(*tuple)) return 1;
479
480    if(!process_hits(*tuple,*hits_times)) return 1;
481    if(!process_digits(*tuple,*digits_time_pe)) return 1;
482
483    irow++;
484  }
485
486  delete treeFactory;
487  delete tree;
488
489  plot(*aida,*hits_times,*digits_time_pe);
490
491  delete aida;
492
493  return 0;
494}
Note: See TracBrowser for help on using the repository browser.