source: trunk/MagicSoft/Mars/mdata/MDataChain.cc@ 2013

Last change on this file since 2013 was 1853, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 18.6 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 04/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MDataChain
28//
29// With this chain you can concatenate simple mathematical operations on
30// members of mars containers.
31//
32// In the constructor you can give rule, like
33// "HillasSource.fDist / MHillas.fLength"
34// Where MHillas/HillasSource is the name of the parameter container in
35// the parameter list and fDist/fLength is the name of the data members
36// in the containers. The result will be fDist divided by fLength.
37//
38// You can also use brackets:
39// "HillasDource.fDist / (MHillas.fLength + MHillas.fWidth)"
40//
41// The allowed operations are: +, -, *, /
42//
43// Warning: There is no priority rule build in. So better use brackets
44// to get correct results. The rule is parsed/evaluated from the left
45// to the right, which means:
46//
47// "MHillas.fWidth + MHillas.fLength / HillasSource.fDist"
48//
49// is parses as
50//
51// "(MHillas.fWidth + MHillas.fLength) / HillasSource.fDist"
52//
53// You can also use mathmatical operators, eg:
54// "5*log10(MMcEvt.fEnergy*MHillas.fSize)"
55//
56// The allowed operators are:
57// exp(x) e^x
58// log(x) natural logarithm of x
59// pow10(x) 10^x
60// log10(x) logarithm of x to base ten
61// cos(x) cosine of x
62// sin(x) sine of x
63// tan(x) tangent of x
64// cosh(x) hyperbolic cosine of x
65// sinh(x) hyperbolic sine of x
66// tanh(x) hyperbolic tangent of x
67// acos(x) arc cosine (inverse cosine) of x
68// asin(x) arc sine (inverse sine) of x
69// atan(x) arc tangent (inverse tangent) of x
70// sqrt(x) square root of x
71// sqr(x) square of x
72// abs(x) absolute value of x, |x|
73// floor(x) round down to the nearest integer (floor(9.9)=9)
74// r2d(x) transform radians to degrees
75// d2r(x) transform degrees to radians
76// rand(x) returns a uniform deviate on the interval ( 0, x ].
77// (gRandom->Uniform(x) is returned)
78// randp(x) returns gRandom->Poisson(x)
79// rande(x) returns gRandom->Exp(x)
80// randi(x) returns gRandom->Integer(x)
81// randg(x) returns gRandom->Gaus(0, x)
82// randl(x) returns gRandom->Landau(0, x)
83//
84// REMARK:
85// - All the random functions are returning 0 if gRandom==0
86// - You may get better results if you are using a TRandom3
87//
88// FIXME: The possibility to use other objects inheriting from MData
89// is missing.
90// Maybe we can use gInterpreter->Calc("") for this.
91// gROOT->ProcessLineFast("line");
92//
93/////////////////////////////////////////////////////////////////////////////
94
95#include "MDataChain.h"
96
97#include <math.h> // fabs on Alpha
98#include <ctype.h> // isalnum, ...
99#include <stdlib.h> // strtod, ...
100
101#include <TRandom.h>
102
103#include "MLog.h"
104#include "MLogManip.h"
105
106#include "MDataList.h"
107#include "MDataValue.h"
108#include "MDataMember.h"
109#include "MDataElement.h"
110
111ClassImp(MDataChain);
112
113// --------------------------------------------------------------------------
114//
115// Constructor which takes a rule and a surrounding operator as argument
116//
117MDataChain::MDataChain(const char *rule, OperatorType_t op)
118 : fOperatorType(op)
119{
120 fName = "MDataChain";
121 fTitle = rule;
122
123 fMember=ParseString(rule, 1);
124}
125
126// --------------------------------------------------------------------------
127//
128// Default constructor
129//
130MDataChain::MDataChain()
131 : fMember(NULL), fOperatorType(kENoop)
132{
133}
134
135// --------------------------------------------------------------------------
136//
137// Constructor taking a rule as an argument. For more details see
138// class description
139//
140MDataChain::MDataChain(const char *rule, const char *name, const char *title)
141 : fOperatorType(kENoop)
142{
143 fName = name ? name : "MDataChain";
144 fTitle = title ? title : rule;
145
146 *fLog << inf << "Trying to resolve rule... " << flush;
147 if (!(fMember=ParseString(rule, 1)))
148 {
149 *fLog << err << dbginf << "Parsing '" << rule << "' failed." << endl;
150 return;
151 }
152 *fLog << inf << "found: " << flush;
153 fMember->Print();
154 *fLog << endl;
155}
156
157// --------------------------------------------------------------------------
158//
159// PreProcesses all members in the list
160//
161Bool_t MDataChain::PreProcess(const MParList *pList)
162{
163 return fMember ? fMember->PreProcess(pList) : kFALSE;
164}
165
166// --------------------------------------------------------------------------
167//
168// Checks whether at least one member has the ready-to-save flag.
169//
170Bool_t MDataChain::IsReadyToSave() const
171{
172 *fLog << all << "fM=" << fMember << "/" << (int)fMember->IsReadyToSave() << " " << endl;
173 return fMember ? fMember->IsReadyToSave() : kFALSE;
174}
175
176// --------------------------------------------------------------------------
177//
178// Destructor. Delete filters.
179//
180MDataChain::~MDataChain()
181{
182 if (fMember)
183 delete fMember;
184}
185
186// --------------------------------------------------------------------------
187//
188// Returns the number of alphanumeric characters (including '.')
189// in the given string
190//
191Int_t MDataChain::IsAlNum(TString txt)
192{
193 int l = txt.Length();
194 for (int i=0; i<l; i++)
195 {
196 if (!isalnum(txt[i]) && txt[i]!='.' && /*txt[i]!='['&&txt[i]!=']'&&*/
197 ((txt[i]!='-' && txt[i]!='+') || i!=0))
198 return i;
199 }
200
201 return l;
202}
203
204Int_t MDataChain::GetBracket(TString txt, char open, char close)
205{
206 Int_t first=1;
207 for (int cnt=0; first<txt.Length(); first++)
208 {
209 if (txt[first]==open)
210 cnt++;
211 if (txt[first]==close)
212 cnt--;
213 if (cnt==-1)
214 break;
215 }
216 return first;
217}
218
219// --------------------------------------------------------------------------
220//
221// Compare a given text with the known operators. If no operator match
222// kENoop is retunred, otherwise the corresponding OperatorType_t
223//
224MDataChain::OperatorType_t MDataChain::ParseOperator(TString txt) const
225{
226 txt.ToLower();
227
228 if (txt=="abs") return kEAbs;
229 if (txt=="fabs") return kEAbs;
230 if (txt=="log") return kELog;
231 if (txt=="log10") return kELog10;
232 if (txt=="sin") return kESin;
233 if (txt=="cos") return kECos;
234 if (txt=="tan") return kETan;
235 if (txt=="sinh") return kESinH;
236 if (txt=="cosh") return kECosH;
237 if (txt=="tanh") return kETanH;
238 if (txt=="asin") return kEASin;
239 if (txt=="acos") return kEACos;
240 if (txt=="atan") return kEATan;
241 if (txt=="sqrt") return kESqrt;
242 if (txt=="sqr") return kESqr;
243 if (txt=="exp") return kEExp;
244 if (txt=="pow10") return kEPow10;
245 if (txt=="sgn") return kESgn;
246 if (txt=="floor") return kEFloor;
247 if (txt=="r2d") return kERad2Deg;
248 if (txt=="d2r") return kEDeg2Rad;
249 if (txt=="rand") return kERandom;
250 if (txt=="randp") return kERandomP;
251 if (txt=="rande") return kERandomE;
252 if (txt=="randi") return kERandomI;
253 if (txt=="randg") return kERandomG;
254 if (txt=="randl") return kERandomL;
255 if (txt[0]=='-') return kENegative;
256 if (txt[0]=='+') return kEPositive;
257
258 return kENoop;
259}
260
261// --------------------------------------------------------------------------
262//
263// Core of the data chain. Here the chain is constructed out of the rule.
264//
265MData *MDataChain::ParseString(TString txt, Int_t level)
266{
267 MData *member0=NULL;
268
269 char type=0;
270 int nlist = 0;
271
272 while (!txt.IsNull())
273 {
274 MData *newmember = NULL;
275
276 txt = txt.Strip(TString::kBoth);
277
278 switch (txt[0])
279 {
280 case '(':
281 {
282 //
283 // Search for the corresponding bracket
284 //
285 Int_t first=GetBracket(txt, '(', ')');
286
287 if (first==txt.Length())
288 {
289 *fLog << err << dbginf << "Syntax Error: ')' missing." << endl;
290 if (member0)
291 delete member0;
292 return NULL;
293 }
294
295 //
296 // Make a copy of the 'interieur' and delete the substringä
297 // including the brackets
298 //
299 TString sub = txt(1, first-1);
300 txt.Remove(0, first+1);
301
302 //
303 // Parse the substring
304 //
305 newmember = ParseString(sub, level+1);
306 if (!newmember)
307 {
308 *fLog << err << dbginf << "Parsing '" << sub << "' failed." << endl;
309 if (member0)
310 delete member0;
311 return NULL;
312 }
313 }
314 break;
315
316 case ')':
317 *fLog << err << dbginf << "Syntax Error: Too many ')'" << endl;
318 if (member0)
319 delete member0;
320 return NULL;
321
322 case '+':
323 case '-':
324 case '*':
325 case '/':
326 if (member0)
327 {
328 //
329 // Check for the type of the symbol
330 //
331 char is = txt[0];
332 txt.Remove(0, 1);
333
334 //
335 // If no filter is available or the available filter
336 // is of a different symbols we have to create a new
337 // data list with the new symbol
338 //
339 if (/*!member0->InheritsFrom(MDataMember::Class()) ||*/ type!=is)
340 {
341 MDataList *list = new MDataList(is);
342 list->SetName(Form("List_%c_%d", is, 10*level+nlist++));
343
344 list->SetOwner();
345 list->AddToList(member0);
346
347 member0 = list;
348
349 type = is;
350 }
351 continue;
352 }
353
354 if (txt[0]!='-' && txt[0]!='+')
355 {
356 *fLog << err << dbginf << "Syntax Error: First argument of symbol '";
357 *fLog << txt[0] << "' missing." << endl;
358 if (member0)
359 delete member0;
360 return NULL;
361 }
362
363 // FALLTHROU
364
365 case '0':
366 case '1':
367 case '2':
368 case '3':
369 case '4':
370 case '5':
371 case '6':
372 case '7':
373 case '8':
374 case '9':
375 if ((txt[0]!='-' && txt[0]!='+') || isdigit(txt[1]) || txt[1]=='.')
376 {
377 char *end;
378 Double_t num = strtod(txt.Data(), &end);
379 if (!end || txt.Data()==end)
380 {
381 *fLog << err << dbginf << "Error trying to convert '" << txt << "' to value." << endl;
382 if (member0)
383 delete member0;
384 return NULL;
385 }
386
387 txt.Remove(0, end-txt.Data());
388
389 newmember = new MDataValue(num);
390 break;
391 }
392
393 // FALLTHROUH
394
395 default:
396 int i = IsAlNum(txt);
397
398 if (i==0)
399 {
400 *fLog << err << dbginf << "Syntax Error: Name of data member missing in '" << txt << "'" << endl;
401 if (member0)
402 delete member0;
403 return NULL;
404 }
405
406 TString text = txt(0, i);
407
408 txt.Remove(0, i);
409 txt = txt.Strip(TString::kBoth);
410
411 if (!txt.IsNull() && txt[0]=='[')
412 {
413 Int_t first = GetBracket(txt, '[', ']');
414 TString op = txt(1, first-1);
415 txt.Remove(0, first+1);
416
417 newmember = new MDataElement(text, atoi(op));
418 break;
419 }
420 if ((txt.IsNull() || txt[0]!='(') && text[0]!='-' && text[0]!='+')
421 {
422 newmember = new MDataMember(text.Data());
423 break;
424 }
425
426 OperatorType_t op = ParseOperator(text);
427 if (op==kENoop)
428 {
429 *fLog << err << dbginf << "Syntax Error: Operator '" << text << "' unknown." << endl;
430 if (member0)
431 delete member0;
432 return NULL;
433 }
434
435 Int_t first = GetBracket(txt, '(', ')');
436 TString sub = op==kENegative || op==kEPositive ? text.Remove(0,1) + txt : txt(1, first-1);
437 txt.Remove(0, first+1);
438
439 newmember = new MDataChain(sub, op);
440 if (!newmember->IsValid())
441 {
442 *fLog << err << dbginf << "Syntax Error: Error parsing contents '" << sub << "' of operator " << text << endl;
443 delete newmember;
444 if (member0)
445 delete member0;
446 return NULL;
447 }
448 }
449
450 if (!member0)
451 {
452 member0 = newmember;
453 continue;
454 }
455
456 if (!member0->InheritsFrom(MDataList::Class()))
457 continue;
458
459 ((MDataList*)member0)->AddToList(newmember);
460 }
461
462 return member0;
463}
464
465// --------------------------------------------------------------------------
466//
467// Returns the value described by the rule
468//
469Double_t MDataChain::GetValue() const
470{
471 if (!fMember)
472 {
473 *fLog << warn << "MDataChain not valid." << endl;
474 return 0;
475 }
476
477 const Double_t val = fMember->GetValue();
478
479 switch (fOperatorType)
480 {
481 case kEAbs: return fabs(val);
482 case kELog: return log(val);
483 case kELog10: return log10(val);
484 case kESin: return sin(val);
485 case kECos: return cos(val);
486 case kETan: return tan(val);
487 case kESinH: return sinh(val);
488 case kECosH: return cosh(val);
489 case kETanH: return tanh(val);
490 case kEASin: return asin(val);
491 case kEACos: return acos(val);
492 case kEATan: return atan(val);
493 case kESqrt: return sqrt(val);
494 case kESqr: return val*val;
495 case kEExp: return exp(val);
496 case kEPow10: return pow(10, val);
497 case kESgn: return val<0 ? -1 : 1;
498 case kENegative: return -val;
499 case kEPositive: return val;
500 case kEFloor: return floor(val);
501 case kERad2Deg: return val*180/TMath::Pi();
502 case kEDeg2Rad: return val*TMath::Pi()/180;
503 case kERandom: return gRandom ? gRandom->Uniform(val) : 0;
504 case kERandomP: return gRandom ? gRandom->Poisson(val) : 0;
505 case kERandomE: return gRandom ? gRandom->Exp(val) : 0;
506 case kERandomI: return gRandom ? gRandom->Integer(val) : 0;
507 case kERandomG: return gRandom ? gRandom->Gaus(0, val) : 0;
508 case kERandomL: return gRandom ? gRandom->Landau(0, val) : 0;
509 case kENoop: return val;
510 }
511
512 *fLog << warn << "No Case for " << fOperatorType << " available." << endl;
513
514 return 0;
515}
516
517 /*
518void MDataChain::Print(Option_t *opt) const
519{
520 *fLog << GetRule() << flush;
521 Bool_t bracket = fOperatorType!=kENoop && !fMember->InheritsFrom(MDataList::Class());
522
523 switch (fOperatorType)
524 {
525 case kEAbs: *fLog << "abs" << flush; break;
526 case kELog: *fLog << "log" << flush; break;
527 case kELog10: *fLog << "log10" << flush; break;
528 case kESin: *fLog << "sin" << flush; break;
529 case kECos: *fLog << "cos" << flush; break;
530 case kETan: *fLog << "tan" << flush; break;
531 case kESinH: *fLog << "sinh" << flush; break;
532 case kECosH: *fLog << "cosh" << flush; break;
533 case kETanH: *fLog << "tanh" << flush; break;
534 case kEASin: *fLog << "asin" << flush; break;
535 case kEACos: *fLog << "acos" << flush; break;
536 case kEATan: *fLog << "atan" << flush; break;
537 case kESqrt: *fLog << "sqrt" << flush; break;
538 case kEExp: *fLog << "exp" << flush; break;
539 case kEPow10: *fLog << "pow10" << flush; break;
540 case kESgn: *fLog << "sgn" << flush; break;
541 case kENegative: *fLog << "-" << flush; break;
542 case kEPositive: *fLog << "+" << flush; break;
543 case kENoop:
544 break;
545 }
546
547 if (bracket)
548 *fLog << "(" << flush;
549
550 fMember->Print();
551
552 if (bracket)
553 *fLog << ")" << flush;
554 }
555 */
556
557// --------------------------------------------------------------------------
558//
559// Builds a rule from all the chain members. This is a rule which could
560// be used to rebuild the chain.
561//
562TString MDataChain::GetRule() const
563{
564 if (!fMember)
565 return "<n/a>";
566
567 TString str;
568
569 Bool_t bracket = fOperatorType!=kENoop && !fMember->InheritsFrom(MDataList::Class());
570
571 switch (fOperatorType)
572 {
573 case kEAbs: str += "abs" ; break;
574 case kELog: str += "log" ; break;
575 case kELog10: str += "log10"; break;
576 case kESin: str += "sin" ; break;
577 case kECos: str += "cos" ; break;
578 case kETan: str += "tan" ; break;
579 case kESinH: str += "sinh" ; break;
580 case kECosH: str += "cosh" ; break;
581 case kETanH: str += "tanh" ; break;
582 case kEASin: str += "asin" ; break;
583 case kEACos: str += "acos" ; break;
584 case kEATan: str += "atan" ; break;
585 case kESqrt: str += "sqrt" ; break;
586 case kESqr: str += "sqr" ; break;
587 case kEExp: str += "exp" ; break;
588 case kEPow10: str += "pow10"; break;
589 case kESgn: str += "sgn" ; break;
590 case kENegative: str += "-" ; break;
591 case kEPositive: str += "+" ; break;
592 case kEFloor: str += "floor"; break;
593 case kERad2Deg: str += "r2d" ; break;
594 case kEDeg2Rad: str += "d2r" ; break;
595 case kERandom: str += "rand" ; break;
596 case kERandomP: str += "randp"; break;
597 case kERandomE: str += "rande"; break;
598 case kERandomI: str += "randi"; break;
599 case kERandomG: str += "randg"; break;
600 case kERandomL: str += "randl"; break;
601 case kENoop:
602 break;
603 }
604
605 if (bracket)
606 str += "(";
607
608 str += fMember->GetRule();
609
610 if (bracket)
611 str += ")";
612
613 return str;
614}
615
616// --------------------------------------------------------------------------
617//
618// Return a comma seperated list of all data members used in the chain.
619// This is mainly used in MTask::AddToBranchList
620//
621TString MDataChain::GetDataMember() const
622{
623 return fMember->GetDataMember();
624}
625
Note: See TracBrowser for help on using the repository browser.