Subversion Repositories f9daq

#include "RTUtil.h"

#include "CGuide.h"

#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"

int show_3d = 1;

int show_data = 1;

int draw_width = 2;

TVector3 center(-2.0, 0.0, 0.0);

double RandIso(int NN = 1e3, double theta = 30.0, int show_rays = 30, int show_rand = 0)

{

  CDetector *detector = new CDetector(center);

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

        //if (show_rand) {

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

        int fate;

        double SiPM = detector->GetSiPM();

        double M = detector->GetM();

        CRay *ray0 = new CRay;

        CRay *ray1 = new CRay;

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

                rx = center.x();

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

                }

                ray0->Set(rx, ry, rz, rl, rm, rn);

                if(i < show_rays)

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

                else

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

        }

        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();

                delete c2rand;

        }

        delete hphi;

        delete htheta;

        delete hcostheta;

        delete hl;

        delete hm;

        delete hn;

        delete ray0;

        delete ray1;

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

}

// a vs. d

void LGI_ad(int NN = 1e4, double min = 2.5, double max = 3.5, double minD = 1, double maxD = 6, const int steps = 10, double theta = 30.0)

{

        show_3d = 0; // don't show every simulation 

        show_data = 1;

        //double d  = detector->GetD();

        double b = 5.0; // upper side of LG

        double R = 0.96;

        TH2F *hAcceptance = new TH2F("hAcceptance","Acceptance",steps-1,minD,maxD,steps-1,min,max);

        printf("   d |   a  | Acceptance\n");

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

                double x = hAcceptance->GetXaxis()->GetBinCenter(i);

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

                  double y = hAcceptance->GetYaxis()->GetBinCenter(j);

                  double M = b/y;

                  // hardcoded n0 and n1

                  SetLG(y, M, x);

                  //Init(); exclude simulation

                  double acceptance = RandIso(NN, theta, 5, 0);

                  //double acceptance = Grid(NN, theta);

                  //double acceptance = -1.0;

                  hAcceptance->Fill(x, y, acceptance);

                  //printf("%.2lf | %.2lf | ", x,y);

                  //PrintGuideStat(NN);

                  }

        }

        TCanvas *cp = new TCanvas("Acceptance", "Acceptance", 0, 0, 640, 480);

        TVirtualPad *pacc = cp->cd(0);

        pacc->SetRightMargin(0.10);

        pacc->SetLeftMargin(0.10);

        pacc->SetTopMargin(0.10);

        //TFile *file = new TFile("acceptance.root","RECREATE");

        //hAcceptance->Write();

        //file->Close();

        //hAcceptance->SetContour(100);

        //gStyle->SetPalette(1,0);

        hAcceptance->SetTitle(";d [mm];a [mm]");

        hAcceptance->GetXaxis()->SetRangeUser(minD,maxD);

        hAcceptance->Draw("COLZ");

 return;

}