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