#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>

#endif
//===============================================================================
//
//
//
//
//===============================================================================
class CaloAxis {

public:
    vector<float> x;// independent coordinates (z)
    vector<float> y;// dependente coordinate (x/y)
    vector<float> w;// fit weight
    vector<float> q;// signal

    Float_t par[2];
    Float_t covpar[2][2];

    TPolyMarker *sel;
    TLine       *axis;
    
    CaloAxis(){
	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;
    };

    CaloAxis(Int_t nin, float* xin, float *yin){
	for(Int_t i=0; i<nin; i++){
	    x.push_back(*xin++);
	    y.push_back(*yin++);
	    w.push_back(1.);    
	    q.push_back(0.);    
	};
	if(nin>2)Fit();
	sel  = 0;
	axis = 0;
    };
    CaloAxis(Int_t nin, float* xin, float *yin, float* qin){
	for(Int_t i=0; i<nin; i++){
	    x.push_back(*xin++);
	    y.push_back(*yin++);
	    q.push_back(*qin++);    
	    w.push_back(1.);    
	};
	if(nin>2)Fit();
	sel  = 0;
	axis = 0;
    };
    CaloAxis(Int_t nin, float* xin, float *yin, float* qin, float* win){
	for(Int_t i=0; i<nin; i++){
	    x.push_back(*xin++);
	    y.push_back(*yin++);
	    w.push_back(*win++);    
	    q.push_back(*qin++);    
	};
	if(nin>2)Fit();
	sel  = 0;
	axis = 0;
    };

    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_t Fit(){

	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 Reset(){
	for(Int_t i=0; i<2;i++){
	    par[i]=0;
	    for(Int_t j=0; j<2;j++)covpar[i][j]=0;
	};
	x.clear();
	y.clear();
	w.clear();
	q.clear();
	if(sel)sel->Delete();
	if(axis)axis->Delete();
	
    }

    void Delete(){
	Reset();
// 	delete x;
//	delete y;
// 	delete w;
//	delete q;
    }

    void Print(){
	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;
	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 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 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(){
	
	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);

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

    Int_t FitAxis( CaloLevel1* calo , Int_t view , Float_t toll ){

	CaloStrip* cstrip = new CaloStrip(calo);

	// ---------------------------
	// first : all hits included
	// ---------------------------
//	CaloAxis *axis = new CaloAxis();    
	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){
//	    cout<<ip<<" "<<is<<" "<< calo->DecodeEstrip(ih,iv,ip,is)<<" "<<endl;
// 	    float x = GetCoordXY(view,ip,is) ;
// 	    float z = GetCoordZ(view,ip) ;
// 	    float q = calo->GetEstrip(view,ip,is);
		cstrip->Set(view,ip,is);
		float x = 0;
		if( !view ) x = cstrip->GetX();
		else        x = cstrip->GetY();
		float z = cstrip->GetZ();
		float q = cstrip->GetE();	    
//	    if(view==0 && ip==21)q=0;
		if(q>0)this->Add(z,x,q);	    
	    }
	
	}
//  cout << "(1) hits "<<this->GetN()<<endl;
	if(this->GetN()<3){
	    delete cstrip;
	    return 0;
	}
	float P0;
	float P1;
	int   dof;
	bool  missing = false;
	int   nmissing = 0;

	// ---------------------------
	// second : iteration
	// ---------------------------
	int niter=0;                      // number of iterations
	do{
	    niter++;
	    float ttoll = toll;           //tolerance (cm)
	    if(niter==1)ttoll = 10.*toll;
	    dof = this->GetN();           //previous fit
	    P0  = this->par[0];
	    P1  = this->par[1];
	    this->Reset();                //fit reset
//	cout << "previous fit "<<dof<<" "<<P0<<" "<<P1<<endl;
//	cout << "(reset) hits "<<this->GetN()<<endl;
	    missing = false;
	    nmissing = 0;
	    for(Int_t ip=0; ip<22; ip++){
		// break the track if more than 2 planes are missing
		if(missing)nmissing++;  
		if(nmissing==3)break;
		missing=true;
//	    float z = GetCoordZ(view,ip) ;
		for(Int_t is=0; is<96; is++){

		    cstrip->Set(view,ip,is);
		    float x = 0;
		    if( !view ) x = cstrip->GetX();
		    else        x = cstrip->GetY();
		    float z = cstrip->GetZ();
		    float q = cstrip->GetE();

		    if( q ){
// 		    float x = GetCoordXY(view,ip,is) ;
// 		    float q = calo->GetEstrip(view,ip,is);
//		    if(view==0 && ip==21)q=0;
			float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1);
//		    cout << d << " "<<ttoll<<endl;
//		    if( q>0 && d < ttoll ){
			if( d < ttoll ){
			    this->Add(z,x,q);
//			cout << "(-) hits "<<this->GetN()<<endl;
			    missing=false;
			    nmissing=0;
			}
		    }
		}
	    }
//	cout << "(2) hits "<<this->GetN()<<endl;
	    if(niter==100)break;
	    if(this->GetN()==0)break;
	}while(niter==1 || this->GetN() != dof);	
	if(this->GetN()<3){
	    delete cstrip;
	    return 0;
	}
	// ---------------------------------------------
	// third : selecting only closest strip-clusters    
	// ---------------------------------------------
	P0  = this->par[0];
	P1  = this->par[1];
	this->Reset();
	missing = false;
	nmissing = 0;
	for(Int_t ip=0; ip<22; ip++){
	    if(missing)nmissing++;
	    if(nmissing==3)break;
	    missing=true;
	    int isclose = -1;
	    float dmin  = 999;
//	float z = GetCoordZ(view,ip) ;
	    //search for closest hit strip
	    for(Int_t is=0; is<96; is++){
		cstrip->Set(view,ip,is);
		float x = 0;
		if( !view ) x = cstrip->GetX();
		else        x = cstrip->GetY();
		float z = cstrip->GetZ();
		float q = cstrip->GetE();
		if( q ){
//		float x = GetCoordXY(view,ip,is) ;
		    float d = fabs(x-P0-P1*z)/sqrt(P1*P1+1);
//		cout << x <<" "<<z<<" "<<P0+P1*z<<" "<<d  << endl;
		    if(d<dmin && d<toll){
			isclose =is;
			dmin    =d;
		    }
		}
	    }
	    if( isclose != -1 ){
		Int_t isfrom = isclose-1;
		if(isfrom<0) isfrom=0;
		Int_t isto = isclose+2;
		if(isto>96) isto=96;
		for(Int_t is=isfrom; is< isto; is++){
		    cstrip->Set(view,ip,is);
		    float x = 0;
		    if( !view ) x = cstrip->GetX();
		    else        x = cstrip->GetY();
		    float z = cstrip->GetZ();
		    float q = cstrip->GetE();
//		float x = GetCoordXY(view,ip,is) ;
//		float q = calo->GetEstrip(view,ip,is);
//		if(view==0 && ip==21)q = 0; //mask last plane
		    if(q){
			this->Add(z,x,q);
			missing  = false;
			nmissing = 0;
		    }
		}
	    }	
	}
    
//     cout << " n. hits :"<<this->GetN()<<endl;
//     cout << " P0 P1   :"<<this->par[0]<< " "<<this->par[1]<<endl;
//     cout << " Chi2    :"<<this->GetChi2()<<endl;
//     cout << " Q axis  :"<<this->GetQaxis()<<endl;
//     cout << " Q strip :"<<this->GetQstrip()<<endl;
	delete cstrip;
	return 1;
    }
};
//===============================================================================
//
//
//
//
//===============================================================================