WebSVN – f9daq – Diff – /lightguide/trunk/src/raySimulator.cpp

 //-----------------------------------------------------------------------------
 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",
+                sprintf(sbuff, "SiPM = %.1lf, L.y. = %.2lf, d = %.1lf | gap = (%.1lf, %.1lf, %.1lf), #theta = %.1lf | acceptance = %.3lf",
-                                parameters.getA(), parameters.getM(), parameters.getD(),
+                                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);
+{
   //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
                                           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();
+    */
-        return (detector->GetHActive())->GetEntries() / (double)((NN)*(NN));
+        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 = 1e3, double theta = 0.0, double phi = 0.0, int show_rays = 30)
+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;
                 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));*/
-        return (detector->GetHActive())->GetEntries() / (double)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 theta = 30.0, int show_rays = 30, int show_rand = 0)
+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;
                 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, 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 ray0;
-          //delete ray1;
+          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(2); hl->Draw();
                 c2rand->cd(4); hm->Draw();
                 c2rand->cd(6); hn->Draw();
+        }
+        double acceptance = 0.0;
+        /*
+        if( !(parameters.getPlateOn()) )
-        double acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
+          acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
+  else
+    acceptance = (detector->GetHActive())->GetEntries() / (detector->GetHPlate())->GetEntries();
+    */
+  acceptance = (detector->GetHActive())->GetEntries() / (double)NN;
         return acceptance;
+}

Subversion Repositories f9daq

(root)/lightguide/trunk/src/raySimulator.cpp – Rev 54 → 70