Subversion Repositories f9daq

Rev

Rev 73 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1. #ifndef CVODNIK_H
  2. #define CVODNIK_H
  3.  
  4. #include "TROOT.h"
  5. #include "TVector3.h"
  6. #include "TMath.h"
  7. #include "TPolyLine3D.h"
  8. #include "TRandom2.h"
  9. #include "TDatime.h"
  10. #include "TH1F.h"
  11. #include "TH2F.h"
  12. #include "TColor.h"
  13.  
  14.  
  15.  
  16. #define MARGIN 1.0e-6
  17. #define DEGREE 180.0/3.14159265358979312
  18. #define MAX_REFLECTIONS 19
  19.  
  20.  
  21. using namespace TMath;
  22.  
  23. const double c_reflectivity = 0.96;
  24. const int REFRACTION = 1;
  25. const int REFLECTION = 2;
  26. const int TOTAL = 3;
  27. const int ABSORBED = -2;
  28. const TVector3 CENTER(-2.0, 0.0, 0.0);
  29. const int nch = 50; // number of "channels" used as number of bins in histograms
  30.  
  31. enum Fate {missed=0, hitExit=1, refracted=-1, enter=2, noreflection=-2, backreflected=-3, rays=3, absorbed=4};
  32.  
  33. int dbg = 0;
  34.  
  35. //=================================================================================
  36. // zarek (premica)
  37. class CRay
  38. {
  39. public:
  40.   CRay() :
  41.     r(TVector3(0,0,0)),
  42.     k(TVector3(1,0,0)),
  43.     p(TVector3(0,1,0)),
  44.     color(kBlack)
  45.   {};
  46.   CRay(TVector3 r0, TVector3 k0) :
  47.     r(r0),
  48.     k(k0.Unit()),
  49.     p(TVector3(0,1,0)),
  50.     //p(TVector3(0,0,1)),
  51.     color(kBlack)
  52.   {};
  53.   CRay(double x0, double y0, double z0, double l0, double m0, double n0) :
  54.     r(TVector3(x0,y0,z0)),
  55.     k(TVector3(l0,m0,n0).Unit()),
  56.     p(TVector3(0,1,0)),
  57.     //p(TVector3(0,0,1)),
  58.     color(kBlack)
  59.   {};
  60.  
  61.   void Set(TVector3 r0, TVector3 k0);
  62.   //void Set(double x0, double y0, double z0, double l0, double m0, double n0);
  63.   void SetColor(int c){color = c;};
  64.   void setPolarization(TVector3 p0) {
  65.     p = p0.Unit();
  66.     if (p.Dot(k)>1e-3) printf("*** ERROR in CRay: E has component || with k\n");
  67.   };
  68.  
  69.   inline CRay & operator = (const CRay &);
  70.  
  71.   TVector3 GetR() const {return r;};
  72.   TVector3 GetK() const {return k;};
  73.   TVector3 GetP() const {return p;};
  74.  
  75.   void Print();
  76.   void Draw();
  77.   void Draw(double x_from, double x_to);
  78.   void DrawS(double x_from, double t);
  79.  
  80.   //r = position, k = unit wave vector, p = polarization
  81. private:
  82.   TVector3 r;
  83.   TVector3 k;
  84.   TVector3 p;
  85.   int color;
  86. };
  87.  
  88.  
  89. //=================================================================================
  90. // ravnina, definirana in omejena s 4 tockami na njej (stirikotni poligon)
  91. // zanasam se na taksno definicijo: (normala kaze noter)
  92. // r1----->r2
  93. // ^        |
  94. // |        |
  95. // |        |
  96. // |
  97. // r0<-----r3
  98. class CPlane4
  99. {
  100.  
  101. public:
  102.   CPlane4();
  103.   CPlane4(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);
  104.   CPlane4(TVector3 *vr);
  105.   ~CPlane4() { };
  106.   void Set(TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4);
  107.   void Set(TVector3 *vr) {Set(vr[0], vr[1], vr[2], vr[3]);};
  108.   void FlipN(){n = -n;};
  109.  
  110.   int GetIntersection(TVector3 *vec, CRay ray);
  111.   int IsInTri(TVector3 vec, TVector3 e1, TVector3 e2, double);
  112.   int IsVectorIn(TVector3 vec);
  113.   int TestIntersection(CRay in); // ray in & ??? plane
  114.   int TestIntersection(TVector3 *vec, CRay in); // plane defined with normal vec and ray in
  115.  
  116.   TVector3 GetN() {return n;};
  117.  
  118.   void Print();
  119.   void Draw(int color = 1, int width = 1);
  120.  
  121. private:
  122.   TVector3 r[4];
  123.   TVector3 n;
  124.   double A, B, C, D;
  125.   TVector3 edge[4]; // vektorji stranic
  126.   double angle_r[4]; // koti ob posameznem vogalu
  127. };
  128. //=================================================================================
  129. // ravnina - krog
  130. class CPlaneR
  131. {
  132.  
  133. public:
  134.   CPlaneR(TVector3 c, TVector3 n0, double R0)
  135.   {center = c; n = n0; _r = R0;};
  136.   CPlaneR(){ n=TVector3(0,0,0); center=TVector3(0,0,0); _r=0; };
  137.   ~CPlaneR() {};
  138.  
  139.   void Set(TVector3 c, TVector3 n0, double R0)
  140.   {center = c; n = n0; _r = R0;};
  141.  
  142.   int TestIntersection(TVector3 *vec, CRay in);
  143.  
  144.   void Draw(int color = 1, int NN = 32);
  145.  
  146. private:
  147.   TVector3 n, center;
  148.   double _r;
  149. };
  150.  
  151. //=================================================================================
  152. // ravna opticna povrsina: refractor, zrcalo ali povrsina s totalnim odbojem
  153. #define SURF_DUMMY 0
  154. #define SURF_REFRA 1
  155. #define SURF_REFLE 2
  156. #define SURF_TOTAL 3
  157. //#define SURF_TOTAL 1
  158. #define SURF_IMPER 4
  159. //#define SURF_DIFUS 5
  160.  
  161. class CSurface : public CPlane4
  162. {
  163.  
  164. public:
  165.   CSurface(int type0 = 0);
  166.   CSurface(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4,
  167.       double n10, double n20, double reflectivity);
  168.   CSurface(int type0, TVector3 *vr,  double n10, double n20, double reflectivity);
  169.   ~CSurface() {};
  170.   CSurface& operator=(const CSurface& rhs) {
  171.     type=rhs.type;
  172.     n1=rhs.n1;
  173.     n2=rhs.n2;
  174.     n1_n2=rhs.n1_n2;
  175.     reflection=rhs.reflection;
  176.     cosTtotal=rhs.cosTtotal;
  177.     fresnel=rhs.fresnel;
  178.     return *this;
  179.   }
  180.  
  181.   void SetV(TVector3 *vr){CPlane4::Set(vr);};
  182.   void SetType(int type0){type = type0;};
  183.   void SetIndex(double n10, double n20);
  184.   void SetReflection(double reflectivity){reflection = reflectivity;};
  185.   void SetFresnel(int f = 1) {fresnel = f;};
  186.  
  187.   void Set(int type0, TVector3 *vr, double n10, double n20, double reflectivity)
  188.   {type = type0; CPlane4::Set(vr); SetIndex(n10, n20); reflection = reflectivity;};
  189.   void Set(int type0, TVector3 r1, TVector3 r2, TVector3 r3, TVector3 r4, double n10, double n20, double reflectivity)
  190.   {type = type0; CPlane4::Set(r1, r2, r3, r4); SetIndex(n10, n20); reflection = reflectivity;};
  191.  
  192.   int PropagateRay(CRay in, CRay& out, TVector3 *intersection);
  193.  
  194.   double N1_N2(int sign) { return ((sign > 0) ? n1/n2 : n2/n1);};
  195.  
  196.   void Print(){printf("Type = %d\n", type); CPlane4::Print();};
  197.  
  198. private:
  199.   int type; //0 = dummy; 1 = refractor; 2 = reflector; 3 = total reflection
  200.   double n1, n2, n1_n2; //index of refraction, n1 @ +normal, n2 @ -normal
  201.   double reflection; // odbojnost stranic
  202.   TRandom2 rand; // za racunanje verjetnosti odboja od zrcala
  203.   double cosTtotal; //cosinus mejnega kota totalnega odboja za dana n1 in n2
  204.  
  205.   int fresnel; // ali naj uposteva Fresnelove enacbe
  206. };
  207.  
  208.  
  209. //=================================================================================
  210. //    +-------   n1                
  211. //    |       -------              
  212. //    |          n2  ----+        
  213. //    |                  |        
  214. //   0| M*SiPM           | SiPM    
  215. //    |                  |        
  216. //    |              ----+        
  217. //    |       ------R              
  218. //    +-------                    
  219. //    <-------- d ------->        
  220. //                                
  221.  
  222.  
  223.  
  224. int RCol(int col, int solid);
  225.  
  226.  
  227. class DetectorParameters
  228. {
  229. public:
  230.   DetectorParameters(double a, double b, double d, double active, double n1, double n2, double n3, TVector3 gap, bool coupling):
  231.     _a(a),
  232.     _b(b),
  233.     _d(d),
  234.     _active(active),
  235.     _n1(n1),
  236.     _n2(n2),
  237.     _n3(n3),
  238.     _gap(gap),
  239.     _fresnel(1),
  240.     _guideOn(1),
  241.     offsetY(b/2.0),
  242.     offsetZ(b/2.0),
  243.     _plateOn(1),
  244.     _plateWidth(1),
  245.     _glassOn(0),
  246.     _glassD(0),
  247.     _coupling(coupling)
  248.   {};
  249.   DetectorParameters(double a, double b, double d):
  250.     _a(a),
  251.     _b(b),
  252.     _d(d),
  253.     _active(a),
  254.     _n1(1.0),
  255.     _n2(1.53),
  256.     _n3(1.46),
  257.     _gap(TVector3(0,0,0)),
  258.     _fresnel(1),
  259.     _guideOn(1),
  260.     offsetY(b/2.0),
  261.     offsetZ(b/2.0),
  262.     _plateOn(1),
  263.     _plateWidth(1),
  264.     _glassOn(0),
  265.     _glassD(0),
  266.     _coupling(false)
  267.   {};
  268.   ~DetectorParameters() {};
  269.  
  270.   void setGuide(double a, double b, double d) {
  271.     _a = a;
  272.     _b = b;
  273.     //_M = b/a;
  274.     _d = d;
  275.   };
  276.   void setGap(double x, double y, double z) { _gap = TVector3(x,y,z); };
  277.   void setFresnel(int fresnel) { _fresnel = fresnel; };
  278.   void setGlass(int glassOn, double glassD) { _glassOn = glassOn; _glassD = glassD; };
  279.   void setGuideOn(int guideOn) { _guideOn = guideOn; };
  280.   void setPlate(int plateOn, double plateWidth) { _plateOn = plateOn; _plateWidth = plateWidth; };
  281.   void setIndices(double n1, double n2, double n3) { _n1 = n1; _n2 = n2; _n3 = n3; };
  282.   void setCoupling(bool coupling) { _coupling = coupling; };
  283.  
  284.   double getLightYield() { return ( Power(_b,2)/Power(_active,2)); };
  285.   double getM() { return (_b/_a); };
  286.   double getA() {return _a;};
  287.   double getB() {return _b;};
  288.   double getD() {return _d;};
  289.   double getN1() {return _n1; };
  290.   double getN2() {return _n2;};
  291.   double getN3() {return _n3;};
  292.   double getActive() {return _active;};
  293.   TVector3 getGap() {return _gap;};
  294.   double getPlateWidth() {return _plateWidth; };
  295.   int getFresnel() { return _fresnel; };
  296.   int getGlassOn() { return _glassOn; };
  297.   double getGlassD() { return _glassD; };
  298.   int getGuideOn() { return _guideOn; };
  299.   int getPlateOn() { return _plateOn; };
  300.   double getOffsetY() { return offsetY; };
  301.   double getOffsetZ() { return offsetZ; };
  302.   // bad coupling in the case the small amount of grease was applied
  303.   bool badCoupling() { return _coupling; };
  304.  
  305. private:
  306.   double _a;
  307.   double _b;
  308.   double _d;
  309.   double _active;
  310.   double _n1;
  311.   double _n2;
  312.   double _n3;
  313.   TVector3 _gap;
  314.  
  315.   int _fresnel;
  316.  
  317.   int _guideOn;
  318.   double offsetY;
  319.   double offsetZ;
  320.  
  321.   int _plateOn;
  322.   double _plateWidth;
  323.  
  324.   int _glassOn;
  325.   double _glassD;
  326.   bool _coupling;
  327.  
  328. };
  329.  
  330. class Guide
  331. {
  332. public:
  333.   Guide(TVector3 center0, DetectorParameters& parameters);
  334.   ~Guide() {
  335.     //for (int jk=0; jk<6; jk++) delete s_side[jk];
  336.     //delete [] s_side;
  337.     //delete hfate;
  338.     //delete hnodb_all;
  339.     //delete hnodb_exit;
  340.     //delete hin;
  341.     //delete hout;
  342.   }
  343.  
  344.   Fate PropagateRay(CRay in, CRay& out, int *n_points, TVector3 *points);
  345.  
  346.   double getD() { return _d; }
  347.   double getN1() { return _n1; };
  348.   double getN2() { return _n2; };
  349.   double getN3() { return _n3; };
  350.   TH1F* GetHFate() {return &hfate;};
  351.   TH1F* GetHNOdb_all() {return &hnodb_all;};
  352.   TH1F* GetHNOdb_exit() {return &hnodb_exit;};
  353.   TH2F* GetHIn() {return &hin;};
  354.   TH2F* GetHOut() {return &hout;};
  355.   int GetExitHits() {return (int)hfate.GetBinContent(5);};
  356.   int GetEnteranceHits() {return (int)hfate.GetBinContent(6);};
  357.   void GetVFate(int *out);
  358.   int GetMAXODB() const {return MAX_REFLECTIONS;};
  359.  
  360.   void Draw(int color = 1, int width = 1);
  361.   void DrawSkel(int color = 1, int width = 1);
  362.  
  363. private:
  364.   Fate fate;
  365.   double _d; // parameters
  366.   double _n1, _n2, _n3; // refractive index: n1 above entry surface, n2 inside, n3 after exit
  367.   double _r;
  368.   double _absorption;
  369.   double _A;
  370.   bool _badCoupling;
  371.   TRandom2 rand; // for material absorption
  372.   CSurface s_side[6];
  373.   CPlaneR grease;
  374.   CSurface noCoupling;
  375.   TVector3 center;
  376.   TVector3 vodnik_edge[8];
  377.  
  378.   TH1F hfate;
  379.   TH1F hnodb_all;
  380.   TH1F hnodb_exit;
  381.   TH2F hin;
  382.   TH2F hout;
  383. };
  384.  
  385.  
  386.  
  387. class Plate
  388. {
  389. public:
  390.   Plate(DetectorParameters& parameters);
  391.   ~Plate() { };
  392.  
  393.   void draw(int color, int width);
  394.   void drawSkel(int color, int width);
  395.   Fate propagateRay(CRay, CRay&, int*, TVector3*);
  396.  
  397. private:
  398.   TVector3 plate_edge[8];
  399.   CSurface sides[6];
  400. };
  401.  
  402.  
  403. class CDetector
  404. {
  405. public:
  406.   CDetector(TVector3 center0, DetectorParameters& parameters);
  407.   ~CDetector() {
  408.     //delete glass;
  409.     //delete glass_circle;
  410.     //delete hglass;
  411.     //delete active;
  412.     //delete hactive;
  413.     //delete hlaser;
  414.     //delete detector;
  415.     //delete hdetector;
  416.     //delete guide;
  417.     //delete plate;
  418.     //delete histoPlate;
  419.     //delete window;
  420.     //delete window_circle;
  421.     //delete hwindow;
  422.   }
  423.  
  424.   //void SetLGType(int in = SURF_REFRA, int side = SURF_REFRA, int out = SURF_REFRA)
  425.   //    {type_in = in; type_side = side; type_out = out;};
  426.  
  427.   /*    void SetLG(double SiPM0, double M0, double d0, double n10, double n20, double n30, double R0)
  428.         { SiPM=SiPM0;
  429.                 M=M0;
  430.                 d=d0;
  431.                 n1=n10;
  432.                 n2=n20;
  433.                 n3=n30;
  434.                 R=R0;
  435.                 }; */
  436.   //void SetR(double R0) {R = R0;};
  437.   //void SetGlass(int glass_on0, double glass_d0)
  438.   //{glass_on = glass_on0; glass_d = glass_d0;};
  439.   //void SetGap(double x_gap0, double y_gap0, double z_gap0)
  440.   //{x_gap = x_gap0; y_gap = y_gap0; z_gap = z_gap0;};
  441.   void SetRCol(int in, int lg, int out, int gla)
  442.   {col_in = in; col_lg = lg; col_out = out; col_rgla = gla;};
  443.   void SetDCol(int LG0, int glass0, int active0)
  444.   {col_LG = LG0; col_glass = glass0; col_active = active0;};
  445.   //void SetFresnel(int b)
  446.   //{fresnel = b;};
  447.   //void SetAbsorption(int b, double A0)
  448.   //{absorption = b; A = A0;};
  449.   //void SetWindow(double wR, double d0) {window_R = wR; window_d = d0;};
  450.   //void SetGuideOn(int b)
  451.   //{guide_on = b;};
  452.  
  453.   //void Init();
  454.  
  455.   int Propagate(CRay, CRay&, int);
  456.  
  457.   Guide* GetLG() {return &guide;};
  458.   //TH2F* GetHWindow() const {return hwindow;};
  459.   TH2F* GetHGlass() {return &hglass;};
  460.   TH2F* GetHActive() {return &hactive;};
  461.   TH2F* GetHLaser() {return &hlaser;};
  462.   TH2F* GetHDetector() {return &hdetector;};
  463.   TH2F* GetHPlate() {return &histoPlate;};
  464.   TH2F* GetGenerated() {return &generated;};
  465.  
  466.   //double GetSiPM() const {return SiPM;};
  467.   //double GetM() const {return M;};
  468.   //double GetD() const {return d;};
  469.   //double GetR() const {return R;};
  470.   //TVector3 GetVGap() const {return (TVector3(x_gap, y_gap, z_gap));};
  471.   //int IsGlass() {return glass_on;};
  472.   //double getGlassWidth() {return glass_d;};
  473.   //double GetWindowR() {return window_R;};
  474.   //double GetWindowD() {return window_d;};
  475.  
  476.   void Draw(int width = 2);
  477.  
  478. private:
  479.   Fate fate;
  480.   TVector3 center;
  481.  
  482.   //int type_in, type_side, type_out;
  483.   //double SiPM, M, d, n1, n2, n3, R;
  484.   //double detectorActive;
  485.  
  486.  
  487.   int glass_on;
  488.   double glass_d;
  489.   CSurface glass;
  490.   CPlaneR glass_circle;
  491.   TH2F hglass;
  492.  
  493.   //double x_gap, y_gap, z_gap;
  494.   CPlane4 active;
  495.   CPlaneR grease;
  496.   TH2F hactive;
  497.   TH2F hlaser;
  498.  
  499.   CPlane4 detector;
  500.   TH2F hdetector;
  501.  
  502.   int col_in, col_lg, col_out, col_rgla;
  503.   int col_LG, col_glass, col_active;
  504.  
  505.   //int fresnel, absorption;
  506.   //double A;
  507.  
  508.   int guide_on;
  509.  
  510.   //double window_R, window_d;
  511.   //CSurface *window;
  512.   //CPlaneR *window_circle;
  513.   //TH2F *hwindow;
  514.  
  515.   Guide guide;
  516.   Plate plate;
  517.  
  518.   double _plateWidth;
  519.   int _plateOn;
  520.  
  521.   TH2F histoPlate;
  522.  
  523.   double offsetY;
  524.   double offsetZ;
  525.  
  526.   TH2F generated;
  527. };
  528.  
  529.  
  530.  
  531. #endif
  532.