PAMELA%20software/PamCut/CaloAxis2.h

#ifndef CALOAXIS2_H_
#define CALOAXIS2_H_

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <TString.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TMath.h>
#include <TLine.h>
#include <TPolyMarker.h>
#include <TSelector.h>
#include <TFile.h>

#include <stdlib.h>
#include <iostream>
using namespace std;

#include <PamLevel2.h>
//#include <CaloEvent.h>

#endif // !defined(__CINT__) || defined(__MAKECINT__)
//===============================================================================
//
//
//
//
//===============================================================================
class CaloEvent {

public:
  //  calo event
  //  ------------------------------
  //  hit strips
  //  ------------------------------
  vector<int> ids; //strip 0-95
  vector<int> idp; //plane 0-21
  vector<float> q; // signal
  vector<float> zcoord;
  vector<float> xycoord;

  //  ------------------------------
  //  methods
  //  ------------------------------
  void Add(int iis, int iip, float fq, float fz, float fxy) {
    ids.push_back(iis);
    idp.push_back(iip);
    q.push_back(fq);
    zcoord.push_back(fz);
    xycoord.push_back(fxy);
  }
  ;

  int GetN() {
    return (int) (q.size());
  }
  ;

  void Reset() {
    ids.clear();
    idp.clear();
    q.clear();
    zcoord.clear();
    xycoord.clear();
  }
  ;

  void Delete() {
    Reset();
  }
  ;

  ~CaloEvent() {
    Delete();
  }
  ;

  CaloEvent() {
  }
  ;

  CaloEvent(CaloLevel1* calo, int view, Bool_t usemechanicalalignement) {
    Set(calo, view, usemechanicalalignement);
  }
  ;

  CaloEvent(CaloLevel1* calo, int view) {
    Set(calo, view);
  }
  ;

  void Set(CaloLevel1* calo, int view, Bool_t usemechanicalalignement) {

    Reset();
    //  cout << " CaloEvent::Set()"<<endl;
    CaloStrip cstrip;
    cstrip = CaloStrip(calo, usemechanicalalignement);

    for (int ih = 0; ih < calo->istrip; ih++) {
      int iv = -1;
      int ip = -1;
      int is = -1;
      calo->DecodeEstrip(ih, iv, ip, is);
      if (iv == view) {
        cstrip.Set(view, ip, is);
        float fxy = 0.;
        if (!view)
          fxy = cstrip.GetX();
        else
          fxy = cstrip.GetY();
        float fz = cstrip.GetZ();
        float fq = cstrip.GetE();
        if (fq > 0)
          Add(is, ip, fq, fz, fxy);
      }
    }

  }
  ;
  void Set(CaloLevel1* calo, int view) {
    Set(calo, view, 0);
  }
  ;

};
//===============================================================================
//
//
//
//
//===============================================================================
class CaloAxis {

public:

  CaloEvent *cevent;

  //  ------------------------------
  //  selected points along the axis
  //  ------------------------------
  //  fitted points
  vector<float> x;// independent coordinates (z)
  vector<float> y;// dependente coordinate (x/y)
  vector<float> w;// fit weight
  vector<float> q;// signal

  //  ------------------------------
  //  axis parameters
  //  ------------------------------
  Float_t par[2];
  Float_t covpar[2][2];

  //  ------------------------------
  //  general variables
  //  ------------------------------
  float cog[22];// baricenter
  float qpl[22]; // charge in each plane

  //  ------------------------------
  //  objects to draw the axis
  //  ------------------------------
  TPolyMarker *sel;
  TLine *axis;

  //  ------------------------------
  //  methods
  //  ------------------------------
  void Init() {
    for (Int_t i = 0; i < 2; i++) {
      par[i] = 0;
      for (Int_t j = 0; j < 2; j++)
        covpar[i][j] = 0;
    };
    sel = 0;
    axis = 0;
    for (Int_t i = 0; i < 22; i++) {
      cog[i] = 0;
      qpl[i] = 0;
    }
    //  if(cevent)cevent->Reset();
  }

  void Set(CaloLevel1* calo, Int_t view) {
    cevent->Set(calo, view, 0);
  }
  ;
  void Set(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement) {
    cevent->Set(calo, view, usemechanicalalignement);
  }
  ;

  CaloAxis() {
    cevent = new CaloEvent();
    Init();
  }
  ;

  CaloAxis(CaloLevel1* calo, Int_t view) {
    cevent = new CaloEvent();
    Init();
    Set(calo, view);
  }
  ;
  CaloAxis(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement) {
    cevent = new CaloEvent();
    Init();
    Set(calo, view, usemechanicalalignement);
  }
  ;

  void Reset() {
    //  cout << " CaloAxis::Reset()"<<endl;
    x.clear();
    y.clear();
    w.clear();
    q.clear();
    //  cout << sel << endl;
    if (sel)
      sel->Delete();
    //  cout << axis << endl;
    if (axis)
      axis->Delete();
    //  if(cevent)cevent->Reset();
    Init();
  }

  void Delete() {
    Reset();
    delete cevent;
  }
  ;

  ~CaloAxis() {
    Delete();
  }
  ;

  void Add(float xin, float yin) {
    x.push_back(xin);
    y.push_back(yin);
    w.push_back(1.);
    q.push_back(0.);
    /*  Int_t nin = x.size(); */
    /*  if(nin>2)Fit(); */
  }
  ;
  void Add(float xin, float yin, float qin) {
    x.push_back(xin);
    y.push_back(yin);
    q.push_back(qin);
    w.push_back(1.);
    /*  Int_t nin = x.size(); */
    /*  if(nin>2)Fit(); */
  }
  ;
  void Add(float xin, float yin, float qin, float win) {
    x.push_back(xin);
    y.push_back(yin);
    w.push_back(win);
    q.push_back(qin);
    /*  Int_t nin = x.size(); */
    /*  if(nin>2)Fit(); */
  }
  ;

  int Fit() {

    //  cout << " CaloAxis::Fit()"<<endl;
    Float_t SSS = 0;
    Float_t SSX = 0;
    Float_t SXX = 0;

    Float_t SSY = 0;
    Float_t SXY = 0;

    Int_t np = 0;
    for (Int_t i = 0; i < (int) (x.size()); i++) {
      SSS += w[i] * w[i];
      SSX += x[i] * w[i] * w[i];
      SXX += x[i] * x[i] * w[i] * w[i];
      SSY += y[i] * w[i] * w[i];
      SXY += x[i] * y[i] * w[i] * w[i];
      if (w[i])
        np++;
    }
    Float_t det = SSS * SXX - SSX * SSX;
    if (det == 0)
      return 0;
    if (np < 3)
      return 0;
    //  cout << np << " points fitted -- ok "<<endl;
    par[0] = (SSY * SXX - SXY * SSX) / det;
    par[1] = (SXY * SSS - SSY * SSX) / det;
    return 1;
  }
  ;

  void Print() {
    // useful for debug for(Int_t i=0; i<(int)(x.size()); i++)if(w[i])cout << x[i] << " - "<<y[i]<<endl;
  }

  float GetPar(Int_t i) {
    return par[i];
  }
  ;

  int GetN() {
    return (int) (x.size());
  }
  ;

  float* GetX() {
    float *xout = new float[(int) (x.size())];
    for (Int_t i = 0; i < (int) (x.size()); i++)
      xout[i] = x[i];
    return xout;
  }

  float* GetY() {
    float *yout = new float[(int) (y.size())];
    for (Int_t i = 0; i < (int) (y.size()); i++)
      yout[i] = y[i];
    return yout;
  }

  float* GetQ() {
    float *qout = new float[(int) (q.size())];
    for (int i = 0; i < (int) (q.size()); i++)
      qout[i] = q[i];
    return qout;
  }

  float GetChi2() {
    float chi2 = 0;
    int nchi = 0;
    for (Int_t i = 0; i < (int) (y.size()); i++) {
      chi2 += w[i] * w[i] * (y[i] - par[0] - par[1] * x[i]) * (y[i] - par[0] - par[1] * x[i]);
      nchi += (int) w[i];
    };

    if (nchi < 3)
      return -999; // original
    //  if(nchi<3 || chi2==0)return -999; // modified by Sergio

    chi2 = chi2 / (nchi - 2);
    return chi2;
  }

  float GetQaxis() {
    float qaxis = 0;
    for (Int_t i = 0; i < (int) (q.size()); i++)
      qaxis += q[i];
    //  for(Int_t i=0; i<(int)(q.size()); i++) cout<<q[i]<<endl;
    return qaxis;
  }

  float GetQaxis(float toll) {
    float qaxis = 0;
    for (int ih = 0; ih < cevent->GetN(); ih++) {
      float x = cevent->xycoord[ih];
      float z = cevent->zcoord[ih];
      float q = cevent->q[ih];
      //int ip = cevent->idp[ih];
      float d = fabs(x - par[0] - par[1] * z) / sqrt(par[1] * par[1] + 1);
      if (d < toll)
        qaxis += q;
    }
    return qaxis;
  }

  float GetCOG(Int_t ip) {
    if (ip >= 0 && ip < 22)
      return cog[ip];
    else
      return 0;
  }

  float GetQ(Int_t ip) {
    if (ip >= 0 && ip < 22)
      return qpl[ip];
    else
      return 0;
  }

  float GetQstrip() {

    float qaxis = 0;
    for (Int_t i = 0; i < (int) (q.size()); i++)
      qaxis += q[i];
    if (q.size() == 0)
      return 0.;
    return qaxis / q.size();
  }

  float GetYfit(float x) {
    return par[0] + x * par[1];
  }

  float GetXmin() {
    float zmin = 9999;
    for (Int_t i = 0; i < (int) (x.size()); i++)
      if (x[i] < zmin)
        zmin = x[i];
    return zmin;
  }
  float GetXmax() {
    float zmax = -9999;
    for (Int_t i = 0; i < (int) (x.size()); i++)
      if (x[i] > zmax)
        zmax = x[i];
    return zmax;

  }

  void Draw() {

    //  cout << " CaloAxis::Draw()"<<endl;
    if (GetN() < 3)
      return;

    sel = new TPolyMarker(GetN(), GetY(), GetX());
    sel->SetMarkerSize(0.3);
    sel->SetMarkerColor(5);

    axis = new TLine(par[0] + GetXmin() * par[1], GetXmin(), par[0] + GetXmax() * par[1], GetXmax());
    axis->SetLineWidth(1);
    axis->SetLineColor(5);

    //  cout << sel << endl;
    //  cout << axis << endl;

    sel->Draw();
    axis->Draw();
  }
  ;

  //-----------------------------------------------------------------
  // fit the axis (optimized for non-interacting patterns)
  //-----------------------------------------------------------------
  int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll, Int_t nPlanes = 22, Bool_t usemechanicalalignement = 0) {

    //  cout << "CaloAxis::FitAxis(...)"<<endl;

    Set(calo, view, usemechanicalalignement);

    // ------------------------------
    // first : all hits included, w=1
    // ------------------------------
    for (int ih = 0; ih < cevent->GetN(); ih++) {
      float x = cevent->xycoord[ih];
      float z = cevent->zcoord[ih];
      float q = cevent->q[ih];
      Add(z, x, q);
    }

    if (GetN() < 3)
      return 0;
    Fit();
    // useful for debug cout << "(1)"<<endl;
    //  cout << " n. hits :"<<GetN()<<endl;
    //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
    //  cout << " Chi2    :"<<GetChi2()<<endl;
    //  cout << " Q axis  :"<<GetQaxis()<<endl;
    //  cout << " Q strip :"<<GetQstrip()<<endl;

    float P0;
    float P1;
    int dof;

    // ---------------------------
    // second : iteration
    // ---------------------------
    int niter = 0; // number of iterations
    int pok[nPlanes];

    // prova
    // pesi lineari
    float W22 = 1.;
    float W1 = 10.;
    float b = (W22 - W1) / ((float) nPlanes - 1.);// (w22 - w1)/(22-1)
    float a = W1 - b;

    for (int i = 0; i < nPlanes; i++)
      pok[i] = 0;
    do {
      niter++;
      float ttoll = toll; //tolerance (cm)
      if (niter == 1)
        ttoll = 10. * toll;
      dof = GetN(); //previous fit
      P0 = par[0];
      P1 = par[1];
      Reset(); //fit reset
      for (int ih = 0; ih < cevent->GetN(); ih++) {//loop over selected hits
        if (cevent->idp[ih] < nPlanes) {
          float x = cevent->xycoord[ih];
          float z = cevent->zcoord[ih];
          float q = cevent->q[ih];
          int ip = cevent->idp[ih];
          float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);//distance to axis
          //            cout<<d<<endl;
          if (d < ttoll && (niter == 1 || (niter > 1 && pok[ip] == 1))) {
            //          if( d < ttoll ){
            float W = a + b * (ip + 1);
            //              cout << ip << " "<<W<<endl;
            Add(z, x, q, W);
            pok[ip] = 1;
          }
        }
      }
      // break the track if more than 3 planes are missing
      int nmissing = 0;
      for (int ip = 0; ip < nPlanes; ip++) {
        if (pok[ip] == 0) {
          nmissing++;
        }
        else {
          nmissing = 0;
        }
        if (nmissing == 3) {
          for (int ipp = ip + 1; ipp < nPlanes; ipp++)
            pok[ipp] = 0;
          break;
        }
      }

      if (niter == 100)
        break;
      if (GetN() < 3)
        break;
      Fit();
    } while (niter == 1 || GetN() != dof);
    //  cout << " >>> "<<GetN()<<endl;
    if (GetN() < 3)
      return 0;
    Fit();
    //  cout << "(2)"<<endl;
    //  cout << " n. hits :"<<GetN()<<endl;
    //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
    //  cout << " Chi2    :"<<GetChi2()<<endl;
    //  cout << " Q axis  :"<<GetQaxis()<<endl;
    //  cout << " Q strip :"<<GetQstrip()<<endl;

    // ---------------------------------------------
    // third : selecting only closest strip-clusters
    // and fit baricenters
    // ---------------------------------------------
    P0 = par[0];
    P1 = par[1];
    Reset();

    float dmin[nPlanes];
    int closest[nPlanes];
    for (int ip = 0; ip < nPlanes; ip++) {
      dmin[ip] = 999;
      closest[ip] = -1;
    }
    for (int ih = 0; ih < cevent->GetN(); ih++) {
      if (cevent->idp[ih] < nPlanes) {
        float x = cevent->xycoord[ih];
        float z = cevent->zcoord[ih];
        //float q = cevent->q[ih];
        int ip = cevent->idp[ih];
        //int is = cevent->ids[ih];
        float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);
        if (d < toll && d < dmin[ip] && pok[ip] == 1) {
          //            Add(z,x,q,q);
          closest[ip] = ih;
          //            cog[ip] += x*q;
          //            qpl[ip] += q;
        }
      }
    }
    for (Int_t ip = 0; ip < nPlanes; ip++) {
      if (closest[ip] != -1) {
        float x = cevent->xycoord[closest[ip]];
        float z = cevent->zcoord[closest[ip]];
        float q = cevent->q[closest[ip]];
        //int is = cevent->ids[closest[ip]];
        Add(z, x, q, q);
        cog[ip] += x * q;
        qpl[ip] += q;
        //              cout << ip << " -o- "<<is<<endl;
      }
    }
    // -----------------
    // add +/- one strip
    // -----------------
    for (int ih = 0; ih < cevent->GetN(); ih++) {
      if (cevent->idp[ih] < nPlanes) {
      float x = cevent->xycoord[ih];
      float z = cevent->zcoord[ih];
      float q = cevent->q[ih];
      int ip = cevent->idp[ih];
      int is = cevent->ids[ih];
      if (closest[ip] != -1) {
        int isc = cevent->ids[closest[ip]];
        if (is == isc + 1 || is == isc - 1) {
          Add(z, x, q, q);
          cog[ip] += x * q;
          qpl[ip] += q;
          //                cout << ip << " -+- "<<is<<endl;
        }
      }}
    }
    Fit();

    if (GetN() < 3)
      return 0;
    for (int ip = 0; ip < nPlanes; ip++) {
      if (qpl[ip] != 0) {
        cog[ip] = cog[ip] / qpl[ip];
        //          Add(z,cog[ip],cog[ip],0.7);
      }
    }

    //  cout << "(3)"<<endl;
    //          cout << " n. hits :"<<GetN()<<endl;
    //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
    //          cout << " Chi2    :"<<GetChi2()<<endl;
    //  cout << " Q axis  :"<<GetQaxis()<<endl;
    //          cout << " Q strip :"<<GetQstrip()<<endl;
    //  cout << " ---------------------------"<<endl;
    return 1;
  }
  ;

  /*int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll, Int_t nPlanes = 22) {
    return FitAxis(calo, view, toll, 0, nPlanes);
  }
  ;*/

  //-----------------------------------------------------------------
  // fit the axis (optimized for interacting patterns?? ...forse)
  //-----------------------------------------------------------------
  int FitShower(CaloLevel1* calo, Int_t view, Float_t toll, Bool_t usemechanicalalignement) {

    //  cout << "CaloAxis::FitAxis(...)"<<endl;

    Set(calo, view, usemechanicalalignement);

    // ------------------------------
    // first :
    // ------------------------------
    for (int ih = 0; ih < cevent->GetN(); ih++) {
      float x = cevent->xycoord[ih];
      //float z = cevent->zcoord[ih];
      float q = cevent->q[ih];
      int ip = cevent->idp[ih];
      cog[ip] += x * q;
      qpl[ip] += q;
      //            Add(z,x,q);
    }
    for (int ip = 0; ip < 22; ip++) {
      if (qpl[ip] != 0)
        cog[ip] = cog[ip] / qpl[ip];
      //            cout << ip << " > "<<qpl[ip]<<" "<<cog[ip]<<endl;
    }
    // ----------------
    // plane of maximum
    // ----------------
    float qplmax = 1;
    int ipmax = -1;
    float dmin[22];
    int closest[22];
    for (int ip = 0; ip < 22; ip++)
      if (qpl[ip] > qplmax) {
        ipmax = ip;
        qplmax = qpl[ip];
        dmin[ip] = 1000.;//init
        closest[ip] = -1;//init
      }
    for (int ih = 0; ih < cevent->GetN(); ih++) {
      float x = cevent->xycoord[ih];
      float z = cevent->zcoord[ih];
      float q = cevent->q[ih];
      int ip = cevent->idp[ih];
      //            if( (ipmax>3 && ip<=ipmax) || ipmax<=3 )Add(z,x,q);
      if (ip <= ipmax || ip <= 3) {
        Add(z, x, q);
      }
    }
    if (GetN() < 3)
      return 0;
    Fit();
    //  cout << " n. hits :"<<GetN()<<endl;
    //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
    //  cout << " Chi2    :"<<GetChi2()<<endl;
    //  cout << " Q axis  :"<<GetQaxis()<<endl;
    //  cout << " Q strip :"<<GetQstrip()<<endl;

    float P0;
    float P1;
    int dof;

    // ---------------------------
    // second : iteration
    // ---------------------------
    int niter = 0; // number of iterations
    int pok[22];
    // prova
    // pesi lineari
    float W22 = 1.;
    float W1 = 10.;
    float b = (W22 - W1) / (22. - 1.);// (w22 - w1)/(22-1)
    float a = W1 - b;

    for (int i = 0; i < 22; i++) {
      pok[i] = 0;
    }
    do {
      niter++;
      float ttoll = toll; //tolerance (cm)
      if (niter == 1)
        ttoll = 10. * toll;
      dof = GetN(); //previous fit
      P0 = par[0];
      P1 = par[1];
      Reset(); //fit reset
      for (int i = 0; i < 22; i++) {
        cog[i] = 0;
        qpl[i] = 0;
      }
      for (int ih = 0; ih < cevent->GetN(); ih++) {
        float x = cevent->xycoord[ih];
        float z = cevent->zcoord[ih];
        float q = cevent->q[ih];
        int ip = cevent->idp[ih];
        float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);
        //              cout<<d<<endl;
        if ((niter == 1 || (niter > 1 && pok[ip] == 1)) && (ip <= ipmax || ip <= 3) &&
        //                  ((ipmax>3 && ip<=ipmax) || ipmax<=3) &&
            true) {

          if (niter == 1 && d < dmin[ip]) {
            dmin[ip] = d;
            closest[ip] = ih;
          }

          float W = a + b * (ip + 1);
          //                cout << ip << " "<<W<<endl;
          if (d < ttoll || (niter == 2 && ih == closest[ip])) {
            /*                          if(ip==ipmax && ipmax>1)Add(z,x,q,W*q);  */
            /*                          else Add(z,x,q,W);  */
            cog[ip] += x * q;
            qpl[ip] += q;
            Add(z, x, q, W);
            pok[ip] = 1;
            cout << ip << " " << qpl[ip] << endl;
          }
        }
      }
      // break the track if more than 3 planes are missing
      int nmissing = 0;
      for (int ip = 0; ip < 22; ip++) {
        if (pok[ip] == 0) {
          nmissing++;
        }
        else {
          nmissing = 0;
        }
        if (nmissing == 6) {
          for (int ipp = ip + 1; ipp < 22; ipp++)
            pok[ipp] = 0;
          break;
        }
      }

      for (int ip = 0; ip < 22; ip++)
        cout << qpl[ip] << " ";
      cout << endl;

      for (int ip = 0; ip < 22; ip++)
        if (qpl[ip] != 0)
          cog[ip] = cog[ip] / qpl[ip];

      if (niter == 100)
        break;
      if (GetN() < 3)
        break;

      Fit();

    } while (niter == 1 || GetN() != dof);
    //  cout << " >>> "<<GetN()<<endl;
    if (GetN() < 3)
      return 0;
    Fit();
    /* //       cout << " n. hits :"<<GetN()<<endl; */
    /* //       cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl; */
    /* //       cout << " Chi2    :"<<GetChi2()<<endl; */
    /* //       cout << " Q axis  :"<<GetQaxis()<<endl; */
    /* //       cout << " Q strip :"<<GetQstrip()<<endl; */

    /*  // --------------------------------------------- */
    /*  // third : selecting only closest strip-clusters */
    /*  // and fit baricenters */
    /*  // --------------------------------------------- */
    /*  P0  = par[0]; */
    /*  P1  = par[1]; */
    /*  Reset(); */

    /*  float dmin[22]; */
    /*  int closest[22]; */
    /*  for(int ip=0; ip<22; ip++){ */
    /*      dmin[ip]=999; */
    /*      closest[ip]=-1; */
    /*  } */
    /*  for(int ih=0;ih<cevent->GetN(); ih++){ */
    /*      float x = cevent->xycoord[ih]; */
    /*      float z = cevent->zcoord[ih]; */
    /*      float q = cevent->q[ih]; */
    /*      int ip = cevent->idp[ih]; */
    /*      int is = cevent->ids[ih]; */
    /*      float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1); */
    /*      if( d < toll && d<dmin[ip] && pok[ip]==1 ){ */
    /* //               Add(z,x,q,q); */
    /*          closest[ip]=ih; */
    /* //               cog[ip] += x*q; */
    /* //               qpl[ip] += q; */
    /*      } */
    /*  } */
    /*  for(Int_t ip=0; ip<22; ip++){ */
    /*      if(closest[ip]!=-1){ */
    /*          float x = cevent->xycoord[closest[ip]]; */
    /*          float z = cevent->zcoord[closest[ip]]; */
    /*          float q = cevent->q[closest[ip]]; */
    /*          int is = cevent->ids[closest[ip]]; */
    /*          Add(z,x,q,q); */
    /*          cog[ip] += x*q; */
    /*          qpl[ip] += q; */
    /* //               cout << ip << " -o- "<<is<<endl; */
    /*      } */
    /*  } */
    /*  // add +/- one strip */
    /*  for(int ih=0;ih<cevent->GetN(); ih++){ */
    /*      float x = cevent->xycoord[ih]; */
    /*      float z = cevent->zcoord[ih]; */
    /*      float q = cevent->q[ih]; */
    /*      int ip = cevent->idp[ih]; */
    /*      int is = cevent->ids[ih]; */
    /*      if( closest[ip]!=-1 ){ */
    /*          int isc = cevent->ids[closest[ip]]; */
    /*          if( is == isc+1 || is == isc-1 ){ */
    /*              Add(z,x,q,q); */
    /*              cog[ip] += x*q; */
    /*              qpl[ip] += q; */
    /* //                   cout << ip << " -+- "<<is<<endl; */
    /*          } */
    /*      } */
    /*  } */
    /*  Fit(); */

    /*  if(GetN()<3)return 0; */
    /*  for(int ip=0; ip<22; ip++){ */
    /*      if(qpl[ip]!=0){ */
    /*          cog[ip]=cog[ip]/qpl[ip]; */
    /* //           Add(z,cog[ip],cog[ip],0.7); */
    /*      } */
    /*  } */

    //  cout << " n. hits :"<<GetN()<<endl;
    //  cout << " P0 P1   :"<<par[0]<< " "<<par[1]<<endl;
    //  cout << " Chi2    :"<<GetChi2()<<endl;
    //  cout << " Q axis  :"<<GetQaxis()<<endl;
    //  cout << " Q strip :"<<GetQstrip()<<endl;
    //  cout << " ---------------------------"<<endl;
    return 1;
  }
  ;

  int FitShower(CaloLevel1* calo, Int_t view, Float_t toll) {
    return FitAxis(calo, view, toll, 0);
  }
  ;

};

//===============================================================================
//
//
//
//
//===============================================================================

#endif // CALOAXIS2_H_