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_

---