Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

KOpticalSystem - source file

// ESAF : Euso Simulation and Analysis Framework
// class KOpticalSystem
// $Id: KOpticalSystem.cc,v 1.37 2005/10/21 13:33:39 moreggia Exp $
// J.Watts created , Sep 21 2003
//

#include "Ktrace_optF1v1.hh"
#include "KEUSO_optF1v1.hh"
#include "KOpticalSystem.hh"
#include "Config.hh"
#include "EsafRandom.hh"
#include "EEvent.hh"

using namespace sou;

ClassImp(KOpticalSystem)

Int_t Kread_para(char*);
Int_t KReadSD(char*);
int Ktrace(double*, double*);

//______________________________________________________________________________
 KOpticalSystem::KOpticalSystem(): OpticalSystem() {
    //
    // Constructor
    //

    fDZdown=0*mm; 
    fDZup=2002*mm;

    // take the size of the system from Ktrace
    fR=R_WALL;
	fPos=EVector(0,0,Conf()->GetNum("KOpticalSystem.fPos.Z")*mm);

    fFirstLensDZ=169*mm;
    fSecondLensDZ=226*mm;

#ifdef DEBUG
/*    cout << "FirstLensTop() "<< FirstLensTop() << endl;
    cout << "FirstLensBottom() "<< FirstLensBottom() << endl;
    cout << "SecondLensTop() "<< SecondLensTop() << endl;
    cout << "SecondLensBottom() "<< SecondLensBottom() << endl;*/
#endif /* DEBUG */
    
    if(Kread_para("config/Optics/KOpticalSystem")<0){
       cout<<"reading surface data files failed\n"<<endl;
       exit(-2);
	}
   if(KReadSD("config/Optics/KOpticalSystem")<0){
       cout<<"reading diffractive data files failed\n"<<endl;
       exit(-2);
	}
}

//______________________________________________________________________________
 Photon *KOpticalSystem::Transport(Photon *p) const {
    double ph_in[8];
    double ph_out[8];
	
	EVector dummy=p->pos-fPos;
    ph_in[0]=dummy[X]/mm;
    ph_in[1]=dummy[Y]/mm;
    ph_in[2]=dummy[Z]/mm;

    ph_in[3]=p->dir[X];
    ph_in[4]=p->dir[Y];
    ph_in[5]=p->dir[Z];
    
    ph_in[6]=p->wl/nm;
    ph_in[7]=p->time/ns;


#ifdef DEBUG/*
    cout<<"KOpticalSystem"<<endl;
    cout<<"position = "<<p->pos<<endl;
    cout<<"direction = "<<p->dir<<endl;
    double theta=p->dir.Angle(EVector(0.,0.,1.));
    cout<<"theta = "<<theta/deg<<endl;
    cout<<"wl = "<<p->wl/nm<<endl;
    cout<<"time = "<<p->time/ns<<endl;*/
#endif /* DEBUG */
    int flag=Ktrace(ph_in, ph_out);
#ifdef DEBUG
    //cout << flag << endl;
#endif /* DEBUG */
	
    dummy[X]=ph_out[0]*mm;
    dummy[Y]=ph_out[1]*mm;
    dummy[Z]=ph_out[2]*mm;
	p->pos=dummy+fPos;
	
    p->dir[X]=ph_out[3];
    p->dir[Y]=ph_out[4];
    p->dir[Z]=ph_out[5];

    p->wl   = ph_out[6]*nm;
    p->time = ph_out[7]*ns;
    p->fate = flag;
    
    if ( flag != 0 )
        return 0;
    
    Double_t DL; 
    // sometime Ktrace returns photons out of the boundaries
    if (p->pos.Perp() > R_WALL || p->pos[Z] > Top() || p->pos[Z] < Bottom() ){
        // the photon is outside the cyl, trace it back to the borders
        EVector savedir = p->dir;
        p->dir = -p->dir;
        if ((DL = CylinderIntersection(p,kTRUE)) > 0){
            p->pos+=p->dir.Unit()*DL;
            p->time-=DL/TMath::C();
        }
        p->dir=savedir;
    } else if ((DL = CylinderIntersection(p,kTRUE)) > 0){
    // all the photon must be returned on the boundaries of the optics

        p->pos+=p->dir.Unit()*DL;
        p->time+=DL/TMath::C();
    }
    
    
#ifdef DEBUG
     //   cout<<"new position = "<<p->pos<<endl;
       // cout<<"new direction = "<<p->dir<<endl;
//        cout<<"new wl = "<<p->wl/nm<<endl;
  //      cout<<"new time = "<<p->time/ns<<endl;
#endif /* DEBUG */
        
    if (IsGoingOut(p))
        return p;
    else
        return 0;
}

//______________________________________________________________________________
 Bool_t KOpticalSystem::IsGoingOut( Photon *p ) const {
    // 
    // helper method to select outward going photons 
    // 

    Bool_t side, face;
    // check if photon is on the lat surface
    side = TMath::Abs((p->pos-fPos).Perp()-fR) < kTolerance  && 
        (p->pos-fPos).XYvector()*p->dir.XYvector() > 0 &&
        ((p->pos[Z] < FirstLensTop() && p->pos[Z] > FirstLensBottom()) ||
         (p->pos[Z] < SecondLensTop() && p->pos[Z] > SecondLensBottom()) );
    // or if it's on the top/bottom face
    face = (TMath::Abs(p->pos[Z]-FirstLensBottom()) < kTolerance && p->dir[Z]<0) ||
        (TMath::Abs(p->pos[Z]-SecondLensTop()) < kTolerance && p->dir[Z]>0); 
    return side || face;
}

//______________________________________________________________________________
 Bool_t KOpticalSystem::HitSide(Photon *p) const {
    //
    // Check if p hist the side of the Cylinder.
    //

    Double_t dt = CylinderIntersection( p, fR, SecondLensBottom(), FirstLensTop() );
    if (dt < 0 || (p->pos+p->dir*dt)[Z]-SecondLensBottom() < kTolerance 
               || (p->pos+p->dir*dt)[Z]-FirstLensTop() < kTolerance )
        return kFALSE;
    return kTRUE;
}
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