OFTAdapter.C

Go to the documentation of this file.

#ifndef GW_GAMMATRACKER_H
#include <GammaTracker.h>
#endif

using namespace Gw;

class OFT : public GammaTracker
{
public:
        OFT(const char *name = "OFT", const char *title = "Orsay's forward tracking program");
        virtual ~OFT() ;

        virtual Int_t DoTracking();

        ClassDef(OFT,0); // Orsay Forward Tracker
};

// real implementation

#ifdef __cplusplus 
extern "C" {
#endif
#include <tracking_define.h>
#include <tracking_data.h>
#include <tracking_data_manip.h>
#include <tracking_utilitaires.h>
#include <tracking_cluster.h>
#include <tracking_events.h>
#ifdef __cplusplus 
}
#endif


namespace {

Double_t r[intmax][intmax], r_ge[intmax][intmax], erfcu[1000], etot[intmax], probtot[intmax], angn[intmax];
Int_t sn[intmax][kmax], numn[intmax], order[intmax], nombre[intmax], flagu[intmax], nbtot[intmax];
Double_t angtheta[intmax], angphi[intmax];
Double_t et[intmax];
Int_t numinter[intmax], numintern[intmax], interaction[intmax][kmax];

Double_t ddu; Int_t nintu;

} // anonymous namespace to protect globals


OFT::OFT(const char *name,const char *title) : GammaTracker(name,title)
{ 

        // Init some constants 
        Init_Tracking_Lib(); 

  /**********************************************************/
  /* mapping random number vs uncertainty in units of sigma */
  /**********************************************************/

        // TODO - should be set through a member function?


        nintu = 1000;   /* initialisation for random number vs uncertainty procedure */
        ddu = 0.005;    /* step in random number vs uncertainty procedure */

        maprandom(erfcu, ddu, nintu);
        
        SetDistFactor(10.0); SetEnerFactor(1000.0); // units are centimeters and MeV 
}

OFT::~OFT() {;}


Int_t OFT::DoTracking()
{
        Int_t n, i, j, itot, itotcurr, mult;

        // symbolic link between local notation and global notation
        Double_t *e = InputE; Double_t *x = InputX, *y = InputY, *z = InputZ;
        Double_t *thetafirst = OutputTheta1, *phifirst = OutputPhi1, *thetasecond = OutputTheta2, *phisecond = OutputPhi2;

        if ( InputN == 0 ) { OutputN = 0; return OutputN; }
        
        Int_t nb_int = InputN;

        mult = 0;
   /****************************************************/
  /******* pack points with r < d_res  ****************/
  /****************************************************/
        if ( gDebug ) printf("pack points\n");
        nb_int = packpoints(nb_int, x, y, z, e, numinter, order);

/*
        Double_t packenergy, totenergy = 1.0; 

        packenergy = 0;
        for (j = 0; j < nb_int; j++)
            packenergy = packenergy + e[j];
        if (fabs(packenergy - totenergy) > 1e-4)
            printf("pb with energy after packing %f %f %d %d \n", packenergy, totenergy, nb_int,InputN); 
        else 
            printf("NO pb with energy after packing %f %f %d %d \n", packenergy, totenergy, nb_int,InputN); */
        
  /*******************************************************/
  /********* smear positions *****************************/
  /*******************************************************/
        if ( gDebug ) printf("smear positions\n");
        smearpoints(nb_int, x, y, z, erfcu, ddu, nintu);


  /*****************************************************/
  /*********** Apply Energy threshold ******************/
  /*****************************************************/
        if ( gDebug ) printf("Apply Energy threshold\n");
        nb_int = energy_thresh(nb_int, x, y, z, e, numinter, order);

        /*
        totenergy = 0;
        for (i = 0; i < InputN; i++) {
            totenergy = totenergy + e[i];
        }
        if (InputN == 1 && totenergy >= emini)
            nsingint++;
        nbinttot = nbinttot + InputN;
        if (totenergy > emini) {
            nbinter = nbinter + InputN;
        } */

  /*****************************************************/
  /******** compute distance source-int and ************/
  /*************theta and phi for each int *************/
  /*****************************************************/
        for (i = 0; i < nb_int; i++) {
            r[i][i] = sqrt(SQ(x[i] - xsource) + SQ(y[i] - ysource) + SQ(z[i] - zsource));
            angtheta[i] = acos((z[i] - zsource) / r[i][i]);
            angphi[i] = atan2((y[i] - ysource), (x[i] - xsource));
            if (angphi[i] < 0)
                angphi[i] = 2 * PI + angphi[i];
        }

  /*****************************************************/
  /********* sort according to increasing theta*********/
  /*****************************************************/
        if ( gDebug ) printf("sort according to increasing theta\n");
        for (i = 0; i < nb_int; i++) {
            for (j = i + 1; j < nb_int; j++) {
                if (angtheta[j] < angtheta[i]) {
                    swap(e, i, j);
                    swap(x, i, j);
                    swap(y, i, j);
                    swap(z, i, j);
                    swapi(numinter, i, j);
                    swapi(order, i, j);
                    swap(angtheta, i, j);
                    swap(angphi, i, j);
                }
            }
        }

  /*****************************************************************/
  /******** computing all distances between interactions   *********/
  /*****************************************************************/
        if ( gDebug ) printf("computing all distances between interactions\n");
        distances(r, r_ge, nb_int, x, y, z);

        /* computes r and r_ge (real distances and effective distances in Ge)
           // given position vectors x,y and z and the number of interaction points InputN
           // NB: distances in Ge are calculated assuming a 4pi ge sphere */

  /****************************************************/
  /*************** find clusters ***********************/
  /****************************************************/
        if ( gDebug ) printf("find clusters\n");
        n = cluster_search(numn, angn, et, sn, nb_int, e, angtheta, angphi);

        /* numn = number of interactions in cluster
           // angn = angle of cluster 
           // et =  total energy of cluster
           // sn[1][2] = identity of 2nd point of cluster 1
           // e = energy vector of interaction points
           // angtheta and angphi = theta and phi vector of interaction points
           // routine returns number of clusters found */

  /*******************************************************/
  /******** compute figure of merit of clusters **********/
  /*******************************************************/
        if ( gDebug ) printf("compute figure of merit of clusters\n");
        cluster_evaluation(n, et, numn, sn, probtot, interaction, e, x, y, z, r, r_ge, flagu);

        /* n = number of clusters
           // probtot = figure of merit of cluster
           // interaction = ordering of interactions in cluster which gives best 
           // figure of merit (set to minprobtrack for singles temporarily)
           // x,y,z,r and e vectors of x,y,z and energy of interaction points
           // r and r_ge = real and effective distances in Ge
           // flagu - is set to 0 for all clusters */

  /***********************************************************************/
  /************ sort clusters according to figure of merit ****************/
  /************* flag worst clusters with same interactions **************/
  /***********************************************************************/
        if ( gDebug ) printf("sort clusters\n");
        cluster_sort_flag(n, et, numn, sn, angn, probtot, interaction, flagu);

        /* flagu is set to 1 if a cluster with a worst figure of merit than another one shares at least one common interaction point */

  /**********************************/
  /* test of direct photopic events */
  /**********************************/
        if ( gDebug ) printf("test of direct photopic events\n");
//      if (nodirect == 0) {    /* treat single interaction points */
            single_interaction(nb_int, r, r_ge, n, et, numn, sn, probtot, flagu);
//      }

        /* routines computes figure of merit for single interaction points */

  /**************************************************************/
  /*********************** cluster validation ******************/
  /*************************************************************/
        if ( gDebug ) printf("cluster validation\n");
        itotcurr = itot;
        mult = cluster_validation(etot, nbtot, n, probtot, numn, interaction, et, flagu, /*e, */ x, y, z, r, thetafirst, phifirst, thetasecond, phisecond);

        /*routine returns number of tracked gamma rays
           // etot - energy of tracked gamma rays
           // nbtot - number of interaction points making up the gamma
           // thetafirst and phifirst are the incident angle direction
           // thetasecond and phisecond are he scattered direction (wrt the incident direction) */

        // set output
        OutputN = mult;  ::memcpy(OutputE,etot,sizeof(Double_t)*OutputN); return OutputN;
}