Subversion Repositories f9daq

Rev

Rev 54 | Go to most recent revision | Details | 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
//-----------------------------------------------------------------------------
156
void PrintGuideStat(CDetector *detector, double n_in = 1)
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
 
171
        double n_active = (detector->GetHActive())->GetEntries();
172
        double r_acc = 100.0 * n_active / n_in;
173
        printf("%7.3lf\n", r_acc);
174
}
175
//-----------------------------------------------------------------------------
176
 
177
//-----------------------------------------------------------------------------
178
// en zarek
179
double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)
180
{
181
  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
182
        theta = Pi()*theta/180.0;
183
        if(theta < 1e-6) theta = 1e-6;
184
        phi = phi*Pi()/180.0;
185
 
186
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
187
 
188
        //detector->Init();
189
        if(show_3d) detector->Draw(draw_width);
190
 
191
        CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),
192
                          TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));
193
        CRay *ray1 = new CRay;
194
 
195
        detector->Propagate(*ray0, ray1, show_3d);
196
 
197
        return (detector->GetHActive())->GetEntries() / (double)(1);
198
}
199
//-----------------------------------------------------------------------------
200
// zarki, razporejeni v mrezi
201
double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)
202
{
203
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
204
        theta = Pi()*theta/180.0;
205
        if(theta < 1e-6) theta = 1e-6;
206
 
207
        //detector->Init();
208
        if(show_3d) detector->Draw(draw_width);
209
 
210
        const double b = parameters.getB();
211
 
212
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
213
 
214
        for(int i = 0; i < NN; i++) {
215
                for(int j = 0; j < NN; j++) {
216
                CRay *ray0 = new CRay(CENTER.x() - offset,
217
                              0.99*(i*b/NN + b/(2.0*NN) - b/2.0),
218
                              0.99*(j*b/NN + b/(2.0*NN) - b/2.0),
219
                                          TMath::Cos(theta),
220
                                          0.0,
221
                                          TMath::Sin(theta));
222
                CRay *ray1 = new CRay;
223
                if(i == (NN/2))
224
                        detector->Propagate(*ray0, ray1, show_3d);
225
                else
226
                        detector->Propagate(*ray0, ray1, 0);   
227
                }
228
        }
229
 
230
        return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
231
}
232
//-----------------------------------------------------------------------------
233
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
234
// vsi pod kotom (theta, phi)
235
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30)
236
{
237
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
238
        theta = theta*3.14159265358979312/180.0;
239
        if(theta < MARGIN) theta = MARGIN;
240
        phi = phi*3.14159265358979312/180.0;
241
        if(phi < MARGIN) phi = MARGIN;
242
 
243
        TDatime now;
244
        TRandom rand;
245
        rand.SetSeed(now.Get());
246
 
247
        //detector->Init();
248
        if(show_3d) detector->Draw(draw_width);
249
 
250
        double SiPM = parameters.getA();
251
        double    M = parameters.getM();
252
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
253
 
254
        for(int i = 0; i < NN; i++) {
255
 
256
                double rx = CENTER.x() - offset;
257
                double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
258
                double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
259
 
260
                double rl = TMath::Cos(theta);
261
                double rm = TMath::Sin(theta)*TMath::Sin(phi);
262
                double rn = TMath::Sin(theta)*TMath::Cos(phi); 
263
 
264
                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
265
                CRay *ray1 = new CRay;
266
 
267
                if(i < show_rays)
268
                        detector->Propagate(*ray0, ray1, show_3d);
269
                else
270
                        detector->Propagate(*ray0, ray1, 0);
271
                delete ray0;
272
                delete ray1;   
273
        }
274
 
275
        return (detector->GetHActive())->GetEntries() / (double)NN;
276
}
277
//-----------------------------------------------------------------------------
278
// zarki, izotropno porazdeljeni znotraj kota theta
279
// = nakljucni vstopni polozaj in kot
280
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)
281
{
282
//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);
283
        double pi = 3.14159265358979312;
284
        theta = theta*3.14159265358979312/180.0;
285
        if(theta < MARGIN) theta = MARGIN;
286
 
287
        TDatime now;
288
        TRandom rand;
289
        rand.SetSeed(now.Get());
290
        double rx, ry, rz, rl, rm, rn;
291
        double rphi, rtheta;
292
        //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);
293
        double theta_min_rad = TMath::Cos(theta);
294
 
295
          TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;
296
          hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;
297
          hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);    
298
          htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;
299
          htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);     
300
          hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;
301
          hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);       
302
          hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;
303
          hl = new TH1F("hl", "hl", 100, 0.0, 1.0);    
304
          hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;
305
          hm = new TH1F("hm", "hm", 100, 0.0, 1.0);    
306
          hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;
307
          hn = new TH1F("hn", "hn", 100, 0.0, 1.0);
308
 
309
 
310
        //detector->Init();
311
        if (show_3d) detector->Draw(draw_width);
312
 
313
        double SiPM = parameters.getA();
314
        double M = parameters.getM();
315
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
316
 
317
        for(int i = 0; i < NN; i++) {
318
 
319
                rx = CENTER.x() - offset;
320
                ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
321
                rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);
322
 
323
                rphi = rand.Uniform(0.0, 2.0*pi);
324
                rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
325
 
326
                rl = TMath::Cos(rtheta);
327
                rm = TMath::Sin(rtheta)*TMath::Sin(rphi);
328
                rn = TMath::Sin(rtheta)*TMath::Cos(rphi);      
329
 
330
                if(show_rand) {
331
                        htheta->Fill(rtheta);
332
                        hcostheta->Fill( TMath::Cos(rtheta) );
333
                        hphi->Fill(rphi);
334
                        hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);
335
                }
336
 
337
                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);
338
                CRay *ray1 = new CRay;
339
 
340
                if(i < show_rays) {
341
                        detector->Propagate(*ray0, ray1, show_3d);
342
                        }
343
                else {
344
                        detector->Propagate(*ray0, ray1, 0);
345
                        }
346
 
347
          delete ray0;
348
          delete ray1; 
349
        }
350
 
351
        if(show_rand) {
352
                TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");
353
                if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);
354
                else c2rand->Clear();
355
                c2rand->Divide(2,3);
356
 
357
                c2rand->cd(1); hphi->Draw();
358
                c2rand->cd(3); htheta->Draw();
359
                c2rand->cd(5); hcostheta->Draw();
360
                c2rand->cd(2); hl->Draw();
361
                c2rand->cd(4); hm->Draw();
362
                c2rand->cd(6); hn->Draw();
363
        }
364
 
365
 
366
        double acceptance = (detector->GetHActive())->GetEntries() / (double)NN;       
367
        return acceptance;
368
}
369