#include <stdlib.h>
#include <iostream.h>
#include <math.h>

// standard definition of NULL and bool
const int null = 0;
#define bool int
#define true 1
#define false 0
#define DEBUG 1


//definition of contants for SiBT geometry
//

const int       nplane = 8;                     // number of planes in SiBT
const float     sbt_ds = 200.;                  // length of SiBT over s-axis
const float     sbt_dx = 40.;                   // width of SiBT over x-axis
const float     sbt_dy = 38.;                   // height of SiBT over y-axis
const float     length = .003;                  // length of strip over s-axis
const int       nstrip = 1024;                  // number of strips
const float     xwidth = .9*sbt_dx/nstrip;      // width of strip in x direction
const float     ywidth = .9*sbt_dy/nstrip;      // width of strip in y direction

const int       hitquot = 3;                    // hit is represented by val/noise > hitquot
const float     linchi  = 10.;                  // chi of accepted line fit

// definition of classes
//

template<class T> class list{
// parametrized class List
// singly linked list of items
// initialized with zero items, with head and cur initialized to null
protected:
        class Link{
        public:
                Link*  next;
                 T*     val;
                Link(Link* n, T* v) {next = n; val =v; } ;
        };
        Link*   head;   // head of singly linked list
        Link*   cur;    // last selected item
public:
        int     number; // number of items in list
        list() {head = null; cur=null; number=0; };
        ~list(){ 
        while ( head != null){
                cur = head;
                delete head->val;
                head = head->next;
                delete cur;}
        };
void    enter(T* item) {
        Link* temp = head;
        if (item != null){
                head = new Link( temp, item);
                number ++;}
        };
T*      first() {cur=head; 
        if (cur != null) {return (cur->val);}
        else { return null;}
        };
T*      next() {
        if (cur != null){ 
                cur = cur->next;
                if (cur != null){return (cur->val);}
                else {return null;}}
        else {return null;}
        };
T*      current() { 
                if (cur != null){ return (cur->val);}
                else {return null;}
         };
bool    sequel(){
        if (cur == null){ return false;}
        else {return cur->next != null;}
        };
void    remove(T* item) {
        Link**    temp = &head;  // temp contains the location of the pointer to the item
                                // temp is initialized with the location of head
        for (Link** temp = &head; *temp == null; temp = &((*temp)->next)){ 
                if ((*temp)->val == item){
                        Link*  rm = *temp;
                        delete item;
                        delete rm;
                        temp = &((*temp)->next);
                        if (cur == *temp) cur = null;}  // cur points to removed item,
                                                        // so cur is set to null
                }
        };      
#ifdef DEBUG
void    print(){
        cout << " list with " << number <<" entries \n" << flush;
        int k =0;
        Link* temp = head; 
        while (temp != null){ 
                if ( cur == temp) { cout << "current ";}
                else { cout << "        ";};
                cout << " item "<< ++k << ":   " << flush;
                temp->val->print();
                temp = temp->next;}
        };
#endif

};

class coord{
public:
float   x,y,z;
};

class channelset{
int  items[50];
int  as;
public:
        channelset(){
                as = 0;
                for( int i =0; i<50; ++i) {items[i] = 0;}
        }

int     strip(int pl, int st) {
                for (int i = 0; i<as ; i = i+3){
                if (items[i] == pl){ 
                        return items[i+st];} // hit of channel pl found
        }
        return 1;}                      // no hit found 
                                        // (minimum noise of 1 to prevent div by zero

void    assign( int strip, int value, int noise){
        items[as] = strip; as++;
        items[as] = value; as++;
        items[as] = noise; as++;}
#ifdef DEBUG
void    print(){
        int i=0;
        while (i < as){
                cout << "          strip " << items[i] <<";  value " << 
items[i+1] <<";   noise " << items[i+2] << "\n" << flush;
                i = i+3;}
        }
#endif
};


class event{
channelset      set[8]; // 1 event: 8 channel sets for 8 planes
public:
        event(){
        for (int pl = 0; pl <8 ; ++pl){

// if (pl != 2) {
                set[pl].assign(40+abs(pl-4),20,4);
                set[pl].assign(41+abs(pl-4),21,5);
                set[pl].assign(42+abs(pl-4),35,6);
//              set[pl].assign(254+abs(pl-4),59,8);
                set[pl].assign(255+abs(pl-4),63,9);}
//           }
/*
if (pl != 2){
                set[pl].assign(40,20,4);
                set[pl].assign(41,21,5);
                set[pl].assign(42,35,6);
//              set[pl].assign(254,59,8);
                set[pl].assign(255,63,9);}
             }
*/
        }

channelset plane(int pl){ return set[pl];}
#ifdef DEBUG
void    print(){
        cout << " EVENT print \n";
        for (int i =0; i<8; ++i){
                cout << "   Plane " << i+1 << " \n" << flush;
                set[i].print();}
        cout<< " --------- End Event \n \n" << flush;
        }
#endif
};

class strip {
public:
int noise;
int value;
float dist;    //distance from origin 
float width;    // width of strip
float length;    // length of strip over s-axis
};

class cluster{
protected:
coord   cc;             //coordinates of cluster
public:
float   size;           // number of hits (>1)
float   weight;         // weight =0 for dummy cluster
float   s_c;            // coordinate along s_axis (z_axis)
float   t_c;            // coordinate along transverse axis (x_ or y_ axis)
        cluster(float cx,float cy,float cs,float sz)
                { cc.x = cx; cc.y = cy; cc.z = cs; weight =1.; size = sz; }
#ifdef DEBUG
virtual void    print(){
        cout << "Cluster " << "t: " << t_c << "   s: " 
                        << s_c << "  size: " << size << " \n" << flush;
        }
#endif
};

class h_cluster: public cluster{
public:
        h_cluster(float cx,float cs,float sz)
                :cluster(cx,0.,cs,sz) {         // y-coordinate =0
                t_c = cx; s_c = cs; }
#ifdef DEBUG
void    print(){
        cout << "Horizontal "  << flush; cluster::print();
        }
#endif
};


class v_cluster: public cluster{
public:
        v_cluster(float cy,float cs,float sz)
                :cluster(0.,cy,cs,sz) {         // x-coordinate = 0
                t_c = cy; s_c = cs; }
#ifdef DEBUG
void    print(){
        cout << "  Vertical "  << flush; cluster::print();
        }
#endif
};

class plane {
protected:
list<cluster>*  clusters;       // cluster set
strip           sl[nstrip];     // array of strips
float           x,y,s;          // coordinates of lower left corner of plane
float           mean;              // statistics
int             nr;             // number of entries

public:
        plane(coord c){
        clusters = new list<cluster>;   // empty cluster list to allow access to clusters
        x = c.x; y=c.y ; s = c.z;       // coordinates of plane
        mean = 0.; nr = 0;              //statistics for calibration
        for (int i = 0; i<nstrip; ++i) { 
                sl[i].value = 0;
                sl[i].noise = 1;}               // to prevent division by zero
        }
cluster* first(){return clusters->first();};
cluster* get(){return clusters->current();};
cluster* next(){return clusters->next();};
bool    sequel(){ return clusters->sequel();}
void    assign(channelset s);
void    clearstatistics(){ mean = 0.; nr = 0;};
float   meandev(){ return mean;};
float   scoord() { return s;};
virtual void    storealignment(float cc){  };
void    adddeviation(float dt){ 
        if (nr == 0) {mean = dt; nr = 1;}
        else { mean = mean*(nr/(nr+1)); nr ++; mean = mean+(dt/nr);} 
        };
virtual void addcluster(float t_c, float s_c, float sz) { };
#ifdef DEBUG
virtual void    print(){
        cout << "PLANE coordinates:: " << "x: " << x << "  y: " << y << 
                                        "  s: " << s << " \n";
        for (int i = 0; i<nstrip; ++i) { 
           if (sl[i].value > 1){
                cout << "     strip: " << i << " value: " << sl[i].value  <<
                         "  noise: " << sl[i].noise << "  distance: " << 
                                sl[i].dist << " \n" << flush;}
           }
        }
#endif
};

 
void plane::assign(channelset ss){
        for (int i = 0; i< nstrip; ++i) {  // strip value and noise from channelset
                sl[i].value = ss.strip(i,1);
                sl[i].noise = ss.strip(i,2);}   // strips initialized

        delete clusters;                        // remove former clusters
        clusters = new list<cluster>;           // empty cluster list
        for (i = 0; i< nstrip; ++i) {  // determine clusters
                if (sl[i].value/sl[i].noise > hitquot){
                        float   total =sl[i].value-sl[i].noise;
                        float   posit = sl[i].dist*total;
                        int  k = i+1;
                        bool b = (k < nstrip);
                        if (b) {b = b & sl[k].value/sl[k].noise > hitquot;} 
                        while (b) {
                            total =total+(sl[k].value -sl[k].noise);
                            posit = posit+sl[k].dist *float( 
                                        sl[k].value-sl[k].noise); 
                            k++;
                            b = (k < nstrip);
                            if (b) {b = b & sl[k].value/sl[k].noise > hitquot;} 
                            }
                        k--;            // k always one too large
                        if (k > i) {
                                //more than two consecutive hits, add cluster
                                addcluster(posit/total, s,
                                        sl[k].dist-sl[i].dist+sl[k].width);}
                        i = k+1;}               //skip already inspected strips
                }
}


class h_plane: public plane{
public:
        h_plane(coord cc): plane(cc) {  
        for (int i = 0; i<nstrip; ++i) { 
                sl[i].dist = xwidth*float(i);   // x-position of strips from x=0
                sl[i].width = xwidth;           // dx width of strips
                sl[i].length = length;}         // thickness of strip over s-axis
        }
void    addcluster(float t_c, float s_c, float sz) {
                // add x coordinate of plane to position within plane
                h_cluster* h_cl = new h_cluster(t_c+x, s_c, sz);
                clusters->enter( h_cl);
        }
void    storealignment(float cc){ x = x+cc;};
#ifdef DEBUG
void    print(){
        cout << " HORIZONTAL ";
        plane::print();
        cout << "     clusters are: \n";
        clusters->print();
        cout << " --------------------- H_PLANE \n \n" << flush;
        }
#endif
};
        


class v_plane: public plane {
public:
        v_plane(coord cc): plane(cc) {
        for (int i = 0; i<nstrip; ++i) { 
                sl[i].dist = ywidth*float(i);   // y-position of strips from y=0
                sl[i].width = ywidth;           // y-width of strips
                sl[i].length = length;}  //length of strip over s-axis
        }
void    addcluster(float t_c, float s_c, float sz) {
                // add y coordinate of plane to position within plane
                v_cluster* v_cl = new v_cluster(t_c+y, s_c, sz);
                clusters->enter( v_cl);
        }       
void    storealignment(float cc){ y = y+cc;};
#ifdef DEBUG
void    print(){
        cout << "   VERTICAL ";
        plane::print();
        cout << "     clusters are: \n";
        clusters->print();
        cout << " --------------------- V_PLANE \n \n" << flush;
        }
#endif

};


class s_line{
public:
float   s_off;
float   ts_tg;
float   chi;
        s_line(float xt,float yt,float of){
                s_off = of; ts_tg = yt; chi = 0;}
bool    fit(list<cluster>* proj) {
        cluster* clpt = proj->first();
        int     ncl =0;                                 // number of clusters
        float   S_s = 0.;
        float   S_ss = 0.;
        float   S_t = 0.;
        float   S_tt = 0.;
        float   S_st = 0.;
        while ( clpt != null){
                if (clpt->weight > 0){
                        ncl = ncl + 1;                  //  number of entries
                        S_s = S_s+clpt->s_c;                    //sum over s
                        S_t = S_t+clpt->t_c;                    //sum over t
                        S_ss = S_ss+clpt->s_c*clpt->s_c;        //sum over s*s
                        S_st = S_st+clpt->s_c*clpt->t_c;        //sum over t*s
                        S_tt = S_tt+clpt->t_c*clpt->t_c;}       //sum over t*t
                clpt = proj->next();
                }
        if (ncl < 3){ return false;}    // at least three coordinates required
        float det = float(ncl)*S_ss - S_s*S_s ;
        if (det <= 0) { return false;}
        s_off = (S_ss*S_t - S_s*S_st)/det;
        ts_tg = (ncl*S_st - S_s*S_t)/det;
        chi = S_tt-s_off*S_t - ts_tg*S_st;
        return abs(chi) < linchi;
        }
#ifdef DEBUG
void    print(){
        cout << "LINE  ts_tg: " << ts_tg << "  s_off: " 
                << s_off << "  chi: " << chi << " \n" << flush;}
#endif
};


class projection{
protected:
list<cluster>*  clp;    // pointer associated cluster
s_line  line;   // fitted line
public:
        projection(): line(0,0,0) { clp = new list<cluster>; }
        ~projection(){ delete clp;}
bool    valid(){ return line.fit(clp);}
float   tline( float sc){ return line.ts_tg*sc + line.s_off;}
void    add(cluster* clref) { if (clref != null) {clp->enter(clref);}}
void    AddDev( plane* planes[]){
        cluster* cp = clp->first();
        int i = (nplane/2)-1;
        while (cp != null) { // assume that clusters come from 4,3,2,1
                while (cp->s_c != planes[i]->scoord() & i >0) { i--;} 
                planes[i]->adddeviation(cp->t_c - tline(cp->s_c));
                cp = clp->next();} // end while over clusters
        }

#ifdef DEBUG
virtual void    print(){ 
        cout << " projection \n" << flush;
        clp->print();
        line.print();}
#endif
};

class v_projection: public projection{
public:

#ifdef DEBUG
void    print () {
        cout <<"   VERTICAL ";
        projection::print();
        }
#endif
};

class h_projection: public projection{
public:

#ifdef DEBUG
void    print () {
        cout <<" HORIZONTAL ";
        projection::print();}
#endif
};

class projlst {
protected:
list<projection>*  lp;
        projlst() { lp = null;}
        ~projlst() { if (lp != null) { delete lp;} }
public:
projection*     first() { return lp->first();}
virtual projection*     newproj() { return new projection();}
int     number() { return lp->number;}
virtual void    FindProj(plane* planes[]){
// find the list of projections in array of planes
        
        bool b = false;
        projection*  pp = newproj();
        for (int i = 0; i<4; ++i){
                pp->add(planes[i]->first());
                b = b | planes[i]->sequel();}  // more than one cluster in plane?
                // b implies there is a plane with more than 1 cluster
                // not b implies all planes have less than 2 clusters
        if (pp->valid()){lp->enter(pp);} // if good projection add to lvp
        else {delete pp;} // else remove contents and memory

// try all possible cluster combinations from vertical planes that have clusters

        cluster* cp;
        while (b){

                pp = newproj();
                b = false;
                bool sw = true;

        // sw implies get next cluster from plane or first one
        // not sw implies get current cluster from plane

                for (i = 0; i<4; ++i){
                        if (sw) {cp = planes[i]->next();}
                        else { cp = planes[i]->get();}
                        if (cp == null) { cp = planes[i]->first();}
                        else { sw = false;}

// sw remains true when when first cluster was taken from former plane
                        pp->add(cp);
                        b = b | planes[i]->sequel();}  // end for (i=0; i<4)

        // not b implies last cluster is found in all planes

                if (pp->valid()){lp->enter(pp);}
                else {delete pp;}
           }                            // end of while(b)
        }                               // end of FindProj

#ifdef DEBUG
virtual void    print() {
        cout << "projection list \n";
        if (lp==null) {cout << " no list present \n";}
        else {lp->print();}
        }
#endif
};

class v_projlst: public projlst {
protected:
projection*     newproj(){ return new v_projection();}
public:
void    FindProj(plane* planes[]){
// find the list of projections in array of vertical planes
        
        delete lp;                      // remove old list
        lp = (list<projection>*) new list<h_projection>; 
        projlst::FindProj( planes);
        }

#ifdef DEBUG
void    print() {
        cout << "\n   Vertical ";
        projlst::print();
        }
#endif
};


class h_projlst: public projlst {
protected:
projection*     newproj(){ return new h_projection();}
public:
void    FindProj(plane* planes[]){
// find the list of projections in array of vertical planes
        
        delete lp;                      // remove old list
        lp = (list<projection>*) new list<h_projection>;
        projlst::FindProj( planes);
        }

#ifdef DEBUG
void    print() {
        cout << "\n Horizontal ";
        projlst::print();
        }
#endif
};


class track{
projection*   hpr;
projection*   vpr;
public:
         track(projection* h, projection* v) {hpr = h; vpr = v;}
         ~track() { delete hpr; delete vpr;}
projection* hproj() {return hpr;}
projection* vproj() {return vpr;}

#ifdef DEBUG
void    print(){
        cout << " TRACK \n" << flush;
        if (hpr != null) {hpr->print();}
        else { cout << " no horizontal projection \n" << flush;}
        if (vpr != null) {vpr->print();}
        else { cout << " no vertical projection \n" << flush;}
        }
#endif
};


class SiBTOO{
plane*  hplanes[nplane/2];      // 4 horizontal planes
plane*  vplanes[nplane/2];      // 4 vertical planes
h_projlst*   lhp;               // list of horizontal projections
v_projlst*   lvp;               // list of vertical projections



public:
        SiBTOO() {      
                
        // initialize with 0 projections
        lhp = new h_projlst();
        lvp = new v_projlst();

        //initialize plane coordinates
        for (int i = 0; i<nplane/2; ++i){
                coord cc;
                cc.z = 2*i*(sbt_ds-0.05)/nplane;// s-coordinate
                cc.x = .1;                      // x-coordinate
                cc.y = cc.x;                    // y coordinate
                vplanes[i] = new v_plane(cc);   // vertical plane 
                cc.z = cc.z+.05;                // s-coordinate
                                                // (slightly different from 
                                                // vertical plane s-coordinate) 
                hplanes[i] = new h_plane(cc);}  // horizontal plane 
        }

track*  FindTracks(event ev){
        for (int i = 0; i<nplane/2; ++i){
                vplanes[i]->assign(ev.plane(i*2));   // vertical plane strips assigned
                hplanes[i]->assign(ev.plane(i*2+1));};// horizontal plane strips assigned
                // clusters of event are found and stored in planes 

        // VERTICAL projections
        lvp->FindProj( vplanes);        // lvp contains vertical projections

        // HORIZONTAL projections
        lhp->FindProj( hplanes);        // lhp contains horizontal projections

        // only one projection from each orientation is allowed
                if(lhp->number() == 1 & lvp->number() == 1) 
                                {return new track(lhp->first(), lvp->first());}
                else { return null;}
        }

void    Calibrate(list<track>* lst){
        for (int i = 0; i<nplane/2; ++i){
                vplanes[i]->clearstatistics();  
                hplanes[i]->clearstatistics();}
        track*  trp = lst->first();     // get first track from list
        while (trp != null) {           // loop over all tracks

        // horizontal projections
                trp->hproj()->AddDev( hplanes);
        // vertical projections
                trp->vproj()->AddDev( vplanes);

                trp = lst->next();} // end while over all tracks
        // realign planes 
        for (i = 0; i<nplane/2; ++i){
                hplanes[i]->storealignment(-hplanes[i]->meandev());
                vplanes[i]->storealignment(-vplanes[i]->meandev());}
        }
        

#ifdef DEBUG
void    print(){
        lhp->print();
        lvp->print();
        for (int i = 0; i<nplane/2; i++){
                vplanes[i]->print();
                hplanes[i]->print();}
        }
#endif
};


int main()
{
SiBTOO  sbt;
event   ev;             // some well chosen values
list<track>*    lst = new list<track>;  
                        //initialize SiBTOO
#ifdef DEBUG
        ev.print();             // print event contents
        sbt.print();            // print all plane contents of SiBTOO
        cout << "\n Track list first time \n\n\n" << flush;
        lst->print();           // print all tracks
#endif
        lst->enter( sbt.FindTracks(ev));        // find one track
#ifdef DEBUG
        sbt.print();
        cout << "\n Track list 2nd time\n\n\n" << flush;
        lst->print();
#endif
        sbt.Calibrate(lst);             // calibrate with one track
#ifdef DEBUG
        sbt.print();
#endif

        lst->enter( sbt.FindTracks(ev));        // find same track in calibrated SiBTOO
#ifdef DEBUG
        sbt.print();
        cout << "\n Track list 3rd time\n\n\n" << flush;
        lst->print();
#endif

cout << " main stopped\n";
}