--- calo/flight/CaloProfile/src/CaloProfile.cpp	2009/06/10 13:00:23	1.4
+++ 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 
@@ -451,6 +652,8 @@
   PKT = 0;
   atime = 0;
   //
+  clp = L2->GetCaloLevel2();
+  //
   sel = true;
   cont = false;
   N = 0;
@@ -459,9 +662,38 @@
   //
   no18x = true;
   debug = false;
+  maskXE = false;
+  maskXO = false;
+  maskYE = false;
+  maskYO = false;
+  //
+  lmax = 0.;
+  umax = 100.;
+  slmax = "";
+  sumax = "";
   //
 };
 
+void CaloLong::MaskSection(TString sec){
+  sec.ToUpper();
+  if ( sec.Contains("XO") ) maskXO = true;
+  if ( sec.Contains("YO") ) maskYO = true;
+  if ( sec.Contains("XE") ) maskXE = true;
+  if ( sec.Contains("YE") ) maskYE = true;
+}
+
+void CaloLong::UnMaskSections(){
+  this->UnMaskSection("XEXOYEYO");
+}
+
+void CaloLong::UnMaskSection(TString sec){
+  sec.ToUpper();
+  if ( sec.Contains("XO") ) maskXO = false;
+  if ( sec.Contains("YO") ) maskYO = false;
+  if ( sec.Contains("XE") ) maskXE = false;
+  if ( sec.Contains("YE") ) maskYE = false;
+}
+
 void CaloLong::Clear(){
   //
   memset(eplane,0, 2*22*sizeof(Float_t));
@@ -469,6 +701,7 @@
   chi2 = 0.;
   ndf = 0.;
   E0 = 0.;
+  defE0 = 0.;
   a = 0.;
   b = 0.;
   errE0 = 0.;
@@ -484,7 +717,7 @@
 
 void CaloLong::Print(){
   //
-  Process();
+  Fit();
   //
   printf("==================== Calorimeter Longitudinal Profile =======================\n");
   printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
@@ -494,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);
@@ -597,17 +835,23 @@
   Int_t plane = 0;
   Int_t strip = 0;
   Float_t mip = 0.;
+  Bool_t gof = true;
   for (Int_t i=0; i < L2->GetCaloLevel1()->istrip; i++){
     mip = L2->GetCaloLevel1()->DecodeEstrip(i,view,plane,strip);
-    eplane[view][plane] += mip;
+    gof = true;
+    if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
+    if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
+    if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
+    if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
+    if ( gof ) eplane[view][plane] += mip;
   };
   //
   // inclination factor (stolen from Daniele's code)
   //
   Float_t ytgx = 0;
   Float_t ytgy = 0;
-  ytgx = 0.76 * L2->GetCaloLevel2()->tanx[0];
-  ytgy = 0.76 * L2->GetCaloLevel2()->tany[0];  
+  ytgx = 0.76 * clp->tanx[0];
+  ytgy = 0.76 * clp->tany[0];  
   X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
   //
   // Find experimental plane of maximum
@@ -635,6 +879,102 @@
   //
 };
 
+void CaloLong::SetEnergies(Float_t myene[][22]){
+  //  
+  if ( !L2 ){
+    printf(" ERROR: cannot find PamLevel2 object, use the correct constructor or check your program!\n");
+    printf(" ERROR: CaloHough variables not filled \n");
+    return;
+  };
+  //
+  Bool_t newentry = false;
+  //
+  if ( L2->IsORB() ){
+    if ( L2->GetOrbitalInfo()->pkt_num != PKT || L2->GetOrbitalInfo()->OBT != OBT || L2->GetOrbitalInfo()->absTime != atime ){
+      newentry = true;
+      OBT = L2->GetOrbitalInfo()->OBT;
+      PKT = L2->GetOrbitalInfo()->pkt_num;
+      atime = L2->GetOrbitalInfo()->absTime;
+    };
+  } else {
+    newentry = true;
+  };
+  //
+  if ( !newentry ) return;
+  //
+  if ( debug ) printf(" SET ENERGIES Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
+  //
+  Clear();
+  //
+  // let's start
+  //
+  if ( cont ){
+    for (Int_t i=0; i<22; i++){
+      if ( i == (18+N) ){
+	mask18b =  18 + N;
+	break;
+      };
+    };
+  };  
+  //  
+  if ( sel ){
+    for (Int_t i=0; i<22; i++){
+      if ( i == (18-N) ){
+	mask18b =  18 - N;
+	break;
+      };
+    };
+  };
+  //
+  //  if ( mask18b == 18 ) mask18b = -1;
+  //
+  Int_t view = 0;
+  Int_t plane = 0;
+  Bool_t gof = true;
+  for (view=0; view < 2; view++){
+    for (plane=0; plane < 22; plane++){      
+      gof = true;
+      if ( maskXE && (plane%2)==0 && view==1 ) gof = false;
+      if ( maskXO && (plane%2)!=0 && view==1 ) gof = false;
+      if ( maskYE && (plane%2)!=0 && view==0 ) gof = false;
+      if ( maskYO && (plane%2)==0 && view==0 ) gof = false;
+      if ( gof ) eplane[view][plane] = myene[view][plane]; 
+      if ( debug ) printf(" XE %i XO %i YE %i YO %i eplane %i %i = %f ........... myene %i %i = %f\n", maskXE, maskXO, maskYE, maskYO,view,plane,eplane[view][plane],view,plane,myene[view][plane]);
+    };
+  };
+  //
+  // inclination factor (stolen from Daniele's code)
+  //
+  Float_t ytgx = 0;
+  Float_t ytgy = 0;
+  ytgx = 0.76 * clp->tanx[0];
+  ytgy = 0.76 * clp->tany[0];  
+  X0pl = sqrt( pow(0.76,2.) + pow(ytgx,2.) + pow(ytgy,2.) );
+  //
+  // Find experimental plane of maximum
+  //
+  Int_t pmax = 0;
+  Int_t vmax = 0;
+  Float_t emax = 0.;
+  for (Int_t v=0; v<2; v++){
+   for (Int_t i=0; i<22; i++){
+     if ( eplane[v][i] > emax ){
+       emax = eplane[v][i];
+       vmax = v;
+       pmax = i;
+     };
+   };
+  };
+  //
+  //
+  //
+  if ( vmax == 0 ) pmax++;
+  etmax = pmax * X0pl;
+  //
+  if ( debug ) this->Print();
+  if ( debug ) printf(" exit \n");
+  //
+};
 
 Double_t ccurve(Double_t *ti,Double_t *par){
   //
@@ -654,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();
@@ -724,13 +1100,29 @@
     xpos = (st - mmin) * X0pl;
     if ( st > mmin && st < mmax ){      
       if ( no18x && ( st == 18+1 || st == mask18b+1 )){
-	enemip = 2. * eplane[1][st]; 
+	if ( !maskYO ){
+	  enemip = 2. * eplane[1][st]; 
+	} else {
+	  enemip = eplane[1][st]; 
+	};
       } else {
 	enemip = eplane[0][st-1] + eplane[1][st]; 
       };
     } else {
-      if ( st == mmin ) enemip = 2. * eplane[1][st]; 
-      if ( st == mmax ) enemip = 2. * eplane[0][st-1]; 
+      if ( st == mmin ){
+	if ( !maskYE ){
+	  enemip = 2. * eplane[1][st]; 
+	} else {
+	  enemip = eplane[1][st]; 
+	};
+      };
+      if ( st == mmax ){
+	if ( !maskXE ){
+	  enemip = 2. * eplane[0][st-1]; 
+	} else {
+	  enemip = eplane[0][st-1]; 
+	};
+      };
     };
     //
     qtotparz += enemip;
@@ -764,9 +1156,9 @@
   };
   //
   TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
-  if ( debug ) printf("qtot %f qtotparz %f \n",L2->GetCaloLevel2()->qtot,qtotparz); 
+  if ( debug ) printf("qtot %f qtotparz %f \n",clp->qtot,qtotparz); 
   E0 = qtotparz;
-  //  E0 = L2->GetCaloLevel2()->qtot;
+  //  E0 = clp->qtot;
   a = 5.;
   b = 0.5;
   if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
@@ -809,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);
@@ -962,3 +1362,23 @@
   Clear();
   //delete this;
 };
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+