flyopunta/src/src_uti62/clusgmp.cpp

/***************************************************************************
 *   Questa File e di test
 ***************************************************************************/

#include "parm.h"
#include "apparato.h"
#include "evento.h"
#include "clusgmp.h"
#include "gmatrix.h"


using namespace GmpConst;
using namespace std;

// CONSTRUCTOR //

ClusterStraw::ClusterStraw ( DevStraw *stwfirst )
{

    Pij.Set_Size(3,3);
    Pjx.Set_Size(3,3);
    Pjy.Set_Size(3,3);
    Pju.Set_Size(3,3);
    Pjv.Set_Size(3,3);        ;


    // Reset some variables
    nclust=-1;
    nx=ny=nu=nv=0;
    ufl=vfl=xfl=yfl=0;

// chicut:  test per la bontà del fit a quattro fili e a tre fili. ( a due chiq è nullo per definizione.)
// la distanza media è d= sqrt (chicut/n) dove n= 3 o 4;
    chicut=0.8;

    // stwfirst is the first plan of the chamber.
    stwu =  stwfirst;
    stwv =  stwu->nextstraw;
    stwx =  stwv->nextstraw;
    stwy =  stwx->nextstraw;
    zy = stwy->Get_Cface().z;
    zx = stwx->Get_Cface().z;
    zv = stwv->Get_Cface().z;
    zu = stwu->Get_Cface().z;


    // Some output
//   Gout<<"\n  ClusterStraw for "<<stwu->Get_Nome() <<" done. ";
//   Gout<<"\n Piani z "<<  zu  <<" "<<  zv  <<" "<< zx  <<" "<<zy;

// calcolo i proiettori verticali degli assi:  ( Pj= I -|u><u|)



    gmatrix Unita(3,1.);
    gvet verso_cw;
    verso_cw.setvn(0.,1.,0.);  //filo nel lab lungo y
    verso_cw>>Pjx;
    Pjx=Unita-Pjx;
    verso_cw.setvn(1.,0.,0.);  //filo nel lab lungo x
    verso_cw>>Pjy;
    Pjy=Unita-Pjy;
    verso_cw.setvn(SINCOS45,-SINCOS45, 0.);
    verso_cw>>Pju;
    Pju=Unita-Pju;
    verso_cw.setvn(SINCOS45,SINCOS45, 0.);
    verso_cw>>Pjv;
    Pjv=Unita-Pjv;
    /*
        Pjx.Print("Pijx");
        Pjy.Print("Pijy");
        Pju.Print("Piju");
        Pjv.Print("Pijv");*/

}


// CONSTRUCT CLUSTERS                                        //
//                                                           //
// The mskcode gives the multiplicity for each value of mask //
// Mask values 0000,0001,0010,0011.... ( da 0 a 15)          //

// Proiettori degli assi dei 4 piani.

int ClusterStraw::Get_Cluster()
{
    nclust=-1;
    single_plans=0;
    ncls=0;
    static int codemask[16]={0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};
    static int coppmask[16]={0,0,0,0,0,1,1,0,0,1,1,0,0,0,0,0};
    static int LowPlan=0; // attenzione la loop sui piani parte da LowPlan+1 up to maxplan;

    int xflag[30],yflag[30],uflag[30],vflag[30];     // flags to kill the used wirecoord.
    for (int i=0;i<30;i++)xflag[i]=yflag[i]=uflag[i]=vflag[i]=1;// flags set to ok

    // Store number of hits per plane and the maximun active number of planes
    int maxplan=0;
    nx=ny=nu=nv=0;
//    if ( (nx = stwx->mhit)>0) maxplan++;
//    if ( (ny = stwy->mhit)>0) maxplan++;
//    if ( (nu = stwu->mhit)>0) maxplan++;
//    if ( (nv = stwv->mhit)>0) maxplan++;


    if ( (nx = stwx->nwire)>0) maxplan++;
    if ( (ny = stwy->nwire)>0) maxplan++;
    if ( (nu = stwu->nwire)>0) maxplan++;
    if ( (nv = stwv->nwire)>0) maxplan++;


// main loop
// prima si analizzano le associazioni a 4 viste, poi a tre ed infine a due



    if ( maxplan<1)  return 0;  // not enough views?  Un

    gvet Xcls;
    double chiq=0.;
    int Nplsvis=0;
    int ix,iy,iu,iv;
    int ixok,iyok,iuok,ivok;
    int id_x,id_y,id_u,id_v;

    id_x=id_y=id_u=id_v=0;
    ixok=iyok=iuok=ivok=0;

    for (int VisAcc=maxplan; VisAcc>LowPlan;VisAcc--)
    {
        int Istart=0;
        if (VisAcc<4)Istart=-1;

        for (iu= Istart;iu<nu;iu++)
        {
            iuok=0;
            id_u=0;
            if (iu>-1&&uflag[iu]>-1) {
                uwc.setvt(stwu->Wirecoor[iu],0.,0.,stwu->Wiredtim[iu]); 
                ufl=uwc.x;
                                                                udt=uwc.t;
                uwc=stwu->Devc2Lab(uwc);
                iuok=1;// primo piano
                id_u=stwu->M_Hits[iu].id;
            }
            else      uwc.Reset();

            for (iv=Istart;iv<nv;iv++)
            {
                ivok=0;
                id_v=0;
                if (iv>-1 && vflag[iv]>-1) {
                    vwc.setvt(stwv->Wirecoor[iv],0.,0.,stwv->Wiredtim[iv]);     
                    vfl=vwc.x;
                                                                                vdt=vwc.t;
                    vwc=stwv->Devc2Lab(vwc);
                    ivok=2;  // secondo piano
                    id_v=stwv->M_Hits[iv].id;
                }
                else      vwc.Reset();

                for (ix= Istart;ix<nx;ix++)
                {
                    ixok=0;
                    id_x=0;
                    if (ix>-1 && xflag[ix]>-1) {
                        xwc.setvt(stwx->Wirecoor[ix],0.,0.,stwx->Wiredtim[ix]); 
                        xfl=xwc.x;
                                                                                                xdt = xwc.t;
                        xwc=stwx->Devc2Lab(xwc);
                        ixok=4; // terzo piano  ok
                        id_x=stwx->M_Hits[ix].id;
                    }
                    else      xwc.Reset();
                    for (iy= Istart;iy<ny;iy++)
                    {
                        iyok=0;
                        id_y=0;
                        if (iy>-1 && yflag[iy]>-1) {
                            ywc.setvt(stwy->Wirecoor[iy],0.,0.,stwy->Wiredtim[iy]);
                            yfl=ywc.x;
                                                                                                                ydt=ywc.t;
                            ywc=stwy->Devc2Lab(ywc);
                            iyok=8;     // quarto piano ok
                            id_y=stwy->M_Hits[iy].id;
                        }
                        else      ywc.Reset();

// end for ==============

                        mask=iuok+ivok+ixok+iyok;
                        Nplsvis=codemask[mask];      // numero di piani visti
                        if (Nplsvis!=VisAcc)  continue; // no  enough active planes

// first special simple  test for cluster with only tho planes
// we accept only the couple u-x, u-y, v-x, v-y or mask = (0101,1001,0110,1010)=(5,9,6,10)

                        if (Nplsvis==2) if (coppmask[mask]==0) continue;


// reset variables
                        Xcls.Reset();
                        double Filo=0.;


                        if (VisAcc==1)  // only one plan is found active,
                        {
                            // special test for cluster with only one plan. It gives back the co0rdinate Xcls reported at Y=0,
                            if ( Filtro_1(Xcls,Filo)<0) continue; //skip if not ok;
                            chiq =-0.2;  // set to -0.2 to flag the cluster with only one plan
                            //                         Xcls.print("One_X");Gout<<" Ev "<< Eventi_Fatti;
                            single_plans++;

                        }

//=================  Calcolo del Punto Cluster con 2,3,4 piani =================
                        else if (VisAcc>1)
                        {
                            chiq =0.0;
                            gvet Noto;
                            Pij.Reset();;
                            //                      Noto.Reset();;

// calcolo della matrice dei coefficienti e il termine noto;
                            if (iyok>0) {
                                Pij+=Pjy;
                                Noto+=Pjy%ywc;
                            }
                            if (ixok>0) {
                                Pij+=Pjx;
                                Noto+=Pjx%xwc;
                            }
                            if (ivok>0) {
                                Pij+=Pjv;
                                Noto+=Pjv%vwc;
                            }
                            if (iuok>0) {
                                Pij+=Pju;
                                Noto+=Pju%uwc;
                            }

                            Pij.Det();


//                     Noto.Print("Noto");
//                     Pij.Print("Pij");

// calcolo le coord. del cluster
                            Xcls=(!Pij)%Noto;

// second  special test for cluster with only tho or three  planes
// the cluster found has to be in the black central region of the other not hit  plans...
// for two plans we accept only the couple u-x, u-y, v-x, v-y or mask = (0101,1001,0110,1010)=(5,9,6,10)


                            if (Nplsvis<4) if ( Filtro_23(Xcls)<0) continue; //skip if not ok;


// calcolo del chiq ( somme quadrate delle distanze del punto trovato  dai fili sul piano xy)
                            //per due piani chiq=0.0 per definizione.
                            if (Nplsvis>2) {
                                if (iuok>0)  chiq+=(Pju%(uwc-Xcls)).XYNorma();
                                if (ivok>0)  chiq+=(Pjv%(vwc-Xcls)).XYNorma();
                                if (ixok>0)  chiq+=(Pjx%(xwc-Xcls)).XYNorma();
                                if (iyok>0)  chiq+=(Pjy%(ywc-Xcls)).XYNorma();

                            }

                        }
//  attenzione end codice per 4-2 piani attivi
//=============================================
                        if (Nplsvis==4) chiq/=2.;

// memorizzo il cluster trovato se chiq/<cut;
                        if (chiq<chicut)
                        {
                            nclust++;
                            Clus[nclust].nome=stwu->Get_Nome();
                             Clus[nclust].Centro=stwu->Get_Centro();                          
                            Clus[nclust].id = nclust;
                            Clus[nclust].mask = mask;
                            Clus[nclust].hitx = id_x;
                            Clus[nclust].hity = id_y;
                            Clus[nclust].hitu = id_u;
                            Clus[nclust].hitv = id_v;
                            Clus[nclust].planes =Nplsvis;
                            if (ixok>0) {
                                Clus[nclust].Xfil = xfl;
                                                                                                                                Clus[nclust].Xdtim = xdt;
                                xflag[ix]=-1;
                                ixok=0;
                                id_x=0;
                            }
                            else  Clus[nclust].Xfil = 0.;
                            if (iyok>0) {
                                Clus[nclust].Yfil = yfl;
                                                                                                                         Clus[nclust].Ydtim = ydt;
                                yflag[iy]=-1;
                                iyok=0;
                                id_y=0;
                            }
                            else  Clus[nclust].Yfil = 0.;
                            if (iuok>0) {
                                Clus[nclust].Ufil = ufl;
                                                                                                                                Clus[nclust].Udtim = udt;
                                uflag[iu]=-1;
                                iuok=0;
                                id_u=0;
                            }
                            else  Clus[nclust].Ufil = 0.;
                            if (ivok>0) {
                                Clus[nclust].Vfil = vfl;
                                                                                                                                Clus[nclust].Vdtim = vdt;
                                vflag[iv]=-1;
                                ivok=0;
                                id_v=0;
                            }
                            else  Clus[nclust].Vfil = 0.;

                            Clus[nclust].chiq = chiq;
                            Clus[nclust].Xcls=Xcls;
                            Clus[nclust].Filo=Filo; //for one plan
                             }
 
                    }

                }
            }
        }
    }


    if (nclust>-1)Cluster2Root();
    
    
    
//     if (nclust>-1)
//     {
//       if(Clus[0].mask==15)
//       { gvet xpar,xxdv;
//            Gout<<"\n Nome "<<stwu->Get_Nome()<<" Centro "<<stwu->Get_Centro().x;
//            xxdv=Clus[0].Xcls;
//            xxdv.print("Xcls");
//            xpar=stwv->Lab2cDev(xxdv);
//            xpar.print("xv");
//             
//              Gout<< " fili  " << Clus[0].Ufil<<" "<<Clus[0].Vfil<<" "<<Clus[0].Xfil<<" "<<Clus[0].Yfil;
//           
//           Gout<<"\n filo x " << Clus[0].Xfil+stwx->Get_Centro().x<<" "<<xxdv.x;
//           Gout<<" filo y " << Clus[0].Yfil<<" "<<-xxdv.y;
//   
//           Gout<<" filo u " << Clus[0].Ufil+stwu->Get_Centro().x*SINCOS45<<" "<<xxdv.x*SINCOS45 +xxdv.y*SINCOS45;
//           
//           Gout<<" filo v " << Clus[0].Vfil+stwu->Get_Centro().x*SINCOS45<<" "<< xxdv.x*SINCOS45 -xxdv.y*SINCOS45 ;
//     
//         
//       } }
//       
//       
//     }
//     if (zu>21000.)
//     {
//
//       for (int i=0;i<nclust+1;i++)Clus[i].Cluster_Print ();
//     }

    return  nclust;

}
//======  Cluster2Root()
void ClusterStraw::Cluster2Root()
{

    ncls=nclust+1;
    nwu=nu;
    nwv=nv;
    nwx=nx;
    nwy=ny;

    for (int in=0;in<nclust+1;in++)
    {

        hitu[in]=    Clus[in].hitu;
        hitv[in]=    Clus[in].hitv;
        hitx[in]=    Clus[in].hitx;
        hity[in]=    Clus[in].hity;
        maskera[in]= Clus[in].mask;
        planes[in]=  Clus[in].planes;
        Ufil[in]=  Clus[in].Ufil;
        Vfil[in]=  Clus[in].Vfil;
        Xfil[in]=  Clus[in].Xfil;
        Yfil[in]=  Clus[in].Yfil;
                        Udtim[in]=  Clus[in].Udtim;
        Vdtim[in]=  Clus[in].Vdtim;
        Xdtim[in]=  Clus[in].Xdtim;
        Ydtim[in]=  Clus[in].Ydtim;
        X[in]=  Clus[in].Xcls.x;
        Y[in]=  Clus[in].Xcls.y;
        Z[in]=  Clus[in].Xcls.z;
        chiq[in]=  Clus[in].chiq;
    }

}

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



void HitClus::Cluster_Print()
{

    Gout.setf ( ios::right,ios::adjustfield );
    Gout<<std::setfill ( ' ' );
    Gout<<std::setprecision ( 3 ) <<fixed;
    Gout<<"\n Ev-clus "<<setw ( 6 ) <<evento_.Gen.Event<<" "
    <<nome
    <<" Id "<<setw ( 2 ) <<id
    <<" Msk "<<setw ( 8 ) <<mask
    << " hts ["<<setw ( 6 ) <<hitu
    <<" "<<setw ( 6 ) <<hitv
    <<" "<<setw ( 6 ) <<hitx
    <<" "<<setw ( 6 ) <<hity
    <<"] W "<<setw ( 8 )<<Ufil
    <<" "<<setw ( 8 ) <<Vfil
    <<" "<<setw ( 8 ) <<Xfil
    <<" "<<setw ( 8 ) <<Yfil;
    Xcls.print("X");
    Gout<<" Chiq "<<setw ( 8 ) <<chiq;

    Gout.setf ( ios::floatfield );


}
//=======================================
int  ClusterStraw::Filtro_23(gvet Xcls)
{
    static int Central_Gap= 13.5/2.;  // attenzione questo parametro genera le righe nel plot dei xy dei cluster....(stringere un po')
    // attenzione vale solo per 2 o 3 camere attive
    static int short msk1=1, msk2=2, msk3=4, msk4=8;
    unsigned short maskera=mask;

    gvet Xcldv;
    if ((maskera&msk1)==0)   // camera u non c'è
    {
        Xcldv=stwu->Lab2cDev(Xcls);  // passo al rif interno al dev
        if ( abs(Xcldv.x)>Central_Gap) return -1; // se la x è fuori buco salto...
    }
    if ((maskera&msk2)==0)   // camera v non c'è
    {
        Xcldv=stwv->Lab2cDev(Xcls);
        if ( abs(Xcldv.x)>Central_Gap) return -1;
    }
    if ((maskera&msk3)==0)   // camera x non c'è
    {
        Xcldv=stwx->Lab2cDev(Xcls);
        if ( abs(Xcldv.x)>Central_Gap) return -1;
    }
    if ((maskera&msk4)==0)   // camera y non c'è
    {
        Xcldv=stwy->Lab2cDev(Xcls);
        if ( abs(Xcldv.x)>Central_Gap) return -1;
    }
    return 1;
}

int  ClusterStraw::Filtro_1(gvet &X,double &Filo)
{// filtro to test the "one plan" case: cut on wire max address
// filtro solo valido per il caso di un solo piano colpito
//restituisce la coordinata del punto nel lab come se nel dev avesse
// coordinate( x=filo,  y=0)
    static int Maxwire_addr= 9.018;
    unsigned short maskera=mask;
    static int short msk1=1, msk2=2, msk3=4, msk4=8;
    DevStraw *pstw=0;

    if ((maskera&msk1)>0) {
        pstw=stwu;
        Filo=ufl;
    }
    else if ((maskera&msk2)>0) {
        pstw=stwv;
        Filo=vfl;
    }
    else if ((maskera&msk3)>0) {
        pstw=stwx;
        Filo=xfl;
    }
    else if ((maskera&msk4)>0) {
        pstw=stwy;
        Filo=yfl;
    }
    else return -1; // se c'è un errore.....

    if (abs(Filo)>Maxwire_addr) return -1; // indirizzo fuori zona....
    X.setvn(Filo,0.,0.);
    X=pstw->Devc2Lab(X);
    X.z=pstw->Get_Centro().z;

    return 1;

}






//=======================================
void  ClusterStraw::Ordina()
{
    // ordina i cluster  in dist decrescente...
    double chiq,chiqb;

    if ( nclust>0 ) // ordina i cluster
    {
        for ( int i=0;i<nclust;i++ )
        {
            chiq=Clus[i].chiq;
            for ( int j=i+1;j<nclust;j++ )
            {
                if ( ( chiqb=Clus[j].chiq ) <chiq )
                {
                    ClusSave=Clus[i];
                    Clus[i]=Clus[j];
                    Clus[j]=ClusSave;
                    chiq=chiqb;
                }
            }
        }
    }
}