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

Last change on this file since 2313 was 2173, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 19.9 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//
85// Constants are implemented in ParseDataMember, namely:
86// kPi: TMath::Pi()
87// kRad2Deg: 180/kPi
88// kDeg2Rad: kPi/180
89//
90// You can also defined constants which are defined in TMath by:
91// kLn10 for static Double_t TMath::Ln10();
92// kLogE for static Double_t TMath::LogE();
93// kRadToDeg for static Double_t TMath::RadToDeg();
94// kDegToRad for static Double_t TMath::DegToRad();
95// Remark: In older root versions only Pi() and E() are implemented in
96// TMath.
97//
98//
99// REMARK:
100// - All the random functions are returning 0 if gRandom==0
101// - You may get better results if you are using a TRandom3
102//
103// FIXME: The possibility to use other objects inheriting from MData
104// is missing.
105// Maybe we can use gInterpreter->Calc("") for this.
106// gROOT->ProcessLineFast("line");
107//
108/////////////////////////////////////////////////////////////////////////////
109
110#include "MDataChain.h"
111
112#include <math.h> // fabs on Alpha
113#include <ctype.h> // isalnum, ...
114#include <stdlib.h> // strtod, ...
115
116#include <TRandom.h>
117#include <TMethodCall.h>
118
119#include "MLog.h"
120#include "MLogManip.h"
121
122#include "MDataList.h"
123#include "MDataValue.h"
124#include "MDataMember.h"
125#include "MDataElement.h"
126
127ClassImp(MDataChain);
128
129using namespace std;
130
131// --------------------------------------------------------------------------
132//
133// Constructor which takes a rule and a surrounding operator as argument
134//
135MDataChain::MDataChain(const char *rule, OperatorType_t op)
136 : fOperatorType(op)
137{
138 fName = "MDataChain";
139 fTitle = rule;
140
141 fMember=ParseString(rule, 1);
142}
143
144// --------------------------------------------------------------------------
145//
146// Constructor taking a rule as an argument. For more details see
147// class description
148//
149MDataChain::MDataChain(const char *rule, const char *name, const char *title)
150 : fMember(NULL), fOperatorType(kENoop)
151{
152 fName = name ? name : "MDataChain";
153 fTitle = title ? title : rule;
154
155 if (TString(rule).IsNull())
156 return;
157
158 *fLog << inf << "Trying to resolve rule... " << flush;
159 if (!(fMember=ParseString(rule, 1)))
160 {
161 *fLog << err << dbginf << "Parsing '" << rule << "' failed." << endl;
162 return;
163 }
164 *fLog << inf << "found: " << GetRule() << endl;
165}
166
167// --------------------------------------------------------------------------
168//
169// PreProcesses all members in the list
170//
171Bool_t MDataChain::PreProcess(const MParList *pList)
172{
173 return fMember ? fMember->PreProcess(pList) : kFALSE;
174}
175
176// --------------------------------------------------------------------------
177//
178// Checks whether at least one member has the ready-to-save flag.
179//
180Bool_t MDataChain::IsReadyToSave() const
181{
182 *fLog << all << "fM=" << fMember << "/" << (int)fMember->IsReadyToSave() << " " << endl;
183 return fMember ? fMember->IsReadyToSave() : kFALSE;
184}
185
186// --------------------------------------------------------------------------
187//
188// Destructor. Delete filters.
189//
190MDataChain::~MDataChain()
191{
192 if (fMember)
193 delete fMember;
194}
195
196// --------------------------------------------------------------------------
197//
198// Returns the number of alphanumeric characters (including '.')
199// in the given string
200//
201Int_t MDataChain::IsAlNum(TString txt)
202{
203 int l = txt.Length();
204 for (int i=0; i<l; i++)
205 {
206 if (!isalnum(txt[i]) && txt[i]!='.' && /*txt[i]!='['&&txt[i]!=']'&&*/
207 ((txt[i]!='-' && txt[i]!='+') || i!=0))
208 return i;
209 }
210
211 return l;
212}
213
214Int_t MDataChain::GetBracket(TString txt, char open, char close)
215{
216 Int_t first=1;
217 for (int cnt=0; first<txt.Length(); first++)
218 {
219 if (txt[first]==open)
220 cnt++;
221 if (txt[first]==close)
222 cnt--;
223 if (cnt==-1)
224 break;
225 }
226 return first;
227}
228
229// --------------------------------------------------------------------------
230//
231// Compare a given text with the known operators. If no operator match
232// kENoop is retunred, otherwise the corresponding OperatorType_t
233//
234MDataChain::OperatorType_t MDataChain::ParseOperator(TString txt) const
235{
236 txt.ToLower();
237
238 if (txt=="abs") return kEAbs;
239 if (txt=="fabs") return kEAbs;
240 if (txt=="log") return kELog;
241 if (txt=="log10") return kELog10;
242 if (txt=="sin") return kESin;
243 if (txt=="cos") return kECos;
244 if (txt=="tan") return kETan;
245 if (txt=="sinh") return kESinH;
246 if (txt=="cosh") return kECosH;
247 if (txt=="tanh") return kETanH;
248 if (txt=="asin") return kEASin;
249 if (txt=="acos") return kEACos;
250 if (txt=="atan") return kEATan;
251 if (txt=="sqrt") return kESqrt;
252 if (txt=="sqr") return kESqr;
253 if (txt=="exp") return kEExp;
254 if (txt=="pow10") return kEPow10;
255 if (txt=="sgn") return kESgn;
256 if (txt=="floor") return kEFloor;
257 if (txt=="r2d") return kERad2Deg;
258 if (txt=="d2r") return kEDeg2Rad;
259 if (txt=="rand") return kERandom;
260 if (txt=="randp") return kERandomP;
261 if (txt=="rande") return kERandomE;
262 if (txt=="randi") return kERandomI;
263 if (txt=="randg") return kERandomG;
264 if (txt=="randl") return kERandomL;
265 if (txt[0]=='-') return kENegative;
266 if (txt[0]=='+') return kEPositive;
267
268 return kENoop;
269}
270
271// --------------------------------------------------------------------------
272//
273// Here the names of data members are interpreted. This can be used to
274// check for constants.
275//
276MData *MDataChain::ParseDataMember(TString txt)
277{
278 //txt.ToLower();
279
280 if (txt=="kRad2Deg") return new MDataValue(kRad2Deg);
281 if (txt=="kDeg2Rad") return new MDataValue(1./kRad2Deg);
282
283 if (!txt.BeginsWith("k"))
284 return new MDataMember(txt.Data());
285
286 const TString name = txt(1, txt.Length());
287 TMethodCall call(TMath::Class(), name, "");
288 switch (call.ReturnType())
289 {
290 case TMethodCall::kLong:
291 Long_t l;
292 call.Execute(l);
293 return new MDataValue(l);
294 case TMethodCall::kDouble:
295 Double_t d;
296 call.Execute(d);
297 return new MDataValue(d);
298 default:
299 break;
300 }
301
302 return new MDataMember(txt.Data());
303}
304
305// --------------------------------------------------------------------------
306//
307// Core of the data chain. Here the chain is constructed out of the rule.
308//
309MData *MDataChain::ParseString(TString txt, Int_t level)
310{
311 MData *member0=NULL;
312
313 char type=0;
314 int nlist = 0;
315
316 while (!txt.IsNull())
317 {
318 MData *newmember = NULL;
319
320 txt = txt.Strip(TString::kBoth);
321
322 switch (txt[0])
323 {
324 case '(':
325 {
326 //
327 // Search for the corresponding bracket
328 //
329 Int_t first=GetBracket(txt, '(', ')');
330
331 if (first==txt.Length())
332 {
333 *fLog << err << dbginf << "Syntax Error: ')' missing." << endl;
334 if (member0)
335 delete member0;
336 return NULL;
337 }
338
339 //
340 // Make a copy of the 'interieur' and delete the substringä
341 // including the brackets
342 //
343 TString sub = txt(1, first-1);
344 txt.Remove(0, first+1);
345
346 //
347 // Parse the substring
348 //
349 newmember = ParseString(sub, level+1);
350 if (!newmember)
351 {
352 *fLog << err << dbginf << "Parsing '" << sub << "' failed." << endl;
353 if (member0)
354 delete member0;
355 return NULL;
356 }
357 }
358 break;
359
360 case ')':
361 *fLog << err << dbginf << "Syntax Error: Too many ')'" << endl;
362 if (member0)
363 delete member0;
364 return NULL;
365
366 case '+':
367 case '-':
368 case '*':
369 case '/':
370 if (member0)
371 {
372 //
373 // Check for the type of the symbol
374 //
375 char is = txt[0];
376 txt.Remove(0, 1);
377
378 //
379 // If no filter is available or the available filter
380 // is of a different symbols we have to create a new
381 // data list with the new symbol
382 //
383 if (/*!member0->InheritsFrom(MDataMember::Class()) ||*/ type!=is)
384 {
385 MDataList *list = new MDataList(is);
386 list->SetName(Form("List_%c_%d", is, 10*level+nlist++));
387
388 list->SetOwner();
389 list->AddToList(member0);
390
391 member0 = list;
392
393 type = is;
394 }
395 continue;
396 }
397
398 if (txt[0]!='-' && txt[0]!='+')
399 {
400 *fLog << err << dbginf << "Syntax Error: First argument of symbol '";
401 *fLog << txt[0] << "' missing." << endl;
402 if (member0)
403 delete member0;
404 return NULL;
405 }
406
407 // FALLTHROU
408
409 case '0':
410 case '1':
411 case '2':
412 case '3':
413 case '4':
414 case '5':
415 case '6':
416 case '7':
417 case '8':
418 case '9':
419 if ((txt[0]!='-' && txt[0]!='+') || isdigit(txt[1]) || txt[1]=='.')
420 {
421 char *end;
422 Double_t num = strtod(txt.Data(), &end);
423 if (!end || txt.Data()==end)
424 {
425 *fLog << err << dbginf << "Error trying to convert '" << txt << "' to value." << endl;
426 if (member0)
427 delete member0;
428 return NULL;
429 }
430
431 txt.Remove(0, end-txt.Data());
432
433 newmember = new MDataValue(num);
434 break;
435 }
436
437 // FALLTHROUH
438
439 default:
440 int i = IsAlNum(txt);
441
442 if (i==0)
443 {
444 *fLog << err << dbginf << "Syntax Error: Name of data member missing in '" << txt << "'" << endl;
445 if (member0)
446 delete member0;
447 return NULL;
448 }
449
450 TString text = txt(0, i);
451
452 txt.Remove(0, i);
453 txt = txt.Strip(TString::kBoth);
454
455 if (!txt.IsNull() && txt[0]=='[')
456 {
457 Int_t first = GetBracket(txt, '[', ']');
458 TString op = txt(1, first-1);
459 txt.Remove(0, first+1);
460
461 newmember = new MDataElement(text, atoi(op));
462 break;
463 }
464 if ((txt.IsNull() || txt[0]!='(') && text[0]!='-' && text[0]!='+')
465 {
466 newmember = ParseDataMember(text.Data());
467 break;
468 }
469
470 OperatorType_t op = ParseOperator(text);
471 if (op==kENoop)
472 {
473 *fLog << err << dbginf << "Syntax Error: Operator '" << text << "' unknown." << endl;
474 if (member0)
475 delete member0;
476 return NULL;
477 }
478
479 Int_t first = GetBracket(txt, '(', ')');
480 TString sub = op==kENegative || op==kEPositive ? text.Remove(0,1) + txt : txt(1, first-1);
481 txt.Remove(0, first+1);
482
483 newmember = new MDataChain(sub, op);
484 if (!newmember->IsValid())
485 {
486 *fLog << err << dbginf << "Syntax Error: Error parsing contents '" << sub << "' of operator " << text << endl;
487 delete newmember;
488 if (member0)
489 delete member0;
490 return NULL;
491 }
492 }
493
494 if (!member0)
495 {
496 member0 = newmember;
497 continue;
498 }
499
500 if (!member0->InheritsFrom(MDataList::Class()))
501 continue;
502
503 ((MDataList*)member0)->AddToList(newmember);
504 }
505
506 return member0;
507}
508
509// --------------------------------------------------------------------------
510//
511// Returns the value described by the rule
512//
513Double_t MDataChain::GetValue() const
514{
515 if (!fMember)
516 {
517 //*fLog << warn << "MDataChain not valid." << endl;
518 return 0;
519 }
520
521 const Double_t val = fMember->GetValue();
522
523 switch (fOperatorType)
524 {
525 case kEAbs: return fabs(val);
526 case kELog: return log(val);
527 case kELog10: return log10(val);
528 case kESin: return sin(val);
529 case kECos: return cos(val);
530 case kETan: return tan(val);
531 case kESinH: return sinh(val);
532 case kECosH: return cosh(val);
533 case kETanH: return tanh(val);
534 case kEASin: return asin(val);
535 case kEACos: return acos(val);
536 case kEATan: return atan(val);
537 case kESqrt: return sqrt(val);
538 case kESqr: return val*val;
539 case kEExp: return exp(val);
540 case kEPow10: return pow(10, val);
541 case kESgn: return val<0 ? -1 : 1;
542 case kENegative: return -val;
543 case kEPositive: return val;
544 case kEFloor: return floor(val);
545 case kERad2Deg: return val*180/TMath::Pi();
546 case kEDeg2Rad: return val*TMath::Pi()/180;
547 case kERandom: return gRandom ? gRandom->Uniform(val) : 0;
548 case kERandomP: return gRandom ? gRandom->Poisson(val) : 0;
549 case kERandomE: return gRandom ? gRandom->Exp(val) : 0;
550 case kERandomI: return gRandom ? gRandom->Integer((int)val) : 0;
551 case kERandomG: return gRandom ? gRandom->Gaus(0, val) : 0;
552 case kERandomL: return gRandom ? gRandom->Landau(0, val) : 0;
553 case kENoop: return val;
554 }
555
556 *fLog << warn << "No Case for " << fOperatorType << " available." << endl;
557
558 return 0;
559}
560
561 /*
562void MDataChain::Print(Option_t *opt) const
563{
564 *fLog << GetRule() << flush;
565 Bool_t bracket = fOperatorType!=kENoop && !fMember->InheritsFrom(MDataList::Class());
566
567 switch (fOperatorType)
568 {
569 case kEAbs: *fLog << "abs" << flush; break;
570 case kELog: *fLog << "log" << flush; break;
571 case kELog10: *fLog << "log10" << flush; break;
572 case kESin: *fLog << "sin" << flush; break;
573 case kECos: *fLog << "cos" << flush; break;
574 case kETan: *fLog << "tan" << flush; break;
575 case kESinH: *fLog << "sinh" << flush; break;
576 case kECosH: *fLog << "cosh" << flush; break;
577 case kETanH: *fLog << "tanh" << flush; break;
578 case kEASin: *fLog << "asin" << flush; break;
579 case kEACos: *fLog << "acos" << flush; break;
580 case kEATan: *fLog << "atan" << flush; break;
581 case kESqrt: *fLog << "sqrt" << flush; break;
582 case kEExp: *fLog << "exp" << flush; break;
583 case kEPow10: *fLog << "pow10" << flush; break;
584 case kESgn: *fLog << "sgn" << flush; break;
585 case kENegative: *fLog << "-" << flush; break;
586 case kEPositive: *fLog << "+" << flush; break;
587 case kENoop:
588 break;
589 }
590
591 if (bracket)
592 *fLog << "(" << flush;
593
594 fMember->Print();
595
596 if (bracket)
597 *fLog << ")" << flush;
598 }
599 */
600
601// --------------------------------------------------------------------------
602//
603// Builds a rule from all the chain members. This is a rule which could
604// be used to rebuild the chain.
605//
606TString MDataChain::GetRule() const
607{
608 if (!fMember)
609 return "<n/a>";
610
611 TString str;
612
613 Bool_t bracket = fOperatorType!=kENoop && !fMember->InheritsFrom(MDataList::Class());
614
615 switch (fOperatorType)
616 {
617 case kEAbs: str += "abs" ; break;
618 case kELog: str += "log" ; break;
619 case kELog10: str += "log10"; break;
620 case kESin: str += "sin" ; break;
621 case kECos: str += "cos" ; break;
622 case kETan: str += "tan" ; break;
623 case kESinH: str += "sinh" ; break;
624 case kECosH: str += "cosh" ; break;
625 case kETanH: str += "tanh" ; break;
626 case kEASin: str += "asin" ; break;
627 case kEACos: str += "acos" ; break;
628 case kEATan: str += "atan" ; break;
629 case kESqrt: str += "sqrt" ; break;
630 case kESqr: str += "sqr" ; break;
631 case kEExp: str += "exp" ; break;
632 case kEPow10: str += "pow10"; break;
633 case kESgn: str += "sgn" ; break;
634 case kENegative: str += "-" ; break;
635 case kEPositive: str += "+" ; break;
636 case kEFloor: str += "floor"; break;
637 case kERad2Deg: str += "r2d" ; break;
638 case kEDeg2Rad: str += "d2r" ; break;
639 case kERandom: str += "rand" ; break;
640 case kERandomP: str += "randp"; break;
641 case kERandomE: str += "rande"; break;
642 case kERandomI: str += "randi"; break;
643 case kERandomG: str += "randg"; break;
644 case kERandomL: str += "randl"; break;
645 case kENoop:
646 break;
647 }
648
649 if (bracket)
650 str += "(";
651
652 str += fMember->GetRule();
653
654 if (bracket)
655 str += ")";
656
657 return str;
658}
659
660// --------------------------------------------------------------------------
661//
662// Return a comma seperated list of all data members used in the chain.
663// This is mainly used in MTask::AddToBranchList
664//
665TString MDataChain::GetDataMember() const
666{
667 return fMember->GetDataMember();
668}
669
Note: See TracBrowser for help on using the repository browser.