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, Bool_t usemechanicalalignement) {

    //  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[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 ih = 0; ih < cevent->GetN(); ih++) {//loop over selected hits
        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 < 22; ip++) {
        if (pok[ip] == 0) {
          nmissing++;
        }
        else {
          nmissing = 0;
        }
        if (nmissing == 3) {
          for (int ipp = ip + 1; ipp < 22; 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[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 << "(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) {
    return FitAxis(calo, view, toll, 0);
  }
  ;

  //-----------------------------------------------------------------
  // 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_
1	pam-fi	1.1	#ifndef CALOAXIS2_H_
2			#define CALOAXIS2_H_
3
4			#if !defined(__CINT__) \|\| defined(__MAKECINT__)
5
6			#include <TString.h>
7			#include <TH1F.h>
8			#include <TH2F.h>
9			#include <TMath.h>
10			#include <TLine.h>
11			#include <TPolyMarker.h>
12			#include <TSelector.h>
13			#include <TFile.h>
14
15			#include <stdlib.h>
16			#include <iostream>
17			using namespace std;
18
19			#include <PamLevel2.h>
20			//#include <CaloEvent.h>
21
22			#endif // !defined(__CINT__) \|\| defined(__MAKECINT__)
23			//===============================================================================
24			//
25			//
26			//
27			//
28			//===============================================================================
29	pam-fi	1.2	class CaloEvent {
30	pam-fi	1.1
31			public:
32	pam-fi	1.2	// calo event
33			// ------------------------------
34			// hit strips
35			// ------------------------------
36			vector<int> ids; //strip 0-95
37			vector<int> idp; //plane 0-21
38			vector<float> q; // signal
39			vector<float> zcoord;
40			vector<float> xycoord;
41
42			// ------------------------------
43			// methods
44			// ------------------------------
45			void Add(int iis, int iip, float fq, float fz, float fxy) {
46			ids.push_back(iis);
47			idp.push_back(iip);
48			q.push_back(fq);
49			zcoord.push_back(fz);
50			xycoord.push_back(fxy);
51			}
52			;
53
54			int GetN() {
55			return (int) (q.size());
56			}
57			;
58
59			void Reset() {
60			ids.clear();
61			idp.clear();
62			q.clear();
63			zcoord.clear();
64			xycoord.clear();
65			}
66			;
67
68			void Delete() {
69			Reset();
70			}
71			;
72
73			~CaloEvent() {
74			Delete();
75			}
76			;
77
78			CaloEvent() {
79			}
80			;
81
82			CaloEvent(CaloLevel1* calo, int view, Bool_t usemechanicalalignement) {
83			Set(calo, view, usemechanicalalignement);
84			}
85			;
86
87			CaloEvent(CaloLevel1* calo, int view) {
88			Set(calo, view);
89			}
90			;
91
92			void Set(CaloLevel1* calo, int view, Bool_t usemechanicalalignement) {
93
94			Reset();
95			// cout << " CaloEvent::Set()"<<endl;
96			CaloStrip cstrip;
97			cstrip = CaloStrip(calo, usemechanicalalignement);
98
99			for (int ih = 0; ih < calo->istrip; ih++) {
100			int iv = -1;
101			int ip = -1;
102			int is = -1;
103			calo->DecodeEstrip(ih, iv, ip, is);
104			if (iv == view) {
105			cstrip.Set(view, ip, is);
106			float fxy = 0.;
107			if (!view)
108			fxy = cstrip.GetX();
109			else
110			fxy = cstrip.GetY();
111			float fz = cstrip.GetZ();
112			float fq = cstrip.GetE();
113			if (fq > 0)
114			Add(is, ip, fq, fz, fxy);
115			}
116			}
117
118			}
119			;
120			void Set(CaloLevel1* calo, int view) {
121			Set(calo, view, 0);
122			}
123			;
124	pam-fi	1.1
125			};
126			//===============================================================================
127			//
128			//
129			//
130			//
131			//===============================================================================
132			class CaloAxis {
133
134			public:
135
136	pam-fi	1.2	CaloEvent *cevent;
137	pam-fi	1.1
138	pam-fi	1.2	// ------------------------------
139			// selected points along the axis
140			// ------------------------------
141			// fitted points
142			vector<float> x;// independent coordinates (z)
143			vector<float> y;// dependente coordinate (x/y)
144			vector<float> w;// fit weight
145			vector<float> q;// signal
146
147			// ------------------------------
148			// axis parameters
149			// ------------------------------
150			Float_t par[2];
151			Float_t covpar[2][2];
152
153			// ------------------------------
154			// general variables
155			// ------------------------------
156			float cog[22];// baricenter
157			float qpl[22]; // charge in each plane
158
159			// ------------------------------
160			// objects to draw the axis
161			// ------------------------------
162			TPolyMarker *sel;
163			TLine *axis;
164
165			// ------------------------------
166			// methods
167			// ------------------------------
168			void Init() {
169			for (Int_t i = 0; i < 2; i++) {
170			par[i] = 0;
171			for (Int_t j = 0; j < 2; j++)
172			covpar[i][j] = 0;
173			};
174			sel = 0;
175			axis = 0;
176			for (Int_t i = 0; i < 22; i++) {
177			cog[i] = 0;
178			qpl[i] = 0;
179			}
180			// if(cevent)cevent->Reset();
181			}
182
183			void Set(CaloLevel1* calo, Int_t view) {
184			cevent->Set(calo, view, 0);
185			}
186			;
187			void Set(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement) {
188			cevent->Set(calo, view, usemechanicalalignement);
189			}
190			;
191
192			CaloAxis() {
193			cevent = new CaloEvent();
194			Init();
195			}
196			;
197
198			CaloAxis(CaloLevel1* calo, Int_t view) {
199			cevent = new CaloEvent();
200			Init();
201			Set(calo, view);
202			}
203			;
204			CaloAxis(CaloLevel1* calo, Int_t view, Bool_t usemechanicalalignement) {
205			cevent = new CaloEvent();
206			Init();
207			Set(calo, view, usemechanicalalignement);
208			}
209			;
210
211			void Reset() {
212			// cout << " CaloAxis::Reset()"<<endl;
213			x.clear();
214			y.clear();
215			w.clear();
216			q.clear();
217			// cout << sel << endl;
218			if (sel)
219			sel->Delete();
220			// cout << axis << endl;
221			if (axis)
222			axis->Delete();
223			// if(cevent)cevent->Reset();
224			Init();
225			}
226
227			void Delete() {
228			Reset();
229			delete cevent;
230			}
231			;
232
233			~CaloAxis() {
234			Delete();
235			}
236			;
237
238			void Add(float xin, float yin) {
239			x.push_back(xin);
240			y.push_back(yin);
241			w.push_back(1.);
242			q.push_back(0.);
243			/* Int_t nin = x.size(); */
244			/* if(nin>2)Fit(); */
245			}
246			;
247			void Add(float xin, float yin, float qin) {
248			x.push_back(xin);
249			y.push_back(yin);
250			q.push_back(qin);
251			w.push_back(1.);
252			/* Int_t nin = x.size(); */
253			/* if(nin>2)Fit(); */
254			}
255			;
256			void Add(float xin, float yin, float qin, float win) {
257			x.push_back(xin);
258			y.push_back(yin);
259			w.push_back(win);
260			q.push_back(qin);
261			/* Int_t nin = x.size(); */
262			/* if(nin>2)Fit(); */
263			}
264			;
265
266			int Fit() {
267
268			// cout << " CaloAxis::Fit()"<<endl;
269			Float_t SSS = 0;
270			Float_t SSX = 0;
271			Float_t SXX = 0;
272
273			Float_t SSY = 0;
274			Float_t SXY = 0;
275
276			Int_t np = 0;
277			for (Int_t i = 0; i < (int) (x.size()); i++) {
278			SSS += w[i] * w[i];
279			SSX += x[i] * w[i] * w[i];
280			SXX += x[i] * x[i] * w[i] * w[i];
281			SSY += y[i] * w[i] * w[i];
282			SXY += x[i] * y[i] * w[i] * w[i];
283			if (w[i])
284			np++;
285			}
286			Float_t det = SSS * SXX - SSX * SSX;
287			if (det == 0)
288			return 0;
289			if (np < 3)
290			return 0;
291			// cout << np << " points fitted -- ok "<<endl;
292			par[0] = (SSY * SXX - SXY * SSX) / det;
293			par[1] = (SXY * SSS - SSY * SSX) / det;
294			return 1;
295			}
296			;
297
298			void Print() {
299			// useful for debug for(Int_t i=0; i<(int)(x.size()); i++)if(w[i])cout << x[i] << " - "<<y[i]<<endl;
300			}
301
302			float GetPar(Int_t i) {
303			return par[i];
304			}
305			;
306
307			int GetN() {
308			return (int) (x.size());
309			}
310			;
311
312			float* GetX() {
313			float *xout = new float[(int) (x.size())];
314			for (Int_t i = 0; i < (int) (x.size()); i++)
315			xout[i] = x[i];
316			return xout;
317			}
318
319			float* GetY() {
320			float *yout = new float[(int) (y.size())];
321			for (Int_t i = 0; i < (int) (y.size()); i++)
322			yout[i] = y[i];
323			return yout;
324			}
325
326			float* GetQ() {
327			float *qout = new float[(int) (q.size())];
328			for (int i = 0; i < (int) (q.size()); i++)
329			qout[i] = q[i];
330			return qout;
331			}
332
333			float GetChi2() {
334			float chi2 = 0;
335			int nchi = 0;
336			for (Int_t i = 0; i < (int) (y.size()); i++) {
337			chi2 += w[i] * w[i] * (y[i] - par[0] - par[1] * x[i]) * (y[i] - par[0] - par[1] * x[i]);
338			nchi += (int) w[i];
339			};
340
341			if (nchi < 3)
342			return -999; // original
343			// if(nchi<3 \|\| chi2==0)return -999; // modified by Sergio
344
345			chi2 = chi2 / (nchi - 2);
346			return chi2;
347			}
348
349			float GetQaxis() {
350			float qaxis = 0;
351			for (Int_t i = 0; i < (int) (q.size()); i++)
352			qaxis += q[i];
353			// for(Int_t i=0; i<(int)(q.size()); i++) cout<<q[i]<<endl;
354			return qaxis;
355			}
356
357			float GetQaxis(float toll) {
358			float qaxis = 0;
359			for (int ih = 0; ih < cevent->GetN(); ih++) {
360			float x = cevent->xycoord[ih];
361			float z = cevent->zcoord[ih];
362			float q = cevent->q[ih];
363			//int ip = cevent->idp[ih];
364			float d = fabs(x - par[0] - par[1] * z) / sqrt(par[1] * par[1] + 1);
365			if (d < toll)
366			qaxis += q;
367			}
368			return qaxis;
369			}
370
371			float GetCOG(Int_t ip) {
372			if (ip >= 0 && ip < 22)
373			return cog[ip];
374			else
375			return 0;
376			}
377
378			float GetQ(Int_t ip) {
379			if (ip >= 0 && ip < 22)
380			return qpl[ip];
381			else
382			return 0;
383			}
384
385			float GetQstrip() {
386
387			float qaxis = 0;
388			for (Int_t i = 0; i < (int) (q.size()); i++)
389			qaxis += q[i];
390			if (q.size() == 0)
391			return 0.;
392			return qaxis / q.size();
393			}
394
395			float GetYfit(float x) {
396			return par[0] + x * par[1];
397			}
398
399			float GetXmin() {
400			float zmin = 9999;
401			for (Int_t i = 0; i < (int) (x.size()); i++)
402			if (x[i] < zmin)
403			zmin = x[i];
404			return zmin;
405			}
406			float GetXmax() {
407			float zmax = -9999;
408			for (Int_t i = 0; i < (int) (x.size()); i++)
409			if (x[i] > zmax)
410			zmax = x[i];
411			return zmax;
412
413			}
414
415			void Draw() {
416
417			// cout << " CaloAxis::Draw()"<<endl;
418			if (GetN() < 3)
419			return;
420
421			sel = new TPolyMarker(GetN(), GetY(), GetX());
422			sel->SetMarkerSize(0.3);
423			sel->SetMarkerColor(5);
424
425			axis = new TLine(par[0] + GetXmin() * par[1], GetXmin(), par[0] + GetXmax() * par[1], GetXmax());
426			axis->SetLineWidth(1);
427			axis->SetLineColor(5);
428
429			// cout << sel << endl;
430			// cout << axis << endl;
431
432			sel->Draw();
433			axis->Draw();
434			}
435			;
436
437			//-----------------------------------------------------------------
438			// fit the axis (optimized for non-interacting patterns)
439			//-----------------------------------------------------------------
440			int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll, Bool_t usemechanicalalignement) {
441
442			// cout << "CaloAxis::FitAxis(...)"<<endl;
443
444			Set(calo, view, usemechanicalalignement);
445
446			// ------------------------------
447			// first : all hits included, w=1
448			// ------------------------------
449			for (int ih = 0; ih < cevent->GetN(); ih++) {
450			float x = cevent->xycoord[ih];
451			float z = cevent->zcoord[ih];
452			float q = cevent->q[ih];
453			Add(z, x, q);
454			}
455
456			if (GetN() < 3)
457			return 0;
458			Fit();
459			// useful for debug cout << "(1)"<<endl;
460			// cout << " n. hits :"<<GetN()<<endl;
461			// cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
462			// cout << " Chi2 :"<<GetChi2()<<endl;
463			// cout << " Q axis :"<<GetQaxis()<<endl;
464			// cout << " Q strip :"<<GetQstrip()<<endl;
465
466			float P0;
467			float P1;
468			int dof;
469
470			// ---------------------------
471			// second : iteration
472			// ---------------------------
473			int niter = 0; // number of iterations
474			int pok[22];
475
476			// prova
477			// pesi lineari
478			float W22 = 1.;
479			float W1 = 10.;
480			float b = (W22 - W1) / (22. - 1.);// (w22 - w1)/(22-1)
481			float a = W1 - b;
482
483			for (int i = 0; i < 22; i++)
484			pok[i] = 0;
485			do {
486			niter++;
487			float ttoll = toll; //tolerance (cm)
488			if (niter == 1)
489			ttoll = 10. * toll;
490			dof = GetN(); //previous fit
491			P0 = par[0];
492			P1 = par[1];
493			Reset(); //fit reset
494			for (int ih = 0; ih < cevent->GetN(); ih++) {//loop over selected hits
495			float x = cevent->xycoord[ih];
496			float z = cevent->zcoord[ih];
497			float q = cevent->q[ih];
498			int ip = cevent->idp[ih];
499			float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);//distance to axis
500			// cout<<d<<endl;
501			if (d < ttoll && (niter == 1 \|\| (niter > 1 && pok[ip] == 1))) {
502			// if( d < ttoll ){
503			float W = a + b * (ip + 1);
504			// cout << ip << " "<<W<<endl;
505			Add(z, x, q, W);
506			pok[ip] = 1;
507			}
508			}
509			// break the track if more than 3 planes are missing
510			int nmissing = 0;
511			for (int ip = 0; ip < 22; ip++) {
512			if (pok[ip] == 0) {
513			nmissing++;
514			}
515			else {
516			nmissing = 0;
517			}
518			if (nmissing == 3) {
519			for (int ipp = ip + 1; ipp < 22; ipp++)
520			pok[ipp] = 0;
521			break;
522			}
523			}
524
525			if (niter == 100)
526			break;
527			if (GetN() < 3)
528			break;
529			Fit();
530			} while (niter == 1 \|\| GetN() != dof);
531			// cout << " >>> "<<GetN()<<endl;
532			if (GetN() < 3)
533			return 0;
534			Fit();
535			// cout << "(2)"<<endl;
536			// cout << " n. hits :"<<GetN()<<endl;
537			// cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
538			// cout << " Chi2 :"<<GetChi2()<<endl;
539			// cout << " Q axis :"<<GetQaxis()<<endl;
540			// cout << " Q strip :"<<GetQstrip()<<endl;
541
542			// ---------------------------------------------
543			// third : selecting only closest strip-clusters
544			// and fit baricenters
545			// ---------------------------------------------
546			P0 = par[0];
547			P1 = par[1];
548			Reset();
549
550			float dmin[22];
551			int closest[22];
552			for (int ip = 0; ip < 22; ip++) {
553			dmin[ip] = 999;
554			closest[ip] = -1;
555			}
556			for (int ih = 0; ih < cevent->GetN(); ih++) {
557			float x = cevent->xycoord[ih];
558			float z = cevent->zcoord[ih];
559			//float q = cevent->q[ih];
560			int ip = cevent->idp[ih];
561			//int is = cevent->ids[ih];
562			float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);
563			if (d < toll && d < dmin[ip] && pok[ip] == 1) {
564			// Add(z,x,q,q);
565			closest[ip] = ih;
566			// cog[ip] += x*q;
567			// qpl[ip] += q;
568			}
569			}
570			for (Int_t ip = 0; ip < 22; ip++) {
571			if (closest[ip] != -1) {
572			float x = cevent->xycoord[closest[ip]];
573			float z = cevent->zcoord[closest[ip]];
574			float q = cevent->q[closest[ip]];
575			//int is = cevent->ids[closest[ip]];
576			Add(z, x, q, q);
577			cog[ip] += x * q;
578			qpl[ip] += q;
579			// cout << ip << " -o- "<<is<<endl;
580			}
581			}
582			// -----------------
583			// add +/- one strip
584			// -----------------
585			for (int ih = 0; ih < cevent->GetN(); ih++) {
586			float x = cevent->xycoord[ih];
587			float z = cevent->zcoord[ih];
588			float q = cevent->q[ih];
589			int ip = cevent->idp[ih];
590			int is = cevent->ids[ih];
591			if (closest[ip] != -1) {
592			int isc = cevent->ids[closest[ip]];
593			if (is == isc + 1 \|\| is == isc - 1) {
594			Add(z, x, q, q);
595			cog[ip] += x * q;
596			qpl[ip] += q;
597			// cout << ip << " -+- "<<is<<endl;
598			}
599			}
600			}
601			Fit();
602
603			if (GetN() < 3)
604			return 0;
605			for (int ip = 0; ip < 22; ip++) {
606			if (qpl[ip] != 0) {
607			cog[ip] = cog[ip] / qpl[ip];
608			// Add(z,cog[ip],cog[ip],0.7);
609			}
610			}
611
612			// cout << "(3)"<<endl;
613			// cout << " n. hits :"<<GetN()<<endl;
614			// cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
615			// cout << " Chi2 :"<<GetChi2()<<endl;
616			// cout << " Q axis :"<<GetQaxis()<<endl;
617			// cout << " Q strip :"<<GetQstrip()<<endl;
618			// cout << " ---------------------------"<<endl;
619			return 1;
620			}
621			;
622
623			int FitAxis(CaloLevel1* calo, Int_t view, Float_t toll) {
624			return FitAxis(calo, view, toll, 0);
625			}
626			;
627
628			//-----------------------------------------------------------------
629			// fit the axis (optimized for interacting patterns?? ...forse)
630			//-----------------------------------------------------------------
631			int FitShower(CaloLevel1* calo, Int_t view, Float_t toll, Bool_t usemechanicalalignement) {
632
633			// cout << "CaloAxis::FitAxis(...)"<<endl;
634
635			Set(calo, view, usemechanicalalignement);
636
637			// ------------------------------
638			// first :
639			// ------------------------------
640			for (int ih = 0; ih < cevent->GetN(); ih++) {
641			float x = cevent->xycoord[ih];
642			//float z = cevent->zcoord[ih];
643			float q = cevent->q[ih];
644			int ip = cevent->idp[ih];
645			cog[ip] += x * q;
646			qpl[ip] += q;
647			// Add(z,x,q);
648			}
649			for (int ip = 0; ip < 22; ip++) {
650			if (qpl[ip] != 0)
651			cog[ip] = cog[ip] / qpl[ip];
652			// cout << ip << " > "<<qpl[ip]<<" "<<cog[ip]<<endl;
653			}
654			// ----------------
655			// plane of maximum
656			// ----------------
657			float qplmax = 1;
658			int ipmax = -1;
659			float dmin[22];
660			int closest[22];
661			for (int ip = 0; ip < 22; ip++)
662			if (qpl[ip] > qplmax) {
663			ipmax = ip;
664			qplmax = qpl[ip];
665			dmin[ip] = 1000.;//init
666			closest[ip] = -1;//init
667			}
668			for (int ih = 0; ih < cevent->GetN(); ih++) {
669			float x = cevent->xycoord[ih];
670			float z = cevent->zcoord[ih];
671			float q = cevent->q[ih];
672			int ip = cevent->idp[ih];
673			// if( (ipmax>3 && ip<=ipmax) \|\| ipmax<=3 )Add(z,x,q);
674			if (ip <= ipmax \|\| ip <= 3) {
675			Add(z, x, q);
676			}
677			}
678			if (GetN() < 3)
679			return 0;
680			Fit();
681			// cout << " n. hits :"<<GetN()<<endl;
682			// cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
683			// cout << " Chi2 :"<<GetChi2()<<endl;
684			// cout << " Q axis :"<<GetQaxis()<<endl;
685			// cout << " Q strip :"<<GetQstrip()<<endl;
686
687			float P0;
688			float P1;
689			int dof;
690
691			// ---------------------------
692			// second : iteration
693			// ---------------------------
694			int niter = 0; // number of iterations
695			int pok[22];
696			// prova
697			// pesi lineari
698			float W22 = 1.;
699			float W1 = 10.;
700			float b = (W22 - W1) / (22. - 1.);// (w22 - w1)/(22-1)
701			float a = W1 - b;
702
703			for (int i = 0; i < 22; i++) {
704			pok[i] = 0;
705			}
706			do {
707			niter++;
708			float ttoll = toll; //tolerance (cm)
709			if (niter == 1)
710			ttoll = 10. * toll;
711			dof = GetN(); //previous fit
712			P0 = par[0];
713			P1 = par[1];
714			Reset(); //fit reset
715			for (int i = 0; i < 22; i++) {
716			cog[i] = 0;
717			qpl[i] = 0;
718			}
719			for (int ih = 0; ih < cevent->GetN(); ih++) {
720			float x = cevent->xycoord[ih];
721			float z = cevent->zcoord[ih];
722			float q = cevent->q[ih];
723			int ip = cevent->idp[ih];
724			float d = fabs(x - P0 - P1 * z) / sqrt(P1 * P1 + 1);
725			// cout<<d<<endl;
726			if ((niter == 1 \|\| (niter > 1 && pok[ip] == 1)) && (ip <= ipmax \|\| ip <= 3) &&
727			// ((ipmax>3 && ip<=ipmax) \|\| ipmax<=3) &&
728			true) {
729
730			if (niter == 1 && d < dmin[ip]) {
731			dmin[ip] = d;
732			closest[ip] = ih;
733			}
734
735			float W = a + b * (ip + 1);
736			// cout << ip << " "<<W<<endl;
737			if (d < ttoll \|\| (niter == 2 && ih == closest[ip])) {
738			/* if(ip==ipmax && ipmax>1)Add(z,x,q,Wq); /
739			/* else Add(z,x,q,W); */
740			cog[ip] += x * q;
741			qpl[ip] += q;
742			Add(z, x, q, W);
743			pok[ip] = 1;
744			cout << ip << " " << qpl[ip] << endl;
745			}
746			}
747			}
748			// break the track if more than 3 planes are missing
749			int nmissing = 0;
750			for (int ip = 0; ip < 22; ip++) {
751			if (pok[ip] == 0) {
752			nmissing++;
753			}
754			else {
755			nmissing = 0;
756			}
757			if (nmissing == 6) {
758			for (int ipp = ip + 1; ipp < 22; ipp++)
759			pok[ipp] = 0;
760			break;
761			}
762			}
763
764			for (int ip = 0; ip < 22; ip++)
765			cout << qpl[ip] << " ";
766			cout << endl;
767
768			for (int ip = 0; ip < 22; ip++)
769			if (qpl[ip] != 0)
770			cog[ip] = cog[ip] / qpl[ip];
771
772			if (niter == 100)
773			break;
774			if (GetN() < 3)
775			break;
776
777			Fit();
778
779			} while (niter == 1 \|\| GetN() != dof);
780			// cout << " >>> "<<GetN()<<endl;
781			if (GetN() < 3)
782			return 0;
783			Fit();
784			/* // cout << " n. hits :"<<GetN()<<endl; */
785			/* // cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl; */
786			/* // cout << " Chi2 :"<<GetChi2()<<endl; */
787			/* // cout << " Q axis :"<<GetQaxis()<<endl; */
788			/* // cout << " Q strip :"<<GetQstrip()<<endl; */
789
790			/* // --------------------------------------------- */
791			/* // third : selecting only closest strip-clusters */
792			/* // and fit baricenters */
793			/* // --------------------------------------------- */
794			/* P0 = par[0]; */
795			/* P1 = par[1]; */
796			/* Reset(); */
797
798			/* float dmin[22]; */
799			/* int closest[22]; */
800			/* for(int ip=0; ip<22; ip++){ */
801			/* dmin[ip]=999; */
802			/* closest[ip]=-1; */
803			/* } */
804			/* for(int ih=0;ih<cevent->GetN(); ih++){ */
805			/* float x = cevent->xycoord[ih]; */
806			/* float z = cevent->zcoord[ih]; */
807			/* float q = cevent->q[ih]; */
808			/* int ip = cevent->idp[ih]; */
809			/* int is = cevent->ids[ih]; */
810			/* float d = fabs(x-P0-P1z)/sqrt(P1P1+1); */
811			/* if( d < toll && d<dmin[ip] && pok[ip]==1 ){ */
812			/* // Add(z,x,q,q); */
813			/* closest[ip]=ih; */
814			/* // cog[ip] += xq; /
815			/* // qpl[ip] += q; */
816			/* } */
817			/* } */
818			/* for(Int_t ip=0; ip<22; ip++){ */
819			/* if(closest[ip]!=-1){ */
820			/* float x = cevent->xycoord[closest[ip]]; */
821			/* float z = cevent->zcoord[closest[ip]]; */
822			/* float q = cevent->q[closest[ip]]; */
823			/* int is = cevent->ids[closest[ip]]; */
824			/* Add(z,x,q,q); */
825			/* cog[ip] += xq; /
826			/* qpl[ip] += q; */
827			/* // cout << ip << " -o- "<<is<<endl; */
828			/* } */
829			/* } */
830			/* // add +/- one strip */
831			/* for(int ih=0;ih<cevent->GetN(); ih++){ */
832			/* float x = cevent->xycoord[ih]; */
833			/* float z = cevent->zcoord[ih]; */
834			/* float q = cevent->q[ih]; */
835			/* int ip = cevent->idp[ih]; */
836			/* int is = cevent->ids[ih]; */
837			/* if( closest[ip]!=-1 ){ */
838			/* int isc = cevent->ids[closest[ip]]; */
839			/* if( is == isc+1 \|\| is == isc-1 ){ */
840			/* Add(z,x,q,q); */
841			/* cog[ip] += xq; /
842			/* qpl[ip] += q; */
843			/* // cout << ip << " -+- "<<is<<endl; */
844			/* } */
845			/* } */
846			/* } */
847			/* Fit(); */
848
849			/* if(GetN()<3)return 0; */
850			/* for(int ip=0; ip<22; ip++){ */
851			/* if(qpl[ip]!=0){ */
852			/* cog[ip]=cog[ip]/qpl[ip]; */
853			/* // Add(z,cog[ip],cog[ip],0.7); */
854			/* } */
855			/* } */
856
857			// cout << " n. hits :"<<GetN()<<endl;
858			// cout << " P0 P1 :"<<par[0]<< " "<<par[1]<<endl;
859			// cout << " Chi2 :"<<GetChi2()<<endl;
860			// cout << " Q axis :"<<GetQaxis()<<endl;
861			// cout << " Q strip :"<<GetQstrip()<<endl;
862			// cout << " ---------------------------"<<endl;
863			return 1;
864			}
865			;
866
867			int FitShower(CaloLevel1* calo, Int_t view, Float_t toll) {
868			return FitAxis(calo, view, toll, 0);
869			}
870			;
871	pam-fi	1.1
872			};
873
874			//===============================================================================
875			//
876			//
877			//
878			//
879			//===============================================================================
880
881			#endif // CALOAXIS2_H_