Subversion Repositories f9daq

// own include

//#include "include/raySimulator.h"

#include "include/RTUtil.h"

#include "include/guide.h"

// general include

#include "TROOT.h"

#include "TSystem.h"

#include "TStyle.h"

#include "TCanvas.h"

#include "TMath.h"

#include "TView3D.h"

#include "TH1F.h"

#include "TH2F.h"

#include "TBenchmark.h"

#include "TPolyMarker.h"

#include "TGraph.h"

#include "TF1.h"

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

TCanvas *c3dview;

TCanvas *cacc;

int show_3d = 0, show_data = 0;

int draw_width = 2;

 // calls default constructor CVodnik.cpp

//void SetCenter(double x = -2.0, double y = 0.0, double z = 0.0)

        //{center.SetXYZ(x,y,z); detector = new CDetector(center);}

//void SetLGType(int in = 1, int side = 1, int out = 0)

        //{detector->SetLGType(in, side, out);}

//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)

        //{detector->SetLG(SiPM0, M0, d0, n10, n20, n30, R0);}

//void SetR(double R0) {detector->SetR(R0);}

        /*

void SetGlass(int glass_on0 = 0, double glass_d0 = 0.3)

        {detector->SetGlass(glass_on0, glass_d0);}

void SetGap(double x_gap0 = 0.3, double y_gap0 = 0.0, double z_gap0 = 0.0)

        {detector->SetGap(x_gap0 , y_gap0, z_gap0);}

void SetRCol(int in = 2, int lg = 8, int out = 4, int gla = 6)

        {detector->SetRCol(in, lg, out, gla);};

void SetDCol(int LG0 = 1, int glass0 = 2, int active0 = 3)

        {detector->SetDCol(LG0, glass0, active0);};

void SetDWidth(int w = 1) {draw_width = w;}

void SetFresnel(int b = 1) {detector->SetFresnel(b);}

void SetAbsorb(int b = 0, double A0 = 1e4) {detector->SetAbsorption(b, A0);}

void SetGuideOn(int b = 0) {detector->SetGuideOn(b);}

void SetTypeU() {SetGlass(1, 0.3);SetGap(0.7);}

void SetTypeC() {SetGlass(0, 0.3);SetGap(1.03);}

*/

//-----------------------------------------------------------------------------

void Init(double rsys_scale = 7.0)

{      

        RTSetStyle(gStyle);

        gStyle->SetOptStat(10);

        if(show_3d) {

                double rsys0[]={-rsys_scale,-rsys_scale,-rsys_scale};

                double rsys1[]={ rsys_scale, rsys_scale, rsys_scale};

                c3dview = (TCanvas*)gROOT->FindObject("c3dview");

                if(!c3dview) {c3dview = new TCanvas("c3dview", "3D View", 0, 0, 700, 723);      c3dview->SetFillColor(0);}

                else c3dview->Clear();

                TView3D *view = new TView3D(1, rsys0, rsys1);

                view->SetRange(rsys0[0], rsys0[1], rsys0[2], rsys1[0], rsys1[1], rsys1[2]);

                //view->SetPerspective();

                view->SetParallel();

                view->FrontView();

                view->Zoom();view->Zoom();view->Zoom();//view->Zoom();

                //view->ToggleRulers();

        }

}

//-----------------------------------------------------------------------------

void DrawData(CDetector *detector, DetectorParameters& parameters, double theta = 0.0, double acc = 0.0, int only2d = 0)

{

        if(show_data) {

                char sbuff[256];

                sprintf(sbuff, "SiPM = %.1lf, M = %.1lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",

                                parameters.getA(), parameters.getM(), parameters.getD(),

                                parameters.getGap().x(), parameters.getGap().y(), parameters.getGap().z(), theta, acc);

                if(!only2d) {

                        RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 950, 950);

                        cdata->Divide(3,3);

                        int cpc = 1;

                        cdata->cd(cpc++); (detector->histoPlate)->Draw("COLZ");

                        cdata->cd(cpc++); ((detector->GetLG())->GetHIn())->Draw("COLZ");

                        cdata->cd(cpc++); ((detector->GetLG())->GetHOut())->Draw("COLZ");

                        cdata->cd(cpc++); (detector->GetHActive())->Draw("COLZ");

                        cdata->cd(cpc++); (detector->GetHLaser())->Draw("COLZ");

                        cdata->cd(cpc++); (detector->GetHDetector())->Draw("COLZ");

                        cdata->cd(cpc++); ((detector->GetLG())->GetHFate())->Draw();

                        cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_all())->Draw();

                        cdata->cd(cpc++); ((detector->GetLG())->GetHNOdb_exit())->Draw();

                }else{         

                        RTCanvas *cdata = new RTCanvas((char*)"Data", sbuff, 350, 0, 500, 500);

                        cdata->cd(0); (detector->GetHLaser())->Draw("COLZ");

                        cdata->SaveAs("cdata.gif");

                }

        }

}

//-----------------------------------------------------------------------------

void DrawAcc(int n, double *x, double *y, char *htitle = (char*)";;Acceptance",

                         double xmin = 1.0, double xmax = 0.0, double ymin = 0.0, double ymax = 1.0)

{

        cacc = new TCanvas((char*)"Acceptance", (char*)"Acceptance", 0, 0, 700, 480);

        (cacc->cd(0))->SetRightMargin(0.05);

        (cacc->cd(0))->SetLeftMargin(0.13);

        (cacc->cd(0))->SetTopMargin(0.08);

        TGraph *gacceptance = new TGraph(n, x, y);

        gacceptance->SetTitle(htitle);

        gacceptance->SetMarkerStyle(8);

        gacceptance->SetLineColor(12);

        gacceptance->SetLineWidth(1);

        if(xmin < xmax)

                (gacceptance->GetXaxis())->SetRangeUser(xmin, xmax);

        (gacceptance->GetYaxis())->SetRangeUser(ymin, ymax);

        gacceptance->Draw("ACP");

        //cacc->Print("acceptance.pdf");

        //delete gacceptance;

        //delete cacc;

}

//-----------------------------------------------------------------------------

void PrintGlassStat(CDetector *detector)

{

        int glass_in = ((detector->GetLG())->GetHOut())->GetEntries();

        int glass_out = (detector->GetHGlass())->GetEntries();

        printf("glass_in = %d, glass_out = %d, frac = %.3lf\n", glass_in, glass_out, glass_out/(double)glass_in);

}

//-----------------------------------------------------------------------------

void PrintGuideHead()

{

        //printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");

        printf("Acceptance\n");

}

//-----------------------------------------------------------------------------

void PrintGuideStat(double acceptance)

{

        /*int n_active = (detector->GetHActive())->GetEntries();

        int n_enter = (detector->GetLG())->GetEnteranceHits();

        double izkoristek = n_active/(double)n_enter;

        int fates[7]; (detector->GetLG())->GetVFate(fates);

//      printf("   gap (x,y,z) | theta   phi | Back | Loss | Total| Miss | Exit | Acceptance\n");

        printf("%4.1lf,%4.1lf,%4.1lf | %4.1lf,%6.1lf |",

                   (detector->GetVGap()).x(), (detector->GetVGap()).y(), (detector->GetVGap()).z(), theta, phi);

        for(int i=0;i<5;i++) printf("%5.1lf |", 100.0*fates[i]/(double)n_enter);       

        printf("%5.1lf\n", 100.0*izkoristek);*/

        //double n_active = (detector->GetHActive())->GetEntries();

        //double r_acc = 100.0 * n_active / n_in;

        double r_acc = 100.0 * acceptance;

        printf("%7.3lf\n", r_acc);

}

//-----------------------------------------------------------------------------

//-----------------------------------------------------------------------------

// en zarek

double Single(CDetector *detector, DetectorParameters& parameters, TVector3 hit = TVector3(0.0, 0.0, 0.0), double theta = 0.0, double phi = 0.0)

{

  //CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);

        theta = Pi()*theta/180.0;

        if(theta < 1e-6) theta = 1e-6;

        phi = phi*Pi()/180.0;

        double offset = (parameters._plateOn ? parameters._plateWidth : 0);

        //detector->Init();

        if(show_3d) detector->Draw(draw_width);

        CRay *ray0 = new CRay(hit.x() - offset, hit.y(), hit.z(),

                          TMath::Cos(theta), TMath::Sin(theta)*TMath::Sin(phi), TMath::Sin(theta)*TMath::Cos(phi));

        CRay *ray1 = new CRay;

        detector->Propagate(*ray0, ray1, show_3d);

        return (detector->GetHActive())->GetEntries() / (double)(1);

}

//-----------------------------------------------------------------------------

// zarki, razporejeni v mrezi

double Grid(CDetector *detector, DetectorParameters& parameters, int NN = 10, double theta = 0.0)

{

//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);

        theta = Pi()*theta/180.0;

        if(theta < 1e-6) theta = 1e-6;

        //detector->Init();

        if(show_3d) detector->Draw(draw_width);

        const double b = parameters.getB();

        double offset = (parameters._plateOn ? parameters._plateWidth : 0);

        for(int i = 0; i < NN; i++) {

                for(int j = 0; j < NN; j++) {

                CRay *ray0 = new CRay(CENTER.x() - offset,

                              0.99*(i*b/NN + b/(2.0*NN) - b/2.0),

                              0.99*(j*b/NN + b/(2.0*NN) - b/2.0),

                                          TMath::Cos(theta),

                                          0.0,

                                          TMath::Sin(theta));

                CRay *ray1 = new CRay;

                if(i == (NN/2))

                        detector->Propagate(*ray0, ray1, show_3d);

                else

                        detector->Propagate(*ray0, ray1, 0);   

                }

        }

        return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));

}

//-----------------------------------------------------------------------------

// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)

// vsi pod kotom (theta, phi)

double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30)

{

//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);

        theta = theta*3.14159265358979312/180.0;

        if(theta < MARGIN) theta = MARGIN;

        phi = phi*3.14159265358979312/180.0;

        if(phi < MARGIN) phi = MARGIN;

        TDatime now;

        TRandom rand;

        rand.SetSeed(now.Get());

        //detector->Init();

        if(show_3d) detector->Draw(draw_width);

        double SiPM = parameters.getA();

        double    M = parameters.getM();

        double offset = (parameters._plateOn ? parameters._plateWidth : 0);

        for(int i = 0; i < NN; i++) {

                double rx = CENTER.x() - offset;

                double ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);

                double rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);

                double rl = TMath::Cos(theta);

                double rm = TMath::Sin(theta)*TMath::Sin(phi);

                double rn = TMath::Sin(theta)*TMath::Cos(phi); 

                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);

                CRay *ray1 = new CRay;

                if(i < show_rays)

                        detector->Propagate(*ray0, ray1, show_3d);

                else

                        detector->Propagate(*ray0, ray1, 0);

                //delete ray0;

                //delete ray1;  

        }

        return (detector->GetHActive())->GetEntries() / (double)NN;

}

//-----------------------------------------------------------------------------

// zarki, izotropno porazdeljeni znotraj kota theta

// = nakljucni vstopni polozaj in kot

double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)

{

//CDetector *detector = new CDetector(center, 3, 1.666, 3, 3, 1, 1.48, 1.48);

        double pi = 3.14159265358979312;

        theta = theta*3.14159265358979312/180.0;

        if(theta < MARGIN) theta = MARGIN;

        TDatime now;

        TRandom rand;

        rand.SetSeed(now.Get());

        double rx, ry, rz, rl, rm, rn;

        double rphi, rtheta;

        //double theta_min_rad = TMath::Power(TMath::Cos(theta), 2);

        double theta_min_rad = TMath::Cos(theta);

          TH1F *hphi, *htheta, *hcostheta, *hl, *hm, *hn;

          hphi = (TH1F*)gROOT->FindObject("hphi"); if(hphi) delete hphi;

          hphi = new TH1F("hphi", "hphi", 100, 0.0, 2.0*pi);    

          htheta = (TH1F*)gROOT->FindObject("htheta"); if(htheta) delete htheta;

          htheta = new TH1F("htheta", "htheta", 100, 0.0, pi/2.0);     

          hcostheta = (TH1F*)gROOT->FindObject("hcostheta"); if(hcostheta) delete hcostheta;

          hcostheta = new TH1F("hcostheta", "hcostheta", 100, 0.0, 1.0);       

          hl = (TH1F*)gROOT->FindObject("hl"); if(hl) delete hl;

          hl = new TH1F("hl", "hl", 100, 0.0, 1.0);    

          hm = (TH1F*)gROOT->FindObject("hm"); if(hm) delete hm;

          hm = new TH1F("hm", "hm", 100, 0.0, 1.0);    

          hn = (TH1F*)gROOT->FindObject("hn"); if(hn) delete hn;

          hn = new TH1F("hn", "hn", 100, 0.0, 1.0);

        //detector->Init();

        if (show_3d) detector->Draw(draw_width);

        double SiPM = parameters.getA();

        double M = parameters.getM();

        double offset = (parameters._plateOn ? parameters._plateWidth : 0);

        for(int i = 0; i < NN; i++) {

                rx = CENTER.x() - offset;

                ry = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);

                rz = rand.Uniform(-M*SiPM/2.0, M*SiPM/2.0);

                rphi = rand.Uniform(0.0, 2.0*pi);

                rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));

                rl = TMath::Cos(rtheta);

                rm = TMath::Sin(rtheta)*TMath::Sin(rphi);

                rn = TMath::Sin(rtheta)*TMath::Cos(rphi);      

                if(show_rand) {

                        htheta->Fill(rtheta);

                        hcostheta->Fill( TMath::Cos(rtheta) );

                        hphi->Fill(rphi);

                        hl->Fill(rl); hm->Fill(rm); hn->Fill(rn);

                }

                CRay *ray0 = new CRay(rx, ry, rz, rl, rm, rn);

                CRay *ray1 = new CRay;

                if(i < show_rays) {

                        detector->Propagate(*ray0, ray1, show_3d);

                        }

                else {

                        detector->Propagate(*ray0, ray1, 0);

                        }

          //delete ray0;

          //delete ray1;        

        }

        if(show_rand) {

                TCanvas *c2rand = (TCanvas*)gROOT->FindObject("c2rand");

                if(!c2rand) c2rand = new TCanvas("c2rand", "Random", 750, 550, 700, 700);

                else c2rand->Clear();

                c2rand->Divide(2,3);

                c2rand->cd(1); hphi->Draw();

                c2rand->cd(3); htheta->Draw();

                c2rand->cd(5); hcostheta->Draw();

                c2rand->cd(2); hl->Draw();

                c2rand->cd(4); hm->Draw();

                c2rand->cd(6); hn->Draw();

        }

        double acceptance = (detector->GetHActive())->GetEntries() / (double)NN;       

        return acceptance;

}