flyopunta/src/src_ana/analisi.cpp

/***************************************************************************
                          analisi.cpp  -  description
                             -------------------
    begin                : Sat Jan 18 2003
    copyright            : (C) 2003 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/

//    Procedura per l'analisi dei dati
// Tirrenia 5.09.00       by GmP

#define EPSR     1.0E-200
#include "flyoh.h"
#include "flyoana.h"

//---------------------------------------------
//
//==================  A n a l i s i  ============
Analisi::Analisi()
{
        static Analisi *last=0;

// links...
        if ( last == 0 ) Primofit = this;
        else             last->next = this;
        up = last;
        last = this;
        this->next = 0;

        count_call=count_decay=getta_gg=getta_en=getta_cog=getta_tau=getta_gf=0;
        getta_trk=getta_ch=getta_fit=count_wnt=getta_char=getta_kabes=0;
        for (int i=0;i<30;i++) counter[i]=0;
        ana_tempo=ana_dtim=0.0;
        
// get the  chamber pointers.....
// occorrono i nome definiti nella file *.epc
//......................
        titol=" Analisi title Not defined! ";

 /*   da definire a run time altrove??????????????????????????????
 // Avo parameters
  Particella *Pr_avo=SelRea->avo;
   msAvo= Pr_avo->Get_Massa(); 
   msqAvo=msAvo*msAvo;
   ctauAvo=Pr_avo->Get_Ctau();
  
//      msAvo=.49767;          // default:  Kl  mass

        msqAvo=msAvo*msAvo;
        reckik=0.0;
        if ( Mag48!=0 )
        {
                reckik=Mag48->Get_Kik();
                Gout<<"\n Analisys momentum kik="<<reckik;
        }
*/



}
Analisi::~Analisi() {  }

void Analisi::Write_titol()
{
        Gout<<"\n Ev "<< evento_.Gen.Event<<"  Fit: "<<nome<<" => Rea: "<<titol<<"  <=";
}
//----------------------------------------------------
void Analisi::Set_tempo_parms()
{
        ana_dtim=evento_.Gen.tempo-ana_tempo;
        ana_tempo=evento_.Gen.tempo;    
}


void Analisi::print_scale()
{

        Gout<<"\n\n ----->  "<<titol<<"   Analysis Summary. <-------\n";
        Gout<<"\n Fit entries          "<<count_call;
        Gout<<"\n No kabes             "<<getta_kabes;
        Gout<<"\n Too few tracks       "<<getta_ch;
        Gout<<"\n Killed by energy     "<<getta_en;
        Gout<<"\n Tau cut              "<<getta_tau;
        Gout<<"\n No good vertex       "<<getta_fit;
        Gout<<"\n Good events          "<<count_wnt;
        Gout<<"\n\n ======>           d o n e     <============ ";
}


double Analisi::Vertice2(gvet &Xvert,gvet P1,gvet V1,gvet P2,gvet V2)
{
        // calcolo del vertice o meglio si calcola il punto più
        // vicino a due rette e si  ritorna come chiq il quadrato delle distanze
        // del punto(vertice) delle rette
        //   si potrebbe facilmente generalizzare a n rette.
static gmatrix Unit(3,1.0);
static gvet Vnoto;
static gmatrix Pro1(3,3),Pro2(3,3),Pros(3,3);

        
    V1>>Pro1;    //matrice Pro1= |V1><V1|
    Pro1=Unit-Pro1;
    V2>>Pro2;    //matrice Pro2= |V2><V2|
    Pro2=Unit-Pro2;
    Pros=Pro1+Pro2;
    Vnoto= Pro1%P1 +Pro2%P2;
    Pros.Det();
    Xvert= (!Pros)%Vnoto;
//chiq
    gvet d1= Pro1%(Xvert-P1);
    gvet d2= Pro2%(Xvert-P2);
    double chiq= d1%d1+d2%d2;
    return chiq;

}

//  calcolo del vertice con due o più tracce

Vtxstruct::Vtxstruct(int n)
  {
   Ntr=n;
   P=new gvet[Ntr];
   V=new gvet[Ntr];
   Proj= new gmatrix[Ntr];
   for(int i=0;i<Ntr;i++)Proj[i].Set_Size(3,3);
   gmatrix U(3,1.0);
      
  }

void Vtxstruct::Vertice()
{
// algebra matriciale 

static gmatrix Unita(3,1.0);
gmatrix Prjtot(3,3);  
gvet noto;  
for(int i=0;i<Ntr;i++)
{
V[i]>>Proj[i];//matrice Proj= |V><V|
Proj[i]=Unita-Proj[i];
Prjtot+=Proj[i];
noto+=Proj[i]%P[i]; 
}
  
Prjtot.Det();
Vtx=(!Prjtot)%noto;

chiq=0.0;
for(int i=0;i<Ntr;i++)
{
  gvet dist= Proj[i]%(Vtx-P[i]);
  dist.Norma();
   chiq+= dist.normaq;
}
  
 
  
}