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