--- calo/flight/CaloProfile/src/CaloProfile.cpp 2009/08/11 14:23:09 1.8 +++ calo/flight/CaloProfile/src/CaloProfile.cpp 2009/08/12 14:57:00 1.10 @@ -4,6 +4,207 @@ ClassImp(CaloLong); ClassImp(Calo2D); + +int isempty(struct stack *s){ + return (s->top == EMPTY) ? 1 : 0; +} + +void emptystack(struct stack* s){ + s->top=EMPTY; + strcpy(s->data," "); +} + +void push(struct stack* s,int item){ + if(s->top == (MAX-1)){ + printf("\n ERROR! STACK FULL (too many digits in SetLower/UpperLimit)\n"); + } else { + s->data[++s->top]=(char)item; + }; +} + +char pop(struct stack* s){ + char ret=(char)EMPTY; + if(!isempty(s)){ + ret= s->data[s->top--]; + }; + return ret; +} + + +int ipop(struct stack* s){ + int ret=EMPTY; + if(s->top == EMPTY) + printf("\n ERROR! STACK EMPTY (too few digits in SetLower/UpperLimit)\n"); + else { + ret= s->data[s->top--]; + }; + return ret; +} + + +void display(struct stack s){ + int ss = s.top; + while(s.top != EMPTY){ + printf("\n%d",s.data[s.top--]); + }; + printf(" s.top %i \n",ss); +} + +int isoperator(char e){ + if(e == '+' || e == '-' || e == '*' || e == '/' || e == '%') + return 1; + else + return 0; +} + + +int priority(char e){ + int pri = 0; + if(e == '*' || e == '/' || e =='%') + pri = 2; + else { + if(e == '+' || e == '-') pri = 1; + }; + return pri; +} + +void infix2postfix(char* infix, char * postfix, int insertspace){ + char *i,*p; + struct stack X; + char n1; + emptystack(&X); + i = &infix[0]; + p = &postfix[0]; + while(*i){ + while(*i == ' ' || *i == '\t' ){ + i++; + }; + TString c=i; + if( isdigit(*i) || isalpha(*i) || c.BeginsWith(".") ){ + while( isdigit(*i) || isalpha(*i) || c.BeginsWith(".") ){ + *p = *i; + p++; + i++; + c=i; + }; + /*SPACE CODE*/ + if(insertspace){ + *p = ' '; + p++; + }; + /*END SPACE CODE*/ + }; + if( *i == '(' ){ + push(&X,*i); + i++; + }; + if( *i == ')'){ + n1 = pop(&X); + while( n1 != '(' ){ + *p = n1; + p++; + /*SPACE CODE*/ + if(insertspace){ + *p = ' '; + p++; + }; + /*END SPACE CODE*/ + n1 = pop(&X); + }; + i++; + }; + if( isoperator(*i) ){ + if(isempty(&X)){ + push(&X,*i); + } else { + n1 = pop(&X); + while( priority(n1) >= priority(*i) ){ + *p = n1; + p++; + /*SPACE CODE*/ + if(insertspace){ + *p = ' '; + p++; + }; + /*END SPACE CODE*/ + n1 = pop(&X); + }; + push(&X,n1); + push(&X,*i); + }; + i++; + }; + }; + while(!isempty(&X)){ + n1 = pop(&X); + *p = n1; + p++; + /*SPACE CODE*/ + if(insertspace){ + *p = ' '; + p++; + }; + /*END SPACE CODE*/ + }; + *p = '\0'; +} + + +Float_t evaluate(char *postfix){ + // + Float_t op1 = 0.; + Float_t op2 = 0.; + Float_t result = 0.; + // + TString e = postfix; + Float_t st[50]; + memset(st,0,50*sizeof(Float_t)); + TObjArray *ae = e.Tokenize(" "); + Int_t o = 0; + Int_t a = 0; + Int_t ap = 0; + // + while ( o < ae->GetEntries() ){ + if ( ((TString)(((TObjString*)ae->At(o))->GetString())).IsFloat() ){ + st[a] = (Float_t)((TString)(((TObjString*)ae->At(o))->GetString())).Atof();; + a++; + } else { + ap = a-1; + op1 = st[ap--]; + op2 = st[ap]; + const char *p=((TString)(((TObjString*)ae->At(o))->GetString())).Data(); + switch(*p){ + case '+': + result = op2 + op1; + break; + case '-': + result = op2 - op1; + break; + case '/': + result = op2 / op1; + break; + case '*': + result = op2 * op1; + break; + case '%': + result = (Int_t)round(op2) % (Int_t)round(op1); + break; + default: + printf("\nInvalid Operator: %s \n",((TString)(((TObjString*)ae->At(o))->GetString())).Data()); + return 0; + }; + st[ap] = result; + a = ap+1; + // + }; + o++; + }; + return result; + // +} + + + //////////////////////////////////////////////////////////////////////// /** * 1-dimension function describing lateral distribution of the @@ -466,6 +667,11 @@ maskYE = false; maskYO = false; // + lmax = 0.; + umax = 100.; + slmax = ""; + sumax = ""; + // }; void CaloLong::MaskSection(TString sec){ @@ -495,6 +701,7 @@ chi2 = 0.; ndf = 0.; E0 = 0.; + defE0 = 0.; a = 0.; b = 0.; errE0 = 0.; @@ -520,6 +727,11 @@ printf(" nchi2 :.. %f\n",chi2/ndf); printf(" E0 :.. %f\n",E0); printf(" E0/260. :.. %f\n",E0/260.); + printf(" defE0 :.. %f\n",defE0); + printf(" lower :.. %f\n",lmax); + printf(" upper :.. %f\n",umax); + printf(" s lower :.. %s\n",slmax.Data()); + printf(" s upper :.. %s\n",sumax.Data()); printf(" a :.. %f\n",a); printf(" b :.. %f\n",b); printf(" errE0 :.. %f\n",errE0); @@ -782,6 +994,42 @@ this->Fit(false); }; +Float_t CaloLong::Evaluate(TString s, Float_t tmax){ + /* SAMPLE OUTPUT: + Enter Infix Expression : A + B + C / (E - F) + Postfix Expression is : A B + C E F - / + + */ + if ( !s.Contains("t") ){ + printf(" ERROR, the input formula must contain \"t\"\n"); + return 0.; + }; + if ( tmax != tmax ){ + printf(" ERROR, tmax is nan! \n"); + return 0.; + }; + TString g=Form("%f",tmax); + TString *ts= new TString(""); + ts->Prepend(s.Data()); + ts->ReplaceAll("t",1,g.Data(),g.Capacity()); + ts->Prepend("("); + ts->Append(")"); + if ( debug ) printf(" ts %s tssize %i capac %i s %s g %s \n",ts->Data(),ts->Sizeof(),ts->Capacity(),s.Data(),g.Data()); + char in[50],post[50]; + strcpy(&in[0]," "); + strcpy(&in[0],ts->Data()); + strcpy(&post[0]," "); + if ( debug ) printf("Infix Expression is : ##%s##\n",&in[0]); + infix2postfix(&in[0],&post[0],1); + if ( debug ) printf("Postfix Expression is : ##%s##\n",&post[0]); + /* SAMPLE OUTPUT: + Enter Postfix Expression : 3 5 + 2 / + 3 5 + 2 / EQUALS 4 + */ + Float_t res = evaluate(&post[0]); + if ( debug ) printf("%s EQUALS %f\n",post,res); + return res; +} + void CaloLong::Fit(Bool_t draw){ // Process(); @@ -953,8 +1201,16 @@ if ( fitresult != 0 ){ if ( debug ) printf(" The fit failed, no integrals calculation and asymm is set to -1. \n"); asymm = -1.; + defE0 = -1.; } else { + if ( slmax.MaybeRegexp() ) lmax = this->Evaluate(slmax,tmax); + if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax); Int_t npp = 1000; + double *xpp=new double[npp]; + double *wpp=new double[npp]; + lfit->CalcGaussLegendreSamplingPoints(npp,xpp,wpp,1e-12); + defE0 = lfit->IntegralFast(npp,xpp,wpp,lmax,umax); + // double *xp=new double[npp]; double *wp=new double[npp]; lfit->CalcGaussLegendreSamplingPoints(npp,xp,wp,1e-12); @@ -1106,3 +1362,23 @@ Clear(); //delete this; }; + + + + + + + + + + + + + + + + + + + +