--- calo/flight/CaloProfile/src/CaloProfile.cpp	2009/06/10 13:00:23	1.4
+++ calo/flight/CaloProfile/src/CaloProfile.cpp	2010/05/13 13:55:37	1.18
@@ -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 
@@ -80,7 +281,9 @@
   PKT = 0;
   atime = 0;
   //
+  suf = "";
   debug = false;
+  usepl18x = false;
   //
 };
 
@@ -96,7 +299,9 @@
   PKT = 0;
   atime = 0;
   //
+  suf = "";
   debug = false;
+  usepl18x = false;
   //
 };
 
@@ -174,7 +379,7 @@
   //
   for (Int_t v=minv; v<maxv;v++){
      for (Int_t p=minp; p<maxp;p++){
-	TString hid = Form("clatv%ip%i",v,p);	
+       TString hid = Form("clatv%ip%i%s",v,p,suf.Data());	
         TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
 	if ( tc ){
 //	 tc->Clear();
@@ -182,7 +387,7 @@
 	 tc = new TCanvas(hid,hid);
 	};
 	//
-	TString thid = Form("hlatv%ip%i",v,p);	
+	TString thid = Form("hlatv%ip%i%s",v,p,suf.Data());	
         TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
 	if ( th ) th->Delete();
 //	 th->Clear();
@@ -225,7 +430,7 @@
   gStyle->SetNdivisions(510,"XY");
   //
   for (Int_t p=minp; p<maxp;p++){
-    TString hid = Form("c2dp%i",p);	
+    TString hid = Form("c2dp%i%s",p,suf.Data());	
     TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
     if ( tc ){
       //	 tc->Clear();
@@ -233,7 +438,7 @@
       tc = new TCanvas(hid,hid);
     };
     //
-    TString thid = Form("h2dp%i",p);	
+    TString thid = Form("h2dp%i%s",p,suf.Data());	
     TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
     if ( th ) th->Delete();
     //	 th->Clear();
@@ -324,6 +529,9 @@
   //
   for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
     mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
+    //
+    if ( !usepl18x && view1==0 && plane1==18 ) mip1 = 0.;
+    //
     estrip[view1][plane1][strip1] = mip1;
   };
   //
@@ -371,6 +579,9 @@
   //
   for (Int_t i=0; i<L2->GetCaloLevel1()->istrip ; i++){
     mip1 = L2->GetCaloLevel1()->DecodeEstrip(i,view1,plane1,strip1);
+    //
+    if ( !usepl18x && view1==0 && plane1==18 ) mip1 = 0.;
+    //
     es[view1][plane1][strip1] = mip1;
   };
   //
@@ -451,6 +662,8 @@
   PKT = 0;
   atime = 0;
   //
+  clp = L2->GetCaloLevel2();
+  //
   sel = true;
   cont = false;
   N = 0;
@@ -458,17 +671,64 @@
   mask18b = -1;
   //
   no18x = true;
+  usepl18x = !no18x;
   debug = false;
+  maskXE = false;
+  maskXO = false;
+  maskYE = false;
+  maskYO = false;
+  //
+  lmax = 0.;
+  umax = 100.;
+  slmax = "";
+  sumax = "";
+  suf = "";
+  xyaverage = true;
+  //
+  letmax = -1;
+  heavytail = false;
+  //  lmipth = 100.;
+  lmipth = 0.;
   //
 };
 
+TF1 *CaloLong::GetFit(){
+  TString fnam=Form("lfit%s",suf.Data());
+  TF1 *lfit  = dynamic_cast<TF1*>(gDirectory->FindObject(fnam));
+  if ( lfit ) return lfit;
+  return NULL;
+}
+
+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(){
   //
+  if ( debug ) printf(" Clear called \n");
+  //
   memset(eplane,0, 2*22*sizeof(Float_t));
   //
   chi2 = 0.;
   ndf = 0.;
   E0 = 0.;
+  defE0 = 0.;
   a = 0.;
   b = 0.;
   errE0 = 0.;
@@ -484,7 +744,7 @@
 
 void CaloLong::Print(){
   //
-  Process();
+  Fit();
   //
   printf("==================== Calorimeter Longitudinal Profile =======================\n");
   printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
@@ -494,6 +754,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 +862,26 @@
   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;
+    //
+    if ( !usepl18x && view==0 && plane==18 ) mip = 0.;
+    //
+    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];  
+  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
@@ -635,6 +909,106 @@
   //
 };
 
+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]);
+    };
+  };
+  //
+  if ( !clp && debug ) printf(" ERROR!! no CaloLevel2 obj!!\n");
+  //
+  // 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.) );
+  if ( debug ) printf(" X0pl %f tanx %f tany %f \n",X0pl,ytgx,ytgy);
+  //
+  // 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;
+  letmax = etmax;
+  //
+  if ( debug ) this->Print();
+  if ( debug ) printf(" exit \n");
+  //
+};
 
 Double_t ccurve(Double_t *ti,Double_t *par){
   //
@@ -654,6 +1028,44 @@
   this->Fit(false);
 };
 
+Float_t CaloLong::Evaluate(TString s, Float_t tmax, Float_t X0pl){
+  /* SAMPLE OUTPUT:
+     Enter Infix Expression : A + B + C / (E - F)
+     Postfix Expression is : A B + C E F - / +  
+  */  
+  if ( !s.Contains("tmax") && !s.Contains("X0pl") ){
+    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 h=Form("%f",X0pl);
+  TString *ts= new TString("");
+  ts->Prepend(s.Data());
+  ts->ReplaceAll("tmax",4,g.Data(),g.Capacity());
+  ts->ReplaceAll("X0pl",4,h.Data(),h.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();
@@ -687,7 +1099,7 @@
     gStyle->SetNdivisions(510,"XY");
   };
   //
-  TString hid = Form("clongfit");	
+  TString hid = Form("clongfit%s",suf.Data());	
   TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
   //  if ( tc ) tc->Delete();
   //  if ( tc ) tc->Close();
@@ -697,14 +1109,27 @@
     if ( tc ) tc->cd();
   };
   //
-  TString thid = Form("hlongfit");	
-  TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
+  TString thid = Form("hlongfit%s",suf.Data());	
+  TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
   if ( th ) th->Delete();
+  //
+  TString ghid = Form("glongfit%s",suf.Data());	
+  TGraphErrors *gh  = dynamic_cast<TGraphErrors*>(gDirectory->FindObject(ghid));
+  if ( gh ) gh->Delete();
+  //
   Float_t xpos = 0.;
   Float_t enemip = 0.;
   Float_t xmax = NC * X0pl + 0.2;
-  //  th = new TH1F(thid,thid,int(NC*1.5),-0.2,xmax);
-  th = new TH1F(thid,thid,100,-0.2,xmax);
+  //
+  Double_t xxx[50];
+  Double_t exx[50];
+  Double_t yyy[50];
+  Double_t eyy[50];
+  Int_t numpo = 0;
+  memset(xxx,0,50*sizeof(Double_t));
+  memset(exx,0,50*sizeof(Double_t));
+  memset(yyy,0,50*sizeof(Double_t));
+  memset(eyy,0,50*sizeof(Double_t));
   //
   // AGH, BUG!
   //
@@ -718,26 +1143,99 @@
     mmax = NC;
   };
   //
+  Float_t emax = 0.;
   Float_t qtotparz = 0.;
   for (Int_t st=mmin;st<mmax+1;st++){
     enemip = 0.;
     xpos = (st - mmin) * X0pl;
-    if ( st > mmin && st < mmax ){      
-      if ( no18x && ( st == 18+1 || st == mask18b+1 )){
-	enemip = 2. * eplane[1][st]; 
+    //
+    if ( xyaverage ){
+      //
+      if ( st > mmin && st < mmax ){      
+	if ( no18x && ( st == 18+1 || st == mask18b+1 )){
+	  if ( !maskYO ){
+	    enemip = 2. * eplane[1][st]; 
+	  } else {
+	    enemip = eplane[1][st]; 
+	  };
+	} else {
+	  enemip = eplane[0][st-1] + eplane[1][st]; 
+	};
       } else {
-	enemip = eplane[0][st-1] + eplane[1][st]; 
+	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;
+      if ( enemip > emax ) emax = enemip;
+      if ( enemip > 0. ){
+	xxx[numpo] = xpos;
+	exx[numpo] = 0.1;
+	yyy[numpo] = enemip;
+	eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
+	//if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
+	if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = sqrt(enemip)/5.;
+	if ( xpos > letmax-1 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 5.;
+	//if ( xpos > letmax-2 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 7.;
+	if ( xpos < 3. && heavytail ) eyy[numpo] /= 5.;
+	numpo++;
+	//	th->Fill(xpos,enemip);
+	if ( debug ) printf(" AVE Filling: st %i xpos %f energy %f error %f letmax %f lmpith %f \n",st,xpos,enemip,eyy[numpo-1],letmax,lmipth);
       };
     } else {
-      if ( st == mmin ) enemip = 2. * eplane[1][st]; 
-      if ( st == mmax ) enemip = 2. * eplane[0][st-1]; 
-    };
-    //
-    qtotparz += enemip;
-    if ( enemip > 0. ){
-      th->Fill(xpos,enemip);
-      if ( debug ) printf(" Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
+      for (Int_t jj=0; jj<2; jj++){
+	if ( st > mmin && st < mmax ){ 
+	  if ( jj == 0 && no18x && ( st == 18+1 || st == mask18b+1 )){
+	    enemip = eplane[1][st]; 
+	  } else {
+	    if ( jj == 0 ){
+	      enemip = eplane[jj][st-1];
+	    } else {
+	      enemip = eplane[jj][st]; 
+	    };	      
+	  };
+	} else {
+	  if ( st == mmin && jj == 1 ){
+	    enemip = eplane[jj][st]; 
+	  };
+	  if ( st == mmax && jj == 0){
+	    enemip = eplane[jj][st-1]; 
+	  };
+	};
+	//
+	qtotparz += enemip;
+      if ( enemip > emax ) emax = enemip;
+	if ( enemip > 0. ){
+	  xxx[numpo] = xpos;
+	  exx[numpo] = 0.1;
+	  yyy[numpo] = enemip;
+	  eyy[numpo] = sqrt(enemip*3.)+sqrt(5.);
+	  //if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
+	  //if ( xpos > letmax && enemip > lmipth && heavytail ) eyy[numpo] = (sqrt(enemip*3.)+sqrt(5.))/numpo;
+	  if ( xpos > letmax && enemip > lmipth && heavytail) eyy[numpo] = sqrt(enemip)/5.;
+	  if ( xpos > letmax-1 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 5.;
+	//if ( xpos > letmax-2 && xpos < letmax+1 && heavytail ) eyy[numpo] /= 7.;
+	  if ( xpos < 3. && heavytail ) eyy[numpo] /= 5.;	 
+	  //	  eyy[numpo] = sqrt(enemip)/(st*0.95);
+	  numpo++;
+	  //	  th->Fill(xpos,enemip);
+	  if ( debug ) printf(" XY Filling: st %i xpos %f energy %f \n",st,xpos,enemip);
+	};
+      };	
     };
+
     //
     //    for (Int_t v=1; v>=0;v--)// {
 //       //
@@ -763,10 +1261,24 @@
 //     };
   };
   //
-  TF1 *lfit = new TF1("lfit",ccurve,0.,xmax,3);
-  if ( debug ) printf("qtot %f qtotparz %f \n",L2->GetCaloLevel2()->qtot,qtotparz); 
+  //  th = new TH2F(thid,thid,int(NC*1.5),-0.2,xmax);
+  th = new TH2F(thid,thid,1000,-0.2,xmax,1000,0.,emax*1.2);
+  gh = new TGraphErrors(numpo,xxx,yyy,exx,eyy);
+  TString fnam=Form("lfit%s",suf.Data());
+  TF1 *lfit  = dynamic_cast<TF1*>(gDirectory->FindObject(fnam));
+  if ( lfit ) lfit->Delete();
+  lfit = new TF1(fnam,ccurve,0.,xmax,3);
+//   if ( !lfit ){
+//     lfit = new TF1(fnam,ccurve,0.,xmax,3);
+//   };
+//   lfit->Clear();
+//   lfit->SetName(fnam);
+//   lfit->SetRange(0.,xmax);
+//  lfit->Print();
+  //
+  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);
@@ -774,23 +1286,22 @@
   //  lfit->SetParLimits(0,0.,1000.);
   //  lfit->SetParLimits(1,-1.,80.);
   //  lfit->SetParLimits(2,-1.,10.);
-  TString optio;
+  // TString optio = "ROW"; // "RO"
+  TString optio = "RO"; // "RO"
   if ( debug ){
-    //    optio = "MERBOV";
-    //    optio = "MEROV";
-    //    optio = "EROV";
-    optio = "RNOV";
-    if ( draw ) optio = "ROV";
+    optio += "NV";
+    if ( draw ) optio += "V";
   } else {
-    //    optio = "MERNOQ";
-    //    optio = "ERNOQ";
-    optio = "RNOQ";
-    if ( draw ) optio = "ROQ";
+    optio += "NQ";
+    if ( draw ) optio += "Q";
   };
   //
-  if ( debug ) printf(" OK, start the fitting procedure...\n");
+  if ( debug ) printf(" OK, start the fitting procedure... options: %s \n",optio.Data());
   //
-  fitresult = th->Fit("lfit",optio);
+  //  fitresult = th->Fit("lfit",optio);
+  //  lfit->Update();
+  fitresult = gh->Fit("lfit",optio);
+  //  fitresult = gh->Fit("lfit","ROQW");
   //
   if ( debug ) printf(" the fit is done! result: %i \n",fitresult);
   //
@@ -809,8 +1320,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,X0pl);
+    if ( sumax.MaybeRegexp() ) umax = this->Evaluate(sumax,tmax,X0pl);
     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);
@@ -826,6 +1345,8 @@
     delete w;
     delete xp;
     delete wp;
+    delete xpp;
+    delete wpp;
     //    Float_t i10max = lfit->Integral(0.,10.*tmax);
     if ( debug ) printf(" Integral: %f \n",i10max); 
     //
@@ -859,7 +1380,12 @@
     // axis titles
     th->SetXTitle("Depth [X0]");
     th->SetYTitle("Energy [MIP]");
-    th->DrawCopy("Perror");
+    //    th->DrawCopy("Perror");
+    th->DrawCopy();
+    gh->SetMarkerSize(1.);
+    gh->SetMarkerStyle(21);
+    //    gh->SetMarkerColor(kRed);
+    gh->Draw("Psame");
     lfit->Draw("same");
     tc->Modified();
     tc->Update();
@@ -871,7 +1397,16 @@
     if ( th ) th->Delete();
   };
   //
-  delete lfit;
+//   gMinuit->SetPrintLevel(1);
+//   gMinuit->mnamin();
+//   //  gMinuit->mnfree(0);
+//   gMinuit->mncler();
+//   gMinuit->mnrset(1);
+  //  gMinuit->Delete();
+  //ROOT::Math::Minimizer::Clear();
+  //  gMinuit->mnhelp("*");
+  //
+  //  delete lfit;
   //
 };
 
@@ -900,7 +1435,7 @@
   //
   if ( maxv != -1 ){
     for (Int_t v=minv; v<maxv;v++){
-      TString hid = Form("clongv%i",v);	
+      TString hid = Form("clongv%i%s",v,suf.Data());	
       TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
       if ( tc ){
 	//	 tc->Clear();
@@ -908,7 +1443,7 @@
 	tc = new TCanvas(hid,hid);
       };
       //
-      TString thid = Form("hlongv%i",v);	
+      TString thid = Form("hlongv%i%s",v,suf.Data());	
       TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
       if ( th ) th->Delete();
       //	 th->Clear();
@@ -927,14 +1462,14 @@
     };
   } else {
     //
-    TString hid = Form("clongvyvx");	
+    TString hid = Form("clongvyvx%s",suf.Data());	
     TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
     if ( tc ){
     } else {
       tc = new TCanvas(hid,hid);
     };
     //
-    TString thid = Form("hlongvyvx");	
+    TString thid = Form("hlongvyvx%s",suf.Data());	
     TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
     if ( th ) th->Delete();
     th = new TH1F(thid,thid,44,-0.5,43.5);
@@ -962,3 +1497,23 @@
   Clear();
   //delete this;
 };
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+