CaloProfile/src/CaloProfile.cpp

#include <CaloProfile.h>
//
ClassImp(CaloLat);
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 
 * shower as viewed by calorimeter
 * (projection of 3d function in one direction)
 *  
 *  xi[0] = x or y coordinate relative to shower axis
 *  parmin[0] = rt
 *  parmin[1] = p
 *  parmin[2] = rc
 *
 */
////////////////////////////////////////////////////////////////////////            
Double_t cfradx(Double_t *xi, Double_t *parmin) {

        double fradxmin2,p,rt,rc,es,x,pig,norm,c;
        x=*xi;
        pig = acos(-1.);
        rt=parmin[0];
        p=parmin[1];
        rc=parmin[2];
        norm=parmin[3];
        c=parmin[4];
        x=x-c;
        es=1.5;
        fradxmin2=p*pig*pow(rc,2)/pow((pow(x,2)+pow(rc,2)),es);
        fradxmin2=fradxmin2+(1-p)*pig*pow(rt,2)/pow((pow(x,2)+pow(rt,2)),es);
        fradxmin2=norm*fradxmin2/(2*pig);     
        //cout<<"x,fradxmin2 "<< x<<" "<<fradxmin2  <<endl;
  return fradxmin2;
}


Double_t cfunc(Double_t *tin,Double_t *par)
{

 Double_t gaa,a,tm,et,value,t,t0;
 t=*tin;
 et=par[0];
 a=par[1];
 tm=par[2];
 t0=par[3];
 gaa=TMath::Gamma(a);
 
 value=et*((a-1)/(tm*gaa))*pow(((a-1)*(t-t0)/tm),(a-1))*exp(-(a-1)*(t-t0)/tm);
 return value;
}


//--------------------------------------
/**
 * Default constructor 
 */
CaloLat::CaloLat(){
  Clear();
};

/**
 * Default constructor 
 */
Calo2D::Calo2D(){
  Clear();
};

CaloLat::CaloLat(PamLevel2 *l2p){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  //
  suf = "";
  debug = false;
  usepl18x = false;
  //
};

Calo2D::Calo2D(PamLevel2 *l2p){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  //
  suf = "";
  debug = false;
  usepl18x = false;
  //
};

void CaloLat::Clear(){
  //
};

void Calo2D::Clear(){
  //
};

void CaloLat::Print(){
  //
  Process();
  //
  printf("==================== Calorimeter Lateral Profile =======================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
//  printf(" nx [number of X combination]:.. %i\n",nx);
//  printf(" ny [number of Y combination]:.. %i\n",ny);
  printf("========================================================================\n");
  //
};

void Calo2D::Print(){
  //
  Process();
  //
  printf("==================== Calorimeter 2D Profile =======================\n");
  printf(" OBT: %u PKT: %u ATIME: %u \n",OBT,PKT,atime);
//  printf(" nx [number of X combination]:.. %i\n",nx);
//  printf(" ny [number of Y combination]:.. %i\n",ny);
  printf("========================================================================\n");
  //
};

void CaloLat::Draw(){
  //
  Process();
  Draw(-1,-1);
};

void Calo2D::Draw(){
  //
  Process();
  Draw(-1);
};

void CaloLat::Draw(Int_t view,Int_t plane){
  //
  Int_t minv = 0;
  Int_t maxv = 0;
  Int_t minp = 0;
  Int_t maxp = 0;
  //
  if ( view == -1 ){
        minv = 0;
        maxv = 2;
  } else {
        minv = view;
        maxv = view+1;
  };

  if ( plane == -1 ){
        minp = 0;
        maxp = 22;
  } else {
        minp = plane;
        maxp = plane+1;
  };
  // 
  Process();
  //
  gStyle->SetLabelSize(0.04);
  gStyle->SetNdivisions(510,"XY");
  //
  for (Int_t v=minv; v<maxv;v++){
     for (Int_t p=minp; p<maxp;p++){
       TString hid = Form("clatv%ip%i%s",v,p,suf.Data());       
        TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
        if ( tc ){
//       tc->Clear();
        } else {
         tc = new TCanvas(hid,hid);
        };
        //
        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();
//       th->Reset();
//      } else {
         th = new TH1F(thid,thid,96,-0.5,95.5);
//      };
        tc->cd();
        //
        for (Int_t st=0;st<96;st++){
         th->Fill(st,estrip[v][p][st]);
        };
        th->Draw();
        tc->Modified();
        tc->Update();
     };
  };
  //
  gStyle->SetLabelSize(0);
  gStyle->SetNdivisions(1,"XY");
  //
};

void Calo2D::Draw(Int_t plane){
  //
  Int_t minp = 0;
  Int_t maxp = 0;
  //
  if ( plane == -1 ){
        minp = 0;
        maxp = 23;
  } else {
        minp = plane;
        maxp = plane+1;
  };
  // 
  Process();
  //
  gStyle->SetLabelSize(0.04);
  gStyle->SetNdivisions(510,"XY");
  //
  for (Int_t p=minp; p<maxp;p++){
    TString hid = Form("c2dp%i%s",p,suf.Data());        
    TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
    if ( tc ){
      //         tc->Clear();
    } else {
      tc = new TCanvas(hid,hid);
    };
    //
    TString thid = Form("h2dp%i%s",p,suf.Data());       
    TH2F *th  = dynamic_cast<TH2F*>(gDirectory->FindObject(thid));
    if ( th ) th->Delete();
    //   th->Clear();
    //   th->Reset();
    //  } else {
    Int_t minx = smax[p] - 10;
    if ( minx < 0 ) minx = 0;
    Int_t maxx = minx + 20;
    if ( maxx > 95 ){
      maxx = 95;
      minx = 75;
    };
    Int_t miny = smay[p] - 10;
    if ( miny < 0 ) miny = 0;
    Int_t maxy = miny + 20;
    if ( maxy > 95 ){
      maxy = 95;
      miny = 75;
    };
    th = new TH2F(thid,thid,20,(Float_t)minx-0.5,(Float_t)maxx-0.5,20,(Float_t)miny-0.5,(Float_t)maxy-0.5);
    //    th = new TH2F(thid,thid,96,-0.5,95.5,96,-0.5,95.5);
    //  };
    tc->cd();
    //
    for (Int_t stx=minx;stx<maxx+1;stx++){
      for (Int_t sty=miny;sty<maxy+1;sty++){
        th->Fill(stx,sty,estrip[p][stx][sty]);
      };
    };
    gStyle->SetPalette(1);
    //    tc->SetLogz();    
    //    th->Draw("colbox");
    th->Draw("cont4");
    tc->Modified();
    tc->Update();
  };
  //
  gStyle->SetLabelSize(0);
  gStyle->SetNdivisions(1,"XY");
  //
};

void CaloLat::Delete(){
  Clear();
  //delete this;
};

void Calo2D::Delete(){
  Clear();
  //delete this;
};


void CaloLat::Process(){
  //  
  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(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();
  //
  // let's start
  //
  memset(estrip,0, 4224*sizeof(Float_t));
  Float_t mip1 = 0.;
  Int_t view1 = 0;
  Int_t plane1 = 0;
  Int_t strip1 = 0;
  //
  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;
  };
  //
  if ( debug ) this->Print();
  if ( debug ) printf(" exit \n");
  //
};

void Calo2D::Process(){
  //  
  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(" Start processing event at OBT %u PKT %u time %u \n",OBT,PKT,atime);
  //
  Clear();
  //
  // let's start
  //
  Float_t es[2][22][96];
  memset(es,0, 4224*sizeof(Float_t));
  memset(estrip,0, 4224*sizeof(Float_t));
  Float_t mip1 = 0.;
  Int_t view1 = 0;
  Int_t plane1 = 0;
  Int_t strip1 = 0;
  //
  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;
  };
  //
  Int_t plane2 = 0;
  Float_t emax[23];
  memset(emax,0,sizeof(Float_t)*23);
  memset(smax,0,sizeof(Int_t)*23);
  memset(smay,0,sizeof(Int_t)*23);
  //
  //
  for (Int_t p=0; p < 23 ; p++){
    //
    plane1 = p-1;
    plane2 = p;
    //
    if ( p == 0 ){
      plane1 = -1;
      plane2 = 0;
    };
    if ( p == 22 ){
      plane1 = 21;
      plane2 = -1;
    };
    //
    for (Int_t s=0; s < 96 ; s++){ // x
      for (Int_t ss=0; ss < 96 ; ss++){ // y
        if ( p == 0 ){
          estrip[p][s][ss] += es[1][plane2][ss];
          if ( (es[1][plane2][ss]) > emax[p] ){
            smax[p] = 45;
            smay[p] = ss;
            emax[p] = es[1][plane2][ss] ;
          };
        };
        if ( p > 0 && p < 22 ){
          estrip[p][s][ss] += es[0][plane1][s] + es[1][plane2][ss];
          if ( (es[0][plane1][s] + es[1][plane2][ss]) > emax[p] ){
            smax[p] = s;
            smay[p] = ss;
            emax[p] = es[0][plane1][s] + es[1][plane2][ss] ;
          };
        };      
        if ( p == 22 ){
          estrip[p][s][ss] += es[0][plane1][s];
          if ( (es[1][plane2][s]) > emax[p] ){
            smax[p] = s;
            smay[p] = 45;
            emax[p] = es[1][plane2][s] ;
          };
        };
      };
    };
    //
  };
  //
  if ( debug ) this->Print();
  if ( debug ) printf(" exit \n");
  //
};


/**
 * Default constructor 
 */
CaloLong::CaloLong(){
  Clear();
};

CaloLong::CaloLong(PamLevel2 *l2p){  
  //
  Clear();
  //
  L2 = l2p;
  //
  if ( !L2->IsORB() ) printf(" WARNING: OrbitalInfo Tree is needed, the plugin could not work properly without it \n");
  //
  OBT = 0;
  PKT = 0;
  atime = 0;
  //
  clp = L2->GetCaloLevel2();
  //
  sel = true;
  cont = false;
  N = 0;
  NC = 22;
  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.;
  //
  lfit = 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.;
  erra = 0.;
  errb = 0.;
  etmax = 0.;
  asymm = 0.;
  fitresult = 0;
  //
  X0pl = 0.76;
  //
};

void CaloLong::Print(){
  //
  Fit();
  //
  printf("==================== Calorimeter Longitudinal Profile =======================\n");
  printf(" OBT: %u PKT: %u ATIME: %u suf %s \n",OBT,PKT,atime,suf.Data());
  printf(" fitresult:.. %i\n",fitresult);
  printf(" chi2     :.. %f\n",chi2);
  printf(" ndf      :.. %f\n",ndf);
  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);
  printf(" erra     :.. %f\n",erra);
  printf(" errb     :.. %f\n",errb);
  printf(" asymm    :.. %f\n",asymm);
  printf(" tmax     :.. %f\n",((a-1.)/b));
  printf(" etmax    :.. %f\n",etmax);
  printf(" X0pl     :.. %f\n",X0pl);
  printf("========================================================================\n");
  //
};

void CaloLong::SetNoWpreSampler(Int_t n){
  //
  if ( NC+n <= 22 && NC+n >= 0 ){
    N = n;
  } else {
    printf(" ERROR! Calorimeter is made of 22 W planes\n");
    printf(" you are giving N presampler = %i and N calo = %i \n",n,NC);
    printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
    NC = 22;
    N = 0;
  };
}

void CaloLong::SetNoWcalo(Int_t n){
  if ( N+n <= 22 && N+n >= 0 ){
    NC = n;
  } else {
    printf(" ERROR! Calorimeter is made of 22 W planes\n");
    printf(" you are giving N W presampler = %i and N W calo = %i \n",N,n);
    printf(" WARNING: using default values NWpre = 0, NWcalo = 22\n");
    NC = 22;
    N = 0;
  };
}

void CaloLong::SplitInto(Int_t NoWpreSampler, Int_t NoWcalo){
  this->SetNoWpreSampler(0);
  this->SetNoWcalo(0);
  if ( NoWpreSampler < NoWcalo ){
          this->SetNoWpreSampler(NoWpreSampler);
          this->SetNoWcalo(NoWcalo);    
  } else {
          this->SetNoWcalo(NoWcalo);    
          this->SetNoWpreSampler(NoWpreSampler);
  };
}

void CaloLong::Process(){
  //  
  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(" 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;
  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);
    //
    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 * 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");
  //
};

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){
  //
  Double_t t = *ti;
  Double_t cE0 = par[0];
  Double_t ca = par[1];
  Double_t cb = par[2];
  Double_t gammaa = TMath::Gamma(ca); 
  //
  Double_t value = cE0 * cb * ( pow((cb * t),(ca - 1)) * exp( -cb * t ) ) / gammaa;
  //
  return value;
  //
}

void CaloLong::Fit(){
  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();
  //
  //  
  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 != fPKT || L2->GetOrbitalInfo()->OBT != fOBT || L2->GetOrbitalInfo()->absTime != fatime ){
      newentry = true;
      fOBT = L2->GetOrbitalInfo()->OBT;
      fPKT = L2->GetOrbitalInfo()->pkt_num;
      fatime = L2->GetOrbitalInfo()->absTime;
    };
  } else {
    newentry = true;
  };
  //
  if ( !newentry ) return;
  //
  if ( debug ) printf(" Start fitting event at OBT %u PKT %u time %u \n",fOBT,fPKT,fatime);
  //
  if ( draw ){
    gStyle->SetLabelSize(0.04);
    gStyle->SetNdivisions(510,"XY");
  };
  //
  TString hid = Form("clongfit%s",suf.Data());  
  TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
  //  if ( tc ) tc->Delete();
  //  if ( tc ) tc->Close();
  if ( !tc && draw ){
    tc = new TCanvas(hid,hid);
  } else {
    if ( tc ) tc->cd();
  };
  //
  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;
  //
  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!
  //
  Int_t mmin = 0;
  Int_t mmax = 0;
  if ( cont ){
    mmin = N;
    mmax = NC+N;
  } else {
    mmin = 0;
    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 ( 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 {
        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 {
      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--)// {
//       //
//       if ( v == 1 ){
//      xpos = (st - N) * X0pl;
//       } else {
//      xpos = (st + 1 - N) * X0pl;
//       };
//       //
//       if ( no18x && st == 18 && v == 0 ){
//      // skip plane 18x
//       } else {
//      if ( v == 1 && st == mask18b ){
//        // emulate plane 18x
//      } else {
//        if ( eplane[v][st] > 0. ){
//          th->Fill(xpos,eplane[v][st]);
//          if ( debug ) printf(" Filling: st %i v %i xpos %f energy %f \n",st,v,xpos,eplane[v][st]);
//        };
//      };
//       };
//       //
//     };
  };
  //
  //  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 = clp->qtot;
  a = 5.;
  b = 0.5;
  if ( debug ) printf(" STARTING PARAMETERS: E0 %f a %f b %f \n",E0,a,b);
  lfit->SetParameters(E0,a,b);
  //  lfit->SetParLimits(0,0.,1000.);
  //  lfit->SetParLimits(1,-1.,80.);
  //  lfit->SetParLimits(2,-1.,10.);
  // TString optio = "ROW"; // "RO"
  TString optio = "RO"; // "RO"
  if ( debug ){
    optio += "NV";
    if ( draw ) optio += "V";
  } else {
    optio += "NQ";
    if ( draw ) optio += "Q";
  };
  //
  if ( debug ) printf(" OK, start the fitting procedure... options: %s \n",optio.Data());
  //
  //  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);
  //
  E0 = lfit->GetParameter(0);
  a = lfit->GetParameter(1);
  b = lfit->GetParameter(2);
  errE0 = lfit->GetParError(0);
  erra = lfit->GetParError(1);
  errb = lfit->GetParError(2);
  chi2 = lfit->GetChisquare();
  ndf = lfit->GetNDF();
  Float_t tmax = 0.;
  if ( debug ) printf(" Parameters are retrieved \n");
  if ( b != 0 ) tmax = (a - 1.)/b;
  //
  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);
    Float_t imax = lfit->IntegralFast(npp,xp,wp,0.,tmax);
    //    Float_t imax = lfit->Integral(0.,tmax);
    if ( debug ) printf(" Integral till maximum (%f): %f \n",tmax,imax); 
    Int_t np = 1000;
    double *x=new double[np];
    double *w=new double[np];
    lfit->CalcGaussLegendreSamplingPoints(np,x,w,1e-12);
    Float_t i10max = lfit->IntegralFast(np,x,w,0.,10.*tmax);
    delete x;
    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); 
    //
    if ( i10max != imax ){
      asymm = imax / (i10max-imax);
    } else {
      if ( debug ) printf(" i10max == imax, asymm undefined\n");
      asymm = -2.;
    };
    if ( asymm != asymm ){
      if ( debug ) printf(" asymm is nan \n");      
      asymm = -3.;
    };
    //lfit->Integral(0.,tmax)/(lfit->Integral(0.,10.*tmax)-lfit->Integral(0.,tmax));
    if ( debug ) printf(" Asymmetry has been calculated \n");
  };
  //
  if ( asymm < 0. || ndf <= 0. || chi2 < 0. || tmax < 0. ){
    if ( debug ) printf(" Funny asymm||ndf||chi2||tmax values, fit failed \n");
    fitresult = 100;
  };
  //
  if ( draw ){
    //
    tc->cd();    
    //    gStyle->SetOptStat(11111);
    tc->SetTitle();
    th->SetTitle("");
    th->SetName("");
    th->SetMarkerStyle(20);
    // axis titles
    th->SetXTitle("Depth [X0]");
    th->SetYTitle("Energy [MIP]");
    //    th->DrawCopy("Perror");
    th->DrawCopy();
    gh->SetMarkerSize(1.);
    gh->SetMarkerStyle(21);
    //    gh->SetMarkerColor(kRed);
    gh->Draw("Psame");
    lfit->Draw("same");
    tc->Modified();
    tc->Update();
    //
    gStyle->SetLabelSize(0);
    gStyle->SetNdivisions(1,"XY");
    //
  } else {
    if ( th ) th->Delete();
  };
  //
//   gMinuit->SetPrintLevel(1);
//   gMinuit->mnamin();
//   //  gMinuit->mnfree(0);
//   gMinuit->mncler();
//   gMinuit->mnrset(1);
  //  gMinuit->Delete();
  //ROOT::Math::Minimizer::Clear();
  //  gMinuit->mnhelp("*");
  //
  //  delete lfit;
  //
};

void CaloLong::Draw(){
  //
  Process();
  Draw(-1);
};

void CaloLong::Draw(Int_t view){
  //
  Int_t minv = 0;
  Int_t maxv = 0;
  //
  if ( view == -1 ){
    maxv = -1;
  } else {
    minv = view;
    maxv = view+1;
  };
  // 
  Process();
  //
  gStyle->SetLabelSize(0.04);
  gStyle->SetNdivisions(510,"XY");
  //
  if ( maxv != -1 ){
    for (Int_t v=minv; v<maxv;v++){
      TString hid = Form("clongv%i%s",v,suf.Data());    
      TCanvas *tc  = dynamic_cast<TCanvas*>(gDirectory->FindObject(hid));
      if ( tc ){
        //       tc->Clear();
      } else {
        tc = new TCanvas(hid,hid);
      };
      //
      TString thid = Form("hlongv%i%s",v,suf.Data());   
      TH1F *th  = dynamic_cast<TH1F*>(gDirectory->FindObject(thid));
      if ( th ) th->Delete();
      //         th->Clear();
      //         th->Reset();
      //        } else {
      th = new TH1F(thid,thid,22,-0.5,21.5);
      //        };
      tc->cd();
      //
      for (Int_t st=0;st<22;st++){
        th->Fill(st,eplane[v][st]);
      };
      th->Draw();
      tc->Modified();
      tc->Update();
    };
  } else {
    //
    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%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);
    tc->cd();
    Int_t pp=0;
    for (Int_t st=0;st<22;st++){
      for (Int_t v=1; v>=0;v--){
        //
        th->Fill(pp,eplane[v][st]);
        //
        pp++;
      };
    };
    th->Draw();
    tc->Modified();
    tc->Update();
  };
  //
  gStyle->SetLabelSize(0);
  gStyle->SetNdivisions(1,"XY");
  //
};

void CaloLong::Delete(){
  Clear();
  //delete this;
};

