Subversion Repositories f9daq

#ifndef CVODNIK_H

#define CVODNIK_H

#include "TROOT.h"

#include "TVector3.h"

#include "TMath.h"

#include "TPolyLine3D.h"

#include "TRandom2.h"

#include "TDatime.h"

#include "TH1F.h"

#include "TH2F.h"

#include "TColor.h"

#define MARGIN 1.0e-6

#define DEGREE 180.0/3.14159265358979312

#define MAX_REFLECTIONS 19

using namespace TMath;

const double c_reflectivity = 0.96;

const int REFRACTION = 1;

const int REFLECTION = 2;

const int TOTAL = 3;

const int ABSORBED = -2;

const TVector3 CENTER(-2.0, 0.0, 0.0);

const int nch = 50; // number of "channels" used as number of bins in histograms

enum Fate {missed=0, hitExit=1, refracted=-1, enter=2, noreflection=-2, backreflected=-3, rays=3, absorbed=4};

int dbg = 0;

//=================================================================================

// zarek (premica)

class CRay

{

public:

  CRay() :

    r(TVector3(0,0,0)),

    k(TVector3(1,0,0)),

    p(TVector3(0,1,0)),

    color(kBlack)

  {};

  CRay(TVector3 r0, TVector3 k0) :

    r(r0),

    k(k0.Unit()),

    p(TVector3(0,1,0)),

    //p(TVector3(0,0,1)),

    color(kBlack)

  {};

  CRay(double x0, double y0, double z0, double l0, double m0, double n0) :

    r(TVector3(x0,y0,z0)),

    k(TVector3(l0,m0,n0).Unit()),

    p(TVector3(0,1,0)),

    //p(TVector3(0,0,1)),

    color(kBlack)

  {};

  void Set(TVector3 r0, TVector3 k0);

  //void Set(double x0, double y0, double z0, double l0, double m0, double n0);

  void SetColor(int c){color = c;};

  void setPolarization(TVector3 p0) {

    p = p0.Unit();

    if (p.Dot(k)>1e-3) printf("*** ERROR in CRay: E has component || with k\n");

  };

  inline CRay & operator = (const CRay &);

  TVector3 GetR() const {return r;};

  TVector3 GetK() const {return k;};

  TVector3 GetP() const {return p;};

  void Print();

  void Draw();

  void Draw(double x_from, double x_to);

  void DrawS(double x_from, double t);

  //r = position, k = unit wave vector, p = polarization

private:

  TVector3 r;

  TVector3 k;

  TVector3 p;

  int color;

};

//=================================================================================

// ravnina, definirana in omejena s 4 tockami na njej (stirikotni poligon)

// zanasam se na taksno definicijo: (normala kaze noter)

// r1----->r2

// ^        |

// |        |

// |        |

// |

// r0<-----r3

class CPlane4

{

public:

  CPlane4();

  CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);

  CPlane4(TVector3 *vr);

  ~CPlane4() { };

  void Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);

  void Set(TVector3 *vr) {Set(vr[0], vr[1], vr[2], vr[3]);};

  void FlipN(){n = -n;};

  int GetIntersection(TVector3 *vec, CRay ray);

  int IsInTri(TVector3 vec, TVector3 e1, TVector3 e2, double);

  int IsVectorIn(TVector3 vec);

  int TestIntersection(CRay in); // ray in & ??? plane

  int TestIntersection(TVector3 *vec, CRay in); // plane defined with normal vec and ray in

  TVector3 GetN() {return n;};

  void Print();

  void Draw(int color = 1, int width = 1);

private:

  TVector3 r[4];

  TVector3 n;

  double A, B, C, D;

  TVector3 edge[4]; // vektorji stranic

  double angle_r[4]; // koti ob posameznem vogalu

};

//=================================================================================

// ravnina - krog

class CPlaneR

{

public:

  CPlaneR(TVector3 c, TVector3 n0, double R0)

  {center = c; n = n0; _r = R0;};

  CPlaneR(){ n=TVector3(0,0,0); center=TVector3(0,0,0); _r=0; };

  ~CPlaneR() {};

  void Set(TVector3 c, TVector3 n0, double R0)

  {center = c; n = n0; _r = R0;};

  int TestIntersection(TVector3 *vec, CRay in);

  void Draw(int color = 1, int NN = 32);

private:

  TVector3 n, center;

  double _r;

};

//=================================================================================

// ravna opticna povrsina: refractor, zrcalo ali povrsina s totalnim odbojem

#define SURF_DUMMY 0

#define SURF_REFRA 1

#define SURF_REFLE 2

#define SURF_TOTAL 3

//#define SURF_TOTAL 1

#define SURF_IMPER 4

//#define SURF_DIFUS 5

class CSurface : public CPlane4

{

public:

  CSurface(int type0 = 0);

  CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4,

      double n10, double n20, double reflectivity);

  CSurface(int type0, TVector3 *vr,  double n10, double n20, double reflectivity);

  ~CSurface() {};

  CSurface& operator=(const CSurface& rhs) {

    type=rhs.type;

    n1=rhs.n1;

    n2=rhs.n2;

    n1_n2=rhs.n1_n2;

    reflection=rhs.reflection;

    cosTtotal=rhs.cosTtotal;

    fresnel=rhs.fresnel;

    return *this;

  }

  void SetV(TVector3 *vr){CPlane4::Set(vr);};

  void SetType(int type0){type = type0;};

  void SetIndex(double n10, double n20);

  void SetReflection(double reflectivity){reflection = reflectivity;};

  void SetFresnel(int f = 1) {fresnel = f;};

  void Set(int type0, TVector3 *vr, double n10, double n20, double reflectivity)

  {type = type0; CPlane4::Set(vr); SetIndex(n10, n20); reflection = reflectivity;};

  void Set(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)

  {type = type0; CPlane4::Set(r1, r2, r3, r4); SetIndex(n10, n20); reflection = reflectivity;};

  int PropagateRay(CRay in, CRay& out, TVector3 *intersection);

  double N1_N2(int sign) { return ((sign > 0) ? n1/n2 : n2/n1);};

  void Print(){printf("Type = %d\n", type); CPlane4::Print();};

private:

  int type; //0 = dummy; 1 = refractor; 2 = reflector; 3 = total reflection

  double n1, n2, n1_n2; //index of refraction, n1 @ +normal, n2 @ -normal

  double reflection; // odbojnost stranic

  TRandom2 rand; // za racunanje verjetnosti odboja od zrcala

  double cosTtotal; //cosinus mejnega kota totalnega odboja za dana n1 in n2

  int fresnel; // ali naj uposteva Fresnelove enacbe

};

//=================================================================================

//    +-------   n1                

//    |       -------              

//    |          n2  ----+         

//    |                  |         

//   0| M*SiPM           | SiPM    

//    |                  |         

//    |              ----+         

//    |       ------R              

//    +-------                     

//    <-------- d ------->         

//                                 

int RCol(int col, int solid);

class DetectorParameters

{

public:

  DetectorParameters(double a, double b, double d, double active, double n1, double n2, double n3, TVector3 gap, bool coupling):

    _a(a),

    _b(b),

    _d(d),

    _active(active),

    _n1(n1),

    _n2(n2),

    _n3(n3),

    _gap(gap),

    _fresnel(1),

    _guideOn(1),

    offsetY(b/2.0),

    offsetZ(b/2.0),

    _plateOn(1),

    _plateWidth(1),

    _glassOn(0),

    _glassD(0),

    _coupling(coupling)

  {};

  DetectorParameters(double a, double b, double d):

    _a(a),

    _b(b),

    _d(d),

    _active(a),

    _n1(1.0),

    _n2(1.53),

    _n3(1.46),

    _gap(TVector3(0,0,0)),

    _fresnel(1),

    _guideOn(1),

    offsetY(b/2.0),

    offsetZ(b/2.0),

    _plateOn(1),

    _plateWidth(1),

    _glassOn(0),

    _glassD(0),

    _coupling(false)

  {};

  ~DetectorParameters() {};

  void setGuide(double a, double b, double d) {

    _a = a;

    _b = b;

    //_M = b/a;

    _d = d;

  };

  void setGap(double x, double y, double z) { _gap = TVector3(x,y,z); };

  void setFresnel(int fresnel) { _fresnel = fresnel; };

  void setGlass(int glassOn, double glassD) { _glassOn = glassOn; _glassD = glassD; };

  void setGuideOn(int guideOn) { _guideOn = guideOn; };

  void setPlate(int plateOn, double plateWidth) { _plateOn = plateOn; _plateWidth = plateWidth; };

  void setIndices(double n1, double n2, double n3) { _n1 = n1; _n2 = n2; _n3 = n3; };

  void setCoupling(bool coupling) { _coupling = coupling; };

  double getLightYield() { return ( Power(_b,2)/Power(_active,2)); };

  double getM() { return (_b/_a); };

  double getA() {return _a;};

  double getB() {return _b;};

  double getD() {return _d;};

  double getN1() {return _n1; };

  double getN2() {return _n2;};

  double getN3() {return _n3;};

  double getActive() {return _active;};

  TVector3 getGap() {return _gap;};

  double getPlateWidth() {return _plateWidth; };

  int getFresnel() { return _fresnel; };

  int getGlassOn() { return _glassOn; };

  double getGlassD() { return _glassD; };

  int getGuideOn() { return _guideOn; };

  int getPlateOn() { return _plateOn; };

  double getOffsetY() { return offsetY; };

  double getOffsetZ() { return offsetZ; };

  // bad coupling in the case the small amount of grease was applied

  bool badCoupling() { return _coupling; };

private:

  double _a;

  double _b;

  double _d;

  double _active;

  double _n1;

  double _n2;

  double _n3;

  TVector3 _gap;

  int _fresnel;

  int _guideOn;

  double offsetY;

  double offsetZ;

  int _plateOn;

  double _plateWidth;

  int _glassOn;

  double _glassD;

  bool _coupling;

};

class Guide

{

public:

  Guide(TVector3 center0, DetectorParameters& parameters);

  ~Guide() {

    //for (int jk=0; jk<6; jk++) delete s_side[jk];

    //delete [] s_side;

    //delete hfate;

    //delete hnodb_all;

    //delete hnodb_exit;

    //delete hin;

    //delete hout;

  }

  Fate PropagateRay(CRay in, CRay& out, int *n_points, TVector3 *points);

  double getD() { return _d; }

  double getN1() { return _n1; };

  double getN2() { return _n2; };

  double getN3() { return _n3; };

  TH1F* GetHFate() {return &hfate;};

  TH1F* GetHNOdb_all() {return &hnodb_all;};

  TH1F* GetHNOdb_exit() {return &hnodb_exit;};

  TH2F* GetHIn() {return &hin;};

  TH2F* GetHOut() {return &hout;};

  int GetExitHits() {return (int)hfate.GetBinContent(5);};

  int GetEnteranceHits() {return (int)hfate.GetBinContent(6);};

  void GetVFate(int *out);

  int GetMAXODB() const {return MAX_REFLECTIONS;};

  void Draw(int color = 1, int width = 1);

  void DrawSkel(int color = 1, int width = 1);

private:

  Fate fate;

  double _d; // parameters

  double _n1, _n2, _n3; // refractive index: n1 above entry surface, n2 inside, n3 after exit

  double _r;

  double _absorption;

  double _A;

  bool _badCoupling;

  TRandom2 rand; // for material absorption

  CSurface s_side[6];

  CPlaneR grease;

  CSurface noCoupling;

  TVector3 center;

  TVector3 vodnik_edge[8];

  TH1F hfate;

  TH1F hnodb_all;

  TH1F hnodb_exit;

  TH2F hin;

  TH2F hout;

};

class Plate

{

public:

  Plate(DetectorParameters& parameters);

  ~Plate() { };

  void draw(int color, int width);

  void drawSkel(int color, int width);

  Fate propagateRay(CRay, CRay&, int*, TVector3*);

private:

  TVector3 plate_edge[8];

  CSurface sides[6];

};

class CDetector

{

public:

  CDetector(TVector3 center0, DetectorParameters& parameters);

  ~CDetector() {

    //delete glass;

    //delete glass_circle;

    //delete hglass;

    //delete active;

    //delete hactive;

    //delete hlaser;

    //delete detector;

    //delete hdetector;

    //delete guide;

    //delete plate;

    //delete histoPlate;

    //delete window;

    //delete window_circle;

    //delete hwindow;

  }

  //void SetLGType(int in = SURF_REFRA, int side = SURF_REFRA, int out = SURF_REFRA)

  //    {type_in = in; type_side = side; type_out = out;};

  /*    void SetLG(double SiPM0, double M0, double d0, double n10, double n20, double n30, double R0)

        { SiPM=SiPM0;

                M=M0;

                d=d0;

                n1=n10;

                n2=n20;

                n3=n30;

                R=R0;

                }; */

  //void SetR(double R0) {R = R0;};

  //void SetGlass(int glass_on0, double glass_d0)

  //{glass_on = glass_on0; glass_d = glass_d0;};

  //void SetGap(double x_gap0, double y_gap0, double z_gap0)

  //{x_gap = x_gap0; y_gap = y_gap0; z_gap = z_gap0;};

  void SetRCol(int in, int lg, int out, int gla)

  {col_in = in; col_lg = lg; col_out = out; col_rgla = gla;};

  void SetDCol(int LG0, int glass0, int active0)

  {col_LG = LG0; col_glass = glass0; col_active = active0;};

  //void SetFresnel(int b)

  //{fresnel = b;};

  //void SetAbsorption(int b, double A0)

  //{absorption = b; A = A0;};

  //void SetWindow(double wR, double d0) {window_R = wR; window_d = d0;};

  //void SetGuideOn(int b)

  //{guide_on = b;};

  //void Init();

  int Propagate(CRay, CRay&, int);

  Guide* GetLG() {return &guide;};

  //TH2F* GetHWindow() const {return hwindow;};

  TH2F* GetHGlass() {return &hglass;};

  TH2F* GetHActive() {return &hactive;};

  TH2F* GetHLaser() {return &hlaser;};

  TH2F* GetHDetector() {return &hdetector;};

  TH2F* GetHPlate() {return &histoPlate;};

  TH2F* GetGenerated() {return &generated;};

  //double GetSiPM() const {return SiPM;};

  //double GetM() const {return M;};

  //double GetD() const {return d;};

  //double GetR() const {return R;};

  //TVector3 GetVGap() const {return (TVector3(x_gap, y_gap, z_gap));};

  //int IsGlass() {return glass_on;};

  //double getGlassWidth() {return glass_d;};

  //double GetWindowR() {return window_R;};

  //double GetWindowD() {return window_d;};

  void Draw(int width = 2);

private:

  Fate fate;

  TVector3 center;

  //int type_in, type_side, type_out;

  //double SiPM, M, d, n1, n2, n3, R;

  //double detectorActive;

  int glass_on;

  double glass_d;

  CSurface glass;

  CPlaneR glass_circle;

  TH2F hglass;

  //double x_gap, y_gap, z_gap;

  CPlane4 active;

  CPlaneR grease;

  TH2F hactive;

  TH2F hlaser;

  CPlane4 detector;

  TH2F hdetector;

  int col_in, col_lg, col_out, col_rgla;

  int col_LG, col_glass, col_active;

  //int fresnel, absorption;

  //double A;

  int guide_on;

  //double window_R, window_d;

  //CSurface *window;

  //CPlaneR *window_circle;

  //TH2F *hwindow;

  Guide guide;

  Plate plate;

  double _plateWidth;

  int _plateOn;

  TH2F histoPlate;

  double offsetY;

  double offsetZ;

  TH2F generated;

};

#endif