flyopunta/src/src_uti/illumina.cpp

/***************************************************************************
 *   Copyright (C) 2009 by giuseppe Pierazzini                                      *
 *   giuseppe@pierazzini.it                                                               *
 *                                                                         *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.                  *
 *                                                                         *
 *  Dipartimento di Fisica  E.Fermi / INFN Pisa Italy                      *
 ***************************************************************************/

#include "illumina.h"
#include "flyoh.h"

//using namespace std;

Illumina::Illumina()
{
    npr=0;
    illumina_buffer_dim=100;  //attention: ok with cuts... without cuts better to keep 300.
    Make_illumina_buffer ( illumina_buffer_dim );
    Gout<<"\n Illumination required: DONE  of "<<illumina_buffer_dim<<"  dimension";
}


Illumina::~Illumina()
{
}
//==================R a d x =============================
void Illumina::Make_illumina_buffer ( int rdim )
{
// make space
    idm    = new int[rdim];
    id     = new int[rdim];
    ifato  = new int[rdim];
    devstry= new unsigned long long[rdim];
    xg     = new float[rdim];
    yg     = new float[rdim];
    zg     = new float[rdim];
    tg     = new float[rdim];
    xf     = new float[rdim];
    yf     = new float[rdim];
    zf     = new float[rdim];
    pxg    = new float[rdim];
    pyg    = new float[rdim];
    pzg    = new float[rdim];
    eg     = new float[rdim];
    pxf    = new float[rdim];
    pyf    = new float[rdim];
    pzf    = new float[rdim];
    ef     = new float[rdim];
    mass   = new float[rdim];
    theta  = new float[rdim];
    thetastar  = new float[rdim];
    pxsg    = new float[rdim];
    pysg    = new float[rdim];
    pzsg    = new float[rdim];
    esg     = new float[rdim];
}
//===================== I l l u m i n a =========================
void Illumina::Store_Illumina()
{
    //  generated data...info for debug  activated by illumina
    // do not  activate if  pipeline is on
    Particella *pr;
    const char *ft;
    static qvet Stesso;

    // calcoli per gas
    pr=SelRea->avo->next;
    qvet Ptot;
    while (pr!=0&&pr->Get_Id()==0) // calcolo del quadriimpulso delle particelle generate da Fluka
    {
        Ptot+=pr->Get_Gp();
        pr=pr->next;

    }
    Ptot.Invarq();
    Ptot.putq(pxtgas,pytgas,pztgas,Etgas);
    mqgas=Ptot.mq;
    Ptot.x=-1.*Ptot.x;
    Ptot.y=-1.*Ptot.y;
    Ptot.z=-1.*Ptot.z;
    qvet Ppart=SelRea->avo->Get_Gp();
    Ppart.Lburst(Ptot);

    pr=SelRea->avo;
    npr=0;
    while ( pr!=0 )
    {
        if (dectubo==0) {               //luce per TUTTE le particelle generate
            if ( ( ft=pr->Get_Fato() ) == Undef //particelle non inizializzate
                    ||  pr->Get_Ifato() == 0     // particelle non inizializzate (superfluo)
               ) {
                pr=pr->next;
                continue;
            }
        }
        else { //evento fluka: luce per le particelle VISIBILI
            // to be adjusted depending on the study
            if ( ( ft=pr->Get_Fato() ) == Undef  //particelle non inizializzate
                    ||  pr->Get_Ifato() == 0     // particelle non inizializzate (superfluo)
                    || pr->Get_See() == 0        //particelle invisibili (es. nu,pi0..)
                    || (pr->Get_DevStory()==0 &&  pr->Get_Ifato()==-2) //particelle che decadono prima di essere rivelabili
                    || pr->Get_Gp().e == 0 // particelle di energia nulla! (circa 300 fotoni su 10e6 eventi gas)
                    || (pr->X-pr->Get_Gx()).Norma() == 0 // particelle "virtuali"
                    // particelle generate fuori dall'apparato sperimentale
                    || pr->Get_Gx().z < 10200 || pr->Get_Gx().z > 25000 //taglio su z
                    //taglio sul tubo
                    || (pr->Get_Gx().XYNorma() > 55 && pr->Get_Gx().z > 10200 && pr->Get_Gx().z < 17000)
                    //taglio sui detectors downstream
                    || (pr->Get_Gx().XYNorma() > 80 && pr->Get_Gx().z > 17000 && pr->Get_Gx().z < 20000)
                    || (pr->Get_Gx().XYNorma() > 115 && pr->Get_Gx().z > 20000 && pr->Get_Gx().z < 25000)
               )
            {
                pr=pr->next;
                continue;
            }
        }
        idm[npr]= pr->Get_Idm();
        id[npr]=  pr->Get_Id();
        ifato[npr]= pr->Get_Ifato();
        devstry[npr] = pr->Get_DevStory();
//    Gout<<"\n ift "<< ifato[npr];
        tg[npr]=evento_.Gen.tempo;
        //
        pr->Get_Gx().putv( xg[npr],yg[npr],zg[npr] );
        pr->X.putv( xf[npr],yf[npr],zf[npr] );
        pr->Get_Gp().putq( pxg[npr],pyg[npr],pzg[npr],eg[npr] );
        pr->P.putq( pxf[npr],pyf[npr],pzf[npr],ef[npr] );

        mass[npr]=float ( pr->Get_Massa() );

        //angolo per di uscita -> gas_eventi
//         if (pr->Get_Gpz()>0.)
//             theta[npr] = asin(pr->Get_Gp().XYNorma()/pr->Get_Gpn());
//
//
//         else if (pr->Get_Gpz()<0.)
//             theta[npr] = GmpConst::PiGreco - asin(pr->Get_Gp().XYNorma()/pr->Get_Gpn());
//         else
//             theta[npr] = asin(1.0);

        theta[npr] =acos(pr->Get_Gpz()/pr->Get_Gpn());
//calcolo del cost* ( rispetto a z)  qui per tutte le particelle generate
        Ppart=pr->Get_Gp();
        Ppart.Lburst(Stesso);
        thetastar[npr]=acos(Ppart.z/Ppart.mom);
        Ppart.putq( pxsg[npr],pysg[npr],pzsg[npr],esg[npr] );



        if ( npr>illumina_buffer_dim-2 ) {
            Gout<<"\nToo many illumina particles ... cut to  "<<  illumina_buffer_dim ;
            break;
        }
        npr++;


        pr=pr->next;
    }

}