Subversion Repositories f9daq

Rev

Rev 25 | Rev 70 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
25 f9daq 1
// own include
2
//#include "include/raySimulator.h"
3
#include "include/RTUtil.h"
4
#include "include/guide.h"
5
 
6
// general include
7
#include "TROOT.h"
8
#include "TSystem.h"
9
#include "TStyle.h"
10
#include "TCanvas.h"
11
#include "TMath.h"
12
#include "TView3D.h"
13
#include "TH1F.h"
14
#include "TH2F.h"
15
#include "TBenchmark.h"
16
#include "TPolyMarker.h"
17
#include "TGraph.h"
18
#include "TF1.h"
19
 
20
//=================================================================================
21
TCanvas *c3dview;
22
TCanvas *cacc;
23
int show_3d = 0, show_data = 0;
24
int draw_width = 2;
25
 
26
 // calls default constructor CVodnik.cpp
27
 
28
//void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0)
29
        //{center.SetXYZ(x,y,z); detector = new CDetector(center);}
30
 
31
 
32
//void SetLGType(int in = 1, int side = 1, int out = 0)
33
        //{detector->SetLGType(in, side, out);}
34
 
35
//void SetLG(double SiPM0 = 3.0, double M0 = 1.666, double d0 = 3.0, double n10 = 1.0, double n20 = 1.48, double n30 = 1.48, double R0 = 0.96)
36
        //{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);}
37
 
38
//void SetR(double R0) {detector->SetR(R0);}
39
        /*
40
void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3)
41
        {detector->SetGlass(glass_on0, glass_d0);}
42
 
43
void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)
44
        {detector->SetGap(x_gap0 , y_gap0, z_gap0);}
45
 
46
void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6)
47
        {detector->SetRCol(in, lg, out, gla);};
48
void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3)
49
        {detector->SetDCol(LG0, glass0, active0);};
50
 
51
 
52
void SetDWidth(int w = 1) {draw_width = w;}
53
 
54
void SetFresnel(int b = 1) {detector->SetFresnel(b);}
55
void SetAbsorb(int b = 0, double A0 = 1e4) {detector->SetAbsorption(b, A0);}
56
 
57
void SetGuideOn(int b = 0) {detector->SetGuideOn(b);}
58
 
59
void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);}
60
void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}
61
*/
62
//-----------------------------------------------------------------------------
63
void Init(double rsys_scale = 7.0)
64
{      
65
        RTSetStyle(gStyle);
66
        gStyle->SetOptStat(10);
67
 
68
        if(show_3d) {
69
                double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};
70
                double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};
71
 
72
                c3dview = (TCanvas*)gROOT->FindObject("c3dview");
73
                if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);      c3dview->SetFillColor(0);}
74
                else c3dview->Clear();
75
 
76
                TView3D *view = new TView3D(1, rsys0, rsys1);
77
                view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);
78
                //view->SetPerspective();
79
                view->SetParallel();
80
                view->FrontView();
81
                view->Zoom();view->Zoom();view->Zoom();//view->Zoom();
82
                //view->ToggleRulers();
83
        }
84
}
85
//-----------------------------------------------------------------------------
86
void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)
87
{
88
        if(show_data) {
89
                char sbuff[256];
90
                sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
91
                                parameters.getA(), parameters.getM(), parameters.getD(),
92
                                parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);
93
 
94
                if(!only2d) {
95
                        RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);
96
                        cdata->Divide(3,3);
97
                        int cpc = 1;
98
                        cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");
99
                        cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");
100
                        cdata->cd(cpc++); ((detector->GetLG())->GetHOut())->Draw("COLZ");
101
 
102
 
103
                        cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");
104
                        cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");
105
                        cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");
106
 
107
                        cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();
108
                        cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();
109
                        cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();
110
                }else{         
111
                        RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);
112
                        cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");
113
                        cdata->SaveAs("cdata.gif");
114
                }
115
        }
116
}
117
//-----------------------------------------------------------------------------
118
void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",
119
                         double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)
120
{
121
 
122
        cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);
123
 
124
        (cacc->cd(0))->SetRightMargin(0.05);
125
        (cacc->cd(0))->SetLeftMargin(0.13);
126
        (cacc->cd(0))->SetTopMargin(0.08);
127
 
128
        TGraph *gacceptance = new TGraph(n, x, y);
129
        gacceptance->SetTitle(htitle);
130
        gacceptance->SetMarkerStyle(8);
131
        gacceptance->SetLineColor(12);
132
        gacceptance->SetLineWidth(1);
133
        if(xmin < xmax)
134
                (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);
135
        (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);
136
 
137
        gacceptance->Draw("ACP");
138
        //cacc->Print("acceptance.pdf");
139
        //delete gacceptance;
140
        //delete cacc;
141
}
142
//-----------------------------------------------------------------------------
143
void PrintGlassStat(CDetector *detector)
144
{
145
        int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();
146
        int glass_out = (detector->GetHGlass())->GetEntries();
147
        printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);
148
}
149
//-----------------------------------------------------------------------------
150
void PrintGuideHead()
151
{
152
        //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
153
        printf("Acceptance\n");
154
}
155
//-----------------------------------------------------------------------------
54 f9daq 156
void PrintGuideStat(double acceptance)
25 f9daq 157
{
158
        /*int n_active = (detector->GetHActive())->GetEntries();
159
        int n_enter = (detector->GetLG())->GetEnteranceHits();
160
        double izkoristek = n_active/(double)n_enter;
161
        int fates[7]; (detector->GetLG())->GetVFate(fates);
162
 
163
//      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");
164
 
165
        printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |",
166
                   (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi);
167
        for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter);       
168
 
169
        printf("%5.1lf\n", 100.0*izkoristek);*/
170
 
54 f9daq 171
        //double n_active = (detector->GetHActive())->GetEntries();
172
        //double r_acc = 100.0 * n_active / n_in;
173
        double r_acc = 100.0 * acceptance;
25 f9daq 174
        printf("%7.3lf\n", r_acc);
175
}
176
//-----------------------------------------------------------------------------
177
 
178
//-----------------------------------------------------------------------------
179
// en zarek
180
double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
181
{
182
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
183
        theta = Pi()*theta/180.0;
184
        if(theta < 1e-6) theta = 1e-6;
185
        phi = phi*Pi()/180.0;
186
 
187
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
188
 
189
        //detector->Init();
190
        if(show_3d) detector->Draw(draw_width);
191
 
192
        CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
193
                          TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));
194
        CRay *ray1 = new CRay;
195
 
196
        detector->Propagate(*ray0, ray1, show_3d);
197
 
198
        return (detector->GetHActive())->GetEntries() / (double)(1);
199
}
200
//-----------------------------------------------------------------------------
201
// zarki, razporejeni v mrezi
202
double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)
203
{
204
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
205
        theta = Pi()*theta/180.0;
206
        if(theta < 1e-6) theta = 1e-6;
207
 
208
        //detector->Init();
209
        if(show_3d) detector->Draw(draw_width);
210
 
211
        const double b = parameters.getB();
212
 
213
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
214
 
215
        for(int i = 0; i < NN; i++) {
216
                for(int j = 0; j < NN; j++) {
217
                CRay *ray0 = new CRay(CENTER.x() - offset,
218
                              0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
219
                              0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
220
                                          TMath::Cos(theta),
221
                                          0.0,
222
                                          TMath::Sin(theta));
223
                CRay *ray1 = new CRay;
224
                if(i == (NN/2))
225
                        detector->Propagate(*ray0, ray1, show_3d);
226
                else
227
                        detector->Propagate(*ray0, ray1, 0);   
228
                }
229
        }
230
 
231
        return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
232
}
233
//-----------------------------------------------------------------------------
234
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
235
// vsi pod kotom (theta, phi)
236
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30)
237
{
238
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
239
        theta = theta*3.14159265358979312/180.0;
240
        if(theta < MARGIN) theta = MARGIN;
241
        phi = phi*3.14159265358979312/180.0;
242
        if(phi < MARGIN) phi = MARGIN;
243
 
244
        TDatime now;
245
        TRandom rand;
246
        rand.SetSeed(now.Get());
247
 
248
        //detector->Init();
249
        if(show_3d) detector->Draw(draw_width);
250
 
251
        double SiPM = parameters.getA();
252
        double    M = parameters.getM();
253
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
254
 
255
        for(int i = 0; i < NN; i++) {
256
 
257
                double rx = CENTER.x() - offset;
258
                double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
259
                double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
260
 
261
                double rl = TMath::Cos(theta);
262
                double rm = TMath::Sin(theta)*TMath::Sin(phi);
263
                double rn = TMath::Sin(theta)*TMath::Cos(phi); 
264
 
265
                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
266
                CRay *ray1 = new CRay;
267
 
268
                if(i < show_rays)
269
                        detector->Propagate(*ray0, ray1, show_3d);
270
                else
271
                        detector->Propagate(*ray0, ray1, 0);
54 f9daq 272
                //delete ray0;
273
                //delete ray1;  
25 f9daq 274
        }
275
 
276
        return (detector->GetHActive())->GetEntries() / (double)NN;
277
}
278
//-----------------------------------------------------------------------------
279
// zarki, izotropno porazdeljeni znotraj kota theta
280
// = nakljucni vstopni polozaj in kot
281
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)
282
{
283
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
284
        double pi = 3.14159265358979312;
285
        theta = theta*3.14159265358979312/180.0;
286
        if(theta < MARGIN) theta = MARGIN;
287
 
288
        TDatime now;
289
        TRandom rand;
290
        rand.SetSeed(now.Get());
291
        double rx, ry, rz, rl, rm, rn;
292
        double rphi, rtheta;
293
        //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
294
        double theta_min_rad = TMath::Cos(theta);
295
 
296
          TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
297
          hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
298
          hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);    
299
          htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
300
          htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);     
301
          hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
302
          hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);       
303
          hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
304
          hl = new TH1F("hl", "hl", 100, 0.0, 1.0);    
305
          hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
306
          hm = new TH1F("hm", "hm", 100, 0.0, 1.0);    
307
          hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
308
          hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
309
 
310
 
311
        //detector->Init();
312
        if (show_3d) detector->Draw(draw_width);
313
 
314
        double SiPM = parameters.getA();
315
        double M = parameters.getM();
316
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
317
 
318
        for(int i = 0; i < NN; i++) {
319
 
320
                rx = CENTER.x() - offset;
321
                ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
322
                rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
323
 
324
                rphi = rand.Uniform(0.0, 2.0*pi);
325
                rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
326
 
327
                rl = TMath::Cos(rtheta);
328
                rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
329
                rn = TMath::Sin(rtheta)*TMath::Cos(rphi);      
330
 
331
                if(show_rand) {
332
                        htheta->Fill(rtheta);
333
                        hcostheta->Fill( TMath::Cos(rtheta) );
334
                        hphi->Fill(rphi);
335
                        hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
336
                }
337
 
338
                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
339
                CRay *ray1 = new CRay;
340
 
341
                if(i < show_rays) {
342
                        detector->Propagate(*ray0, ray1, show_3d);
343
                        }
344
                else {
345
                        detector->Propagate(*ray0, ray1, 0);
346
                        }
347
 
54 f9daq 348
          //delete ray0;
349
          //delete ray1;        
25 f9daq 350
        }
351
 
352
        if(show_rand) {
353
                TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
354
                if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
355
                else c2rand->Clear();
356
                c2rand->Divide(2,3);
357
 
358
                c2rand->cd(1); hphi->Draw();
359
                c2rand->cd(3); htheta->Draw();
360
                c2rand->cd(5); hcostheta->Draw();
361
                c2rand->cd(2); hl->Draw();
362
                c2rand->cd(4); hm->Draw();
363
                c2rand->cd(6); hn->Draw();
364
        }
365
 
366
 
367
        double acceptance = (detector->GetHActive())->GetEntries() / (double)NN;       
368
        return acceptance;
369
}
370