Subversion Repositories f9daq

Rev

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

// 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, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
                                parameters.getA(), parameters.getLightYield()*acc, 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");
                        (detector->GetHGlass())->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;
        if(phi < 1e-6) phi = 1e-6;
       
        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);
       
        delete ray0;
        delete ray1;
       
        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);
                delete ray0;
                delete ray1;   
                }
               
        }
        double acceptance = 0.0;
        /*
        if( !(parameters.getPlateOn()) )
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN*NN;
  else
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
    */

        acceptance = (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
        return acceptance;
}
//-----------------------------------------------------------------------------
// zarki z nakljucnim polozajem vpada (na vstopni pov. vodnika)
// vsi pod kotom (theta, phi)
double RandYZ(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phi, int show_rays)
{
//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; 
        }
       
        double acceptance = 0.0;
        /*
        if( !(parameters.getPlateOn()) )
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  else
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();       
        //return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));*/

        acceptance = (detector->GetHActive())->GetEntries() / (double) NN;
        return acceptance;     
}
//-----------------------------------------------------------------------------
// zarki, izotropno porazdeljeni znotraj kota theta
// = nakljucni vstopni polozaj in kot
// NN - number of rays to be simulated with angles theta distributed uniformly
//      in cos theta and phi uniformly:
//      theta [0, theta]
//      phi [0,360]
double RandIso(CDetector *detector, DetectorParameters& parameters, int NN = 1e3, double thetaMin=0.0, 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);
        double theta_max_rad = TMath::Cos(thetaMin);

          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);
                //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
                //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
                rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
               
                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 = 0.0;
        /*
        if( !(parameters.getPlateOn()) )
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  else
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
    */

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

// Beamtest distribution
// with fixed theta and phi in interval phiMin, phiMax
double beamtest(CDetector *detector, DetectorParameters& parameters, int NN, double theta, double phiMin, double phiMax, int show_rays, int show_rand)
{
        double pi = 3.14159265358979312;
        theta *= pi/180.0;
        if(theta < MARGIN) theta = MARGIN;
        phiMin *= pi/180.0;
        phiMax *= pi/180.0;
       
        TDatime now;
        TRandom rand;
        rand.SetSeed(now.Get());
        double rx, ry, rz, rl, rm, rn;
        double rphi;
        //double 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, -pi, 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 b = parameters.getB();
        double offset = (parameters._plateOn ? parameters._plateWidth : 0);
       
        for(int i = 0; i < NN; i++) {
               
                rx = CENTER.x() - offset;
                ry = rand.Uniform(-b/2.0, b/2.0);
                rz = rand.Uniform(-b/2.0, b/2.0);
                               
                rphi = rand.Uniform(phiMin, phiMax);
                //double theta_max_rad = TMath::Cos(17.0*pi/180.0);
                //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, 1.0)));
                //rtheta = TMath::ACos((rand.Uniform(theta_min_rad, theta_max_rad)));
               
                rl = TMath::Cos(theta);
                rm = TMath::Sin(theta)*TMath::Sin(rphi);
                rn = TMath::Sin(theta)*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 = 0.0;
        /*
        if( !(parameters.getPlateOn()) )
          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
  else
    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
    */

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