flyopunta/src/src_uti62/trgmpstraw.cpp

/***************************************************************************
 *   Copyright (C) 2005 by Giuseppe Pierazzini          Pisa 20.02.05      *
 *   pierazzini@unipi.it                                                   *
 *                                     F l y o                             *
 *   Epsi/NA48                                                             *
 ***************************************************************************/
/* Authors G.Pierazzini - G.Ruggiero */

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


using namespace GmpConst;

using namespace std;

int errordebug=-1;

HitClus Hitp[4];

// DEFINITION OF CONSTANT VARIABLES

//const double TrgmpStraw::HIGHLEVEL = 5.;

// CONSTRUCTOR & DESTRUCTOR //
TrgmpStraw::TrgmpStraw()
{

    ntrk=-1;
    ncls0=ncls1=ncls2=ncls3=0;

    if (Mnp33!=0) z_centro_mag= Mnp33->Get_Centro().z;
    else  {cout<<"\n Error... Mnp33 not defined ....exit "<<endl; 
         Gout<<"\n Error... Mnp33 not defined ....exit "<<endl; 
         exit(0);}
    
    Gout<<"\n\n Class TrgmpStraw done with  Z_mag "<< z_centro_mag;


    // Pointers to the cluster in each group of chamber
    pcl0=ChambClus[0];
    pcl1=ChambClus[1];
    pcl2=ChambClus[2];
    pcl3=ChambClus[3];
}
TrgmpStraw::~TrgmpStraw()
{}



// MAKE TRACK                                                //
// NOTA
void TrgmpStraw::Make_Trk()
{

    static int codemask[16]={0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4};

    int c0ok,c1ok,c2ok,c3ok, nchamb,nsingle;
    int c0flag[30],c1flag[30],c2flag[30],c3flag[30];     // flags to kill the used cluster
    for (int i=0;i<30;i++)c0flag[i]=c1flag[i]=c2flag[i]=c3flag[i]=1;// flags set

    double chicut=3.;
    ntrk=-1;


    // Find the  Cluster number  in straw chambers..

    ncls0=ncls1=ncls2=ncls3=0;
    if ( pcl0!=0 ) ncls0=pcl0->nclust+1;
    if ( pcl1!=0 ) ncls1=pcl1->nclust+1;
    if ( pcl2!=0 ) ncls2=pcl2->nclust+1;
    if ( pcl3!=0 ) ncls3=pcl3->nclust+1;


    int maxcamere=0;// camere coinvolte.
    if (ncls0 >-1) maxcamere++;
    if (ncls1 >-1) maxcamere++;
    if (ncls2 >-1) maxcamere++;
    if (ncls3 >-1) maxcamere++;
    if (maxcamere <3) return; // minimum 3 chambers

    int Single_esiste=pcl0->single_plans+pcl1->single_plans+pcl2->single_plans+pcl3->single_plans;


    for (int VisCam=maxcamere; VisCam>2;VisCam--)
    {
        for (nsingle=0;nsingle<VisCam-1;nsingle++)   // single plan loop
        {
            if (Single_esiste==0&&nsingle>0) continue;

            int Istart=0;
            if (VisCam <4)Istart=-1;

            for ( int ic0=Istart;ic0<ncls0;ic0++)
            {
                c0ok=0;
                if (ic0>-1&&c0flag[ic0]>-1)
                {
                    Hitp[0]= pcl0->Clus[ic0];
                    c0ok=1;// prima camera
                }

                else     Hitp[0].mask=Hitp[0].planes=0;

                for ( int ic1=Istart;ic1<ncls1;ic1++)
                {
                    c1ok=0;
                    if (ic1>-1&&c1flag[ic1]>-1)  {
                        Hitp[1]= pcl1->Clus[ic1];
                        c1ok=2;  // seconda camera

                    }
                    else   Hitp[1].mask=Hitp[1].planes=0;


                    for ( int ic2=Istart;ic2<ncls2;ic2++)
                    {
                        c2ok=0;
                        if (ic2>-1&&c2flag[ic2]>-1)  {
                            Hitp[2]= pcl2->Clus[ic2];
                            c2ok=4; // terza camera
                        }
                        else    Hitp[2].mask=Hitp[2].planes=0;

                        for ( int ic3=Istart;ic3<ncls3;ic3++)
                        {
                            c3ok=0;
                            if (ic3>-1&&c3flag[ic3]>-1)  {
                                Hitp[3]= pcl3->Clus[ic3];
                                c3ok=8;
                            }
                            else     Hitp[3].mask=Hitp[3].planes=0;

                            short int mask_cam;
                            mask_cam=c0ok+c1ok+c2ok+c3ok;
                            nchamb=codemask[mask_cam];
                            if (nchamb!=VisCam)  continue; // no  enough active chamber

//verifica della configurazione per trattare i cluster con un piano
                             TrkNow.used_plans=0;
                            if (Single_esiste>0)
                                if (Verify(mask_cam,nsingle)<0) continue;

// fitting  data

                            // double chiq= Fitting(mask_cam,Hitp);

                            double chiq= Fitting(mask_cam);

                            if (Eventi_Fatti==errordebug)
                                Gout<<"\n Ev "<< Eventi_Fatti <<" chiq "<<chiq
                                <<" camere "<< VisCam
                                <<" Single_esite "<<Single_esiste <<" Single "<<nsingle

                                <<" msk "<< Hitp[0].mask  <<" "<< Hitp[1].mask  <<" "<< Hitp[2].mask  <<" "<< Hitp[3].mask
                                <<" Ic.. "<<ic0<<" "<<ic1 <<" "<<ic2 <<" "<<ic3<<flush;


//save track
                            if (chiq>-0.001&&chiq<chicut)
                            {
                                ntrk++;
                                // save traccia

                                Traccia[ntrk]=TrkNow;
                                Traccia[ntrk].id=ntrk;
                                Traccia[ntrk].camere=nchamb;
                                Traccia[ntrk].single=nsingle;
                                Traccia[ntrk].ic0=ic0;
                                Traccia[ntrk].ic1=ic1;
                                Traccia[ntrk].ic2=ic2;
                                Traccia[ntrk].ic3=ic3;
                                if (ic0>-1) {
                                    c0flag[ic0]=-1;
                                    c0ok=0;
                                }
                                if (ic1>-1) {
                                    c1flag[ic1]=-1;
                                    c1ok=0;
                                }
                                if (ic2>-1) {
                                    c2flag[ic2]=-1;
                                    c2ok=0;
                                }
                                if (ic3>-1) {
                                    c3flag[ic3]=-1;
                                    c3ok=0;
                                }
                            }
                        }
                    }
                }
            }
        }
    }


    charge_all=0;

    for (int i=0;i<ntrk+1;i++)
    {
        charge_all+=Traccia[i].charge;

        //  Traccia[i].Print_Trk();

    }

    Track2Root();


}

//=====================================================================
//================================F I T T I N G =======================

double TrgmpStraw::Fitting(short int mask_cam)   //,HitClus *//Hitp)
{

//    Gout<<"\n\n Ev "<< Eventi_Fatti<<" Fitting  mask " << mask;

    int uno=0;
    double z_zmg[4];
    double pcam=0;
    unsigned short mskc=1; 
   
    gmatrix F(8,5),Y(8,1);


    for (int icm=0;icm<4;icm++)
    {
        if ( (mask_cam&mskc)>0)  {  // la camera  esiste...
            pcam++;
            int idc=icm*2;
            double dzmg=Hitp[icm].Xcls.z-z_centro_mag;
            z_zmg[icm]=dzmg;
            int piani=Hitp[icm].planes;

            // code to fill F and Y  when there is more than 1 plan
            if (piani>1) {
                dzmg=Hitp[icm].Xcls.z-z_centro_mag;
                Y(idc,0)  =Hitp[icm].Xcls.x;
                Y(idc+1,0)=Hitp[icm].Xcls.y;
                F(idc,2)=F(idc+1,0)=1.;

                if (icm<2) {
                    F(idc,3)=F(idc+1,1)=dzmg;
                }
                else {
                    F(idc,4)=F(idc+1,1)=dzmg;

                }
            }

            // code to fill F and Y when there is only one plan
            else
            {
                short int mask_plan=Hitp[icm].mask;
                double cosa=1.,sena=0.;
                uno++;


                if (mask_plan&1) {
                    cosa=SINCOS45;
                    sena=SINCOS45;
                }
                else if (mask_plan&2) {
                    cosa=SINCOS45;
                    sena=-SINCOS45;
                }
                else if (mask_plan&4) {
                    cosa=1.;
                    sena=0.;
                }
                else if (mask_plan&8) {
                    cosa=0.;
                    sena=-1.;
                }
                F(idc,0)=sena;
                F(idc,1)=dzmg*sena;
                F(idc,2)=cosa;
                if (icm<2)F(idc,3)=dzmg*cosa;
                else F(idc,4)=dzmg*cosa;
                Y(idc,0)  =Hitp[icm].Filo+Hitp[icm].Centro.x*SINCOS45 ;
            }
        }
        mskc*=2;

    }
// F and Y  matrix done

    if (pcam<3.0) return -1.0; // at least tree chambers...

//=============== calcolo dei Y-values

    gmatrix M(5,5), Noto(5,1), X(5,1);

// debug =====

    if (Eventi_Fatti==errordebug)
    {
        Gout<<"\n Fitting Ev "<<Eventi_Fatti<<" mask_cam "<<mask_cam
        <<" msk "<< Hitp[0].mask  <<" "<< Hitp[1].mask  <<" "<< Hitp[2].mask  <<" "<< Hitp[3].mask;
        F.Print("F");
        Y.Print("Y");
    }
// end debug =====

    M=(~F)*F;
    M.Det();
    if (Eventi_Fatti==errordebug) M.Print("M");// debug =====

    if (fabs(M.detrm)<100000.) {
        Gout<<"\n Fitting mtx ERROR Ev Deteminante nullo!!!! "<<Eventi_Fatti<<" mask_cam "<<mask_cam
        <<" msk "<< Hitp[0].mask  <<" "<< Hitp[1].mask  <<" "<< Hitp[2].mask  <<" "<< Hitp[3].mask;
        return -50;
    }

    Noto=(~F)*Y;
    X=  (!M)*Noto;

//====== get the X-parameters by name
    double ymg,tgyz,xmg,tgxz_before,tgxz_after;
    ymg   =X(0,0);  // at magnet
    xmg   =X(2,0);  // at magnet
    tgyz  =X(1,0);
    tgxz_before=X(3,0);
    tgxz_after =X(4,0);
//
    double chiq=0;
    gmatrix DY(8,1),sumq(1,1);

    DY=Y-(F*X);
    sumq= (~DY)*DY;
    chiq=sumq(0,0);

// debug =====
    if (Eventi_Fatti==errordebug) {
        Y.Print("Y");
        DY.Print("DY");
        Gout<<"\n Ev "<<Eventi_Fatti<<" chiq "<<chiq;
    }
// end debug =====


    if (pcam==4.0&&chiq> 1.) return -2.;
    if (pcam==3.0&&chiq> 0.7) return -3.;

    // test di convergenza in y verso la zona di decadimento
    double zby_y= z_centro_mag-ymg/tgyz;
    if (zby_y<9000.||zby_y>18700.)return -10;

// test di convergenza in x verso la zona di decadimento

    double xb=0.0,zb=10000.;   // il beam passa per x =0 a z=10000
    double tgxb=-0.00124;       // inclinazione del beam (nota negativo)
    double zby_x= ((xmg-xb)-tgxz_before*z_centro_mag +tgxb*zb)/(tgxb-tgxz_before);
//    if (uno>0)Gout<<"\n Zby_x "<<zby_x<<" chiq "<<chiqy;
    if (zby_x<9000.||zby_x>18700.)return -11;


// ==================== calcolo del momento ===========
    gvet Versor_in(tgxz_before,tgyz,1.);
    Versor_in=Versor_in.NVerso();
    double cdip=sqrt(Versor_in.x*Versor_in.x+Versor_in.z*Versor_in.z);
    double kik=Mnp33->Get_Kik();
    double pinc= (kik/ (atan(tgxz_after)-atan(tgxz_before )) )/cdip;
    double charge=1.;
    if (pinc<0.0)charge=-1.;


// Gout<<"\nEv fit "<<Eventi_Fatti<<" parms "<<ymg<<" "<<tgyz<<" "<<xmg<<" " <<tgxz_before<<" " <<tgxz_after;

// define ============     TrkNow    ====================
    pinc=abs(pinc);

    TrkNow.mask=mask_cam;

    TrkNow.msk0= Hitp[0].mask;
    TrkNow.msk1= Hitp[1].mask;
    TrkNow.msk2= Hitp[2].mask;
    TrkNow.msk3= Hitp[3].mask;
    TrkNow.zby_y=zby_y;
    TrkNow.zby_x=zby_x;
    TrkNow.chiq_fit= chiq;
    TrkNow.Xmg.setv(xmg,ymg,z_centro_mag);
    TrkNow.P_in=Versor_in*pinc;
    gvet Versor_out(tgxz_after,tgyz,1.);
    Versor_out=Versor_out.NVerso();
    TrkNow.P_out=Versor_out*pinc;
    TrkNow.charge=charge;
    
// estrapolazione traccia all'odoscopio X
    double dkam= (hodox->Get_Cface().z-z_centro_mag)/Versor_out.z;
    TrkNow.Xhodx=TrkNow.Xmg+dkam*Versor_out;
    
  // estrapolazione traccia al calorimetro
    double dlkr= (lkry->Get_Cface().z-z_centro_mag)/Versor_out.z;
    TrkNow.Xlkr=TrkNow.Xmg+dlkr*Versor_out;  
    
    
    
    

    return chiq;  // return chiq;
}
//====================================================================================

int TrgmpStraw::Verify(short int mask_cam,int single)

{

    //test the number of single plans accepted
    // for 4 active chambers  we accept up to two different single planes but both  two  == to Y (8)
    // for 3 active chambers we accept up to two different single planes but not == to Y (8)

    unsigned short mskc=1;
    unsigned short msklast=0;
    unsigned short mskplan=0;
    TrkNow.used_plans=0;
    int singlepiani=0;
    for (int icm=0;icm<4;icm++)
    {
        if ( (mask_cam&mskc)>0)  {  // la camera esiste...
            short int mskpl= Hitp[icm].mask;
            if (single==1&&mskpl==8) return -1; // the y plan is not usefull in only trhee plans
            if (mskpl==1||mskpl==2||mskpl==4||mskpl==8)
            {
                if (msklast==8 && mskpl==8) return -1;
                singlepiani++;
                msklast=mskpl;
                mskplan=(mskplan|mskpl);
            }
            if (singlepiani>single) return -1;
        }
        mskc*=2;
    }
    TrkNow.used_plans=mskplan;

//    if (single==0&&msklast!=0) return -1;

    return 1;

}



//=================
void Track::Print_Trk()
{
    Gout<<"\nEv "<<Eventi_Fatti<<" id "<<id<<" ch "<<charge<<" msk "<<mask<<" cam "<<camere<<" SINGL "<<single
    <<" msk: "<<msk0<<" "<<msk1<<" "<<msk2<<" "<<msk3
    <<" chiq "<<chiq_fit;
    Xmg.print("Xmg");
    Gout<<" Ic.. "<<ic0<<" "<<ic1 <<" "<<ic2 <<" "<<ic3 ;
    P_in.print("P_in");
    P_out.print("P_out");
    Gout<<" Kik "<<(P_out-P_in).px;

}
// ================

void TrgmpStraw::Track2Root()
{
//int ntrk,mask[MAXTRACK],camere[MAXTRACK],ic0[MAXTRACK],ic1[MAXTRACK],ic2[MAXTRACK],ic3[MAXTRACK];
//  float chiq[MAXTRACK], charge[MAXTRACK], xmg[MAXTRACK], ymg[MAXTRACK], px[MAXTRACK],py[MAXTRACK],pz[MAXTRACK];
    ntrk_done=ntrk+1;

    for (int i=0;i<ntrk_done;i++)
    {
        Trk2Root.id[i]=  Traccia[i].id;
        Trk2Root.mask[i]=  Traccia[i].mask;
        Trk2Root.camere[i]=Traccia[i].camere;
        Trk2Root.single[i]=Traccia[i].single;
        Trk2Root.usepl[i]=Traccia[i].used_plans;
        Trk2Root.ic0[i]=   Traccia[i].ic0;
        Trk2Root.ic1[i]=   Traccia[i].ic1;
        Trk2Root.ic2[i]=   Traccia[i].ic2;
        Trk2Root.ic3[i]=   Traccia[i].ic3;
        Trk2Root.zbyy[i]=  Traccia[i].zby_y;
        Trk2Root.zbyx[i]=  Traccia[i].zby_x;
        Trk2Root.chiq[i]=  Traccia[i].chiq_fit;
        Trk2Root.charge[i]=Traccia[i].charge;
        Trk2Root.xmg[i]=   Traccia[i].Xmg.x;
        Trk2Root.ymg[i]=   Traccia[i].Xmg.y;
        Trk2Root.p[i]=     Traccia[i].P_in.mom;
        Trk2Root.py[i]=    Traccia[i].P_in.py;
        Trk2Root.pxin[i]=  Traccia[i].P_in.px;
        Trk2Root.pzin[i]=  Traccia[i].P_in.pz;
        Trk2Root.pxout[i]= Traccia[i].P_out.px;
        Trk2Root.pzout[i]= Traccia[i].P_out.pz;
        Trk2Root.xhodx[i]=   Traccia[i].Xhodx.x;
        Trk2Root.yhodx[i]=   Traccia[i].Xhodx.y;
        Trk2Root.xlkr[i]=   Traccia[i].Xlkr.x;
        Trk2Root.ylkr[i]=   Traccia[i].Xlkr.y;       
        
        

    }


}