Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

ConicBaffle - source file

// $Id: ConicBaffle.cc,v 1.5 2005/10/21 13:33:39 moreggia Exp $
// Author: Alessandro Thea   2005/04/08

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: ConicBaffle                                                           *
 *  Package: <packagename>                                                   *
 *  Coordinator: <coordinator>                                               *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// ConicBaffle
//
// <extensive class description>
//
//   Config file parameters
//   ======================
//
//   fTopRadius [mm]: radius of the optics the baffle is attached to.
//   -Valid options: <available options>
//
//   fAlpha [deg]: baffle's field of view.
//   -Valid options: 0 to 60 degrees.
//

#include "ConicBaffle.hh"
#include "EConst.hh"
#include "Etypes.hh"

using namespace TMath;
using namespace EConst;

ClassImp(ConicBaffle)

//_____________________________________________________________________________
 ConicBaffle::ConicBaffle() {
    //
    // Constructor
    //

    fRTop  = Conf()->GetNum("ConicBaffle.fTopRadius");
    fAlpha = Conf()->GetNum("ConicBaffle.fAlpha");

    if ( fAlpha < 0 || fAlpha > 60 )
        FatalError("fAlpha must be larger than 0 and smaller than 60 deg.");
    fAlpha *= DegToRad();

    fDZdown = Conf()->GetNum("ConicBaffle.DZdown"); 
    fPos    = EVector(0,0,Conf()->GetNum("ConicBaffle.fPos.Z"));

    fR = fRTop + fDZdown*Tan(fAlpha); 

    if ( Abs(TopRadius()-BottomRadius()) < kTolerance ) {
        fCotAlpha = 0;
        fVPosZ = 0;
    } else {
        fVPosZ = -fDZdown+(fDZdown/(BottomRadius()-TopRadius()))*BottomRadius();
        fCotAlpha = BottomRadius()/(-fDZdown-fVPosZ);
    }
}

//_____________________________________________________________________________
 ConicBaffle::~ConicBaffle() {
    //
    // Destructor
    //
}

//______________________________________________________________________________
 Double_t ConicBaffle::NextInteraction( Photon* ph ) const{
    //
    //  finds the next interaction of the photon inside the cylinder
    //
    
    Bool_t fgDebug = kFALSE;

    Double_t dist[4] = {-kHuge,-kHuge,-kHuge,-kHuge};
    // array of intersections, 0, 1 basis, 3, 4 lateral surface
    
    if (fgDebug) Printf("pos(%.8f,%.8f,%.8f)   dir(%.8f,%.8f,%.8f)",
            ph->pos[X],ph->pos[Y],ph->pos[Z],ph->dir[X],ph->dir[Y],ph->dir[Z]);
    
    // improve read
    TVector3& dir = ph->dir;
    TVector3& pos = ph->pos;

    // int with the basis
    if ( Abs(dir[Z]) > kTolerance ) {
        dist[0] = -pos[Z]/dir[Z]; // top basis
        dist[1] = (-fDZdown-pos[Z])/dir[Z]; // bottom basis
    }

    Double_t a, b, c;
    if ( fCotAlpha == 0 ) {
        a = dir[X]*dir[X]+dir[Y]*dir[Y];
        b = 2*(dir[X]*pos[X]+dir[Y]*pos[Y]);
        c = pos[X]*pos[X]+pos[Y]*pos[Y]-TopRadius()*TopRadius();
    } else {
        // vertex
        a = dir.Perp2()-fCotAlpha*fCotAlpha*dir[Z]*dir[Z];
        b = 2*(dir[X]*pos[X]+dir[Y]*pos[Y]-fCotAlpha*fCotAlpha*dir[Z]*(pos[Z]-fVPosZ));
        c = pos[X]*pos[X]+pos[Y]*pos[Y]-fCotAlpha*fCotAlpha*(pos[Z]-fVPosZ)*(pos[Z]-fVPosZ);
    }


    Double_t delta = (b*b)-4*a*c;

    if (fgDebug) Printf("a=%.7e, b=%.7e, c=%.7e, delta=%.7e",a,b,c,delta);
    if ( Abs(a) > kTolerance ) 
        if ( Abs(delta) < kTolerance ) {
            dist[2] = -b/(2*a);
        } else if ( delta >= kTolerance ) {
            dist[2] = (-b-Sqrt(delta))/(2*a);
            dist[3] = (-b+Sqrt(delta))/(2*a);
        }


    Double_t dt = MaxElement(4,dist); 
    for(Int_t i(0); i<4; i++) 
        if ( dist[i] > kTolerance && dist[i] < dt ) dt = dist[i];

    if (fgDebug) Printf("d(%e,%e,%e,%e)   and   dt=%.3f",dist[0],dist[1],dist[2],dist[3],dt); 
    return dt;
}

//______________________________________________________________________________
 Photon* ConicBaffle::Transport( Photon* ph ) const {
    //
    // Transport
    //

    // to local
    ph->pos -= fPos;

    Double_t dt = NextInteraction( ph );

    ph->pos += ph->dir*dt;
    ph->time+= dt/Clight();
    ph->pos += fPos;
    
    if ( Abs(ph->pos[Z] - Bottom()) < kTolerance ||
         Abs(ph->pos[Z] - Top()) < kTolerance )
        return ph;
    else
        return 0;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005