TestTorus.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "TestTorus.h"
#include "RTFunctions.h"

#include <TCanvas.h>
#include <TPolyMarker3D.h>

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/Math.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Point.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/Trace.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/Photon.h>
#include <utl/RandomEngine.h>
#include <utl/UTMPoint.h>

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/RandomEngineRegistry.h>

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/Camera.h>
#include <fdet/Pixel.h>
#include <fdet/Mirror.h>
#include <fdet/Filter.h>
#include <fdet/Corrector.h>

#include <evt/Event.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/Telescope.h>
#include <fevt/PixelSimData.h>
#include <fevt/Pixel.h>

#include <boost/tuple/tuple.hpp>

#include <TDirectory.h>
#include <TFile.h>
#include <TTree.h>

#include <CLHEP/Random/Randomize.h>

#include <iomanip>
#include <fstream>
#include <sstream>

using namespace TestTorusKG;
using namespace std;
using namespace utl;
using namespace evt;
using namespace fwk;
using namespace det;
using namespace fdet;

using CLHEP::RandFlat;



TestTorus::TestTorus() 
{
}


VModule::ResultFlag 
TestTorus::Init()
{
  cout << " Init TestTorus " << endl;
  
  return eSuccess;
}


VModule::ResultFlag
TestTorus::Run(evt::Event& /*event*/)
{
  cout << " Run TestTorus " << endl;

  TimeStamp simTime(900000000, 3);   
  
  Detector& detector = Detector::GetInstance();
  detector.Update (simTime);
  const FDetector& theFDetector = detector.GetFDetector();
  
  const fdet::Eye &detEye = *(theFDetector.EyesBegin ());
  const fdet::Telescope &detTel = *(detEye.TelescopesBegin ());
  const CoordinateSystemPtr& telCS = detTel.GetTelescopeCoordinateSystem ();
  
  utl::RandomEngine* theRandomEngine =
    &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);

  const double tubeRadius = 524351.7*mm;//6130*mm;//8380*mm;
  const double torusRadius = 750.*mm;//850.*mm;//804.1*mm;
  const double torusZ = tubeRadius + 9.39*mm;//tubeRadius + 6.9*mm;//8386.775*mm;

  const double phiMin = 0;
  const double phiMax = 360*deg;
  const double distMin = 0.*m;
  const double distMax = 1.1*(tubeRadius+torusRadius);
  
  //const double thetaMin = 0;
  //const double thetaMax = 30*deg;

  const Point origin(0,0,torusZ ,telCS);
  const Vector inwards(0,0,-1,telCS);

  TCanvas* canvas = new TCanvas("torus");

  const int n = 100000;
  int count = 0;
  TPolyMarker3D* poly = new TPolyMarker3D(n*4);


  for (int i=0; i<n; ++i) {
       
    const double phi = RandFlat::shoot(&theRandomEngine->GetEngine(), phiMin, phiMax);
    const double dist =  RandFlat::shoot(&theRandomEngine->GetEngine(), distMin, distMax);
    
    const utl::Point point(dist*cos(phi), dist*sin(phi), 3.*tubeRadius, telCS);

    //const double phi2 = RandFlat::shoot(&theRandomEngine->GetEngine(), phiMin, phiMax);
    //const double theta = RandFlat::shoot(&theRandomEngine->GetEngine(), thetaMin, thetaMax);
    
    //const utl::Vector dir(sin(theta)*cos(phi2), sin(theta)*sin(phi2), -cos(theta), telCS);
    const utl::Vector dir(0, 0, -1., telCS);
    
    Photon phot(point, dir, 1, 1, 1);
    
    TelescopeSimulatorKG2::RTFunctions::IntersectionList intersections;
    TelescopeSimulatorKG2::RTFunctions::Torus(origin, inwards, torusRadius, tubeRadius, phot, intersections);
    
    if (intersections.size()==0)
      continue;

    if (intersections.size() != 2 && intersections.size() != 4 ) {
      cout << " s=" << intersections.size() << endl;
    }

    unsigned int indexOfSmallestZ = 4; // maximum number of intersections is 4 (0 to 3)
    double smallestZ = 10*tubeRadius/m; // definitely outside the torus and outside the telescope
    for (unsigned int j=0; j<intersections.size(); ++j) {
      
      const double tmp_x = intersections[j].first.GetPosition().GetX(telCS);
      const double tmp_y = intersections[j].first.GetPosition().GetY(telCS);
      const double tmp_r = sqrt(tmp_x*tmp_x + tmp_y*tmp_y);
      const double tmp_z = intersections[j].first.GetPosition().GetZ(telCS);

      if( tmp_r/m < 1.1*m && tmp_r/m > 0.85*m && tmp_z/m < smallestZ/m ) {
        smallestZ = tmp_z/m;
        indexOfSmallestZ = j;
      }
      poly->SetPoint(count++, tmp_x, tmp_y, tmp_z);

    }
    if (indexOfSmallestZ < 4) {
      const double tmpx = intersections[indexOfSmallestZ].first.GetPosition().GetX(telCS);
      const double tmpy = intersections[indexOfSmallestZ].first.GetPosition().GetY(telCS);
      const double tmpz = intersections[indexOfSmallestZ].first.GetPosition().GetZ(telCS);
      const double tmpr = sqrt(tmpx*tmpx + tmpy*tmpy);

      cout << "lenspoint " << tmpx/mm << "  "
                           << tmpy/mm << "  "
                           << tmpz/mm << "  "
                           << tmpr/mm << endl;
     /* poly->SetPoint(count++,
         intersections[indexOfSmallestZ].first.GetPosition().GetX(telCS),
         intersections[indexOfSmallestZ].first.GetPosition().GetY(telCS),
         intersections[indexOfSmallestZ].first.GetPosition().GetZ(telCS));
         */
    }
    
  }
  
  canvas->cd();
  poly->Draw();

  canvas->SaveAs("torus.C");
  canvas->SaveAs("torus.eps");

  return eSuccess;
} // end of Run


VModule::ResultFlag
TestTorus::Finish()
{
  cout << " Finish TestTorus " << endl;

  return eSuccess;
}

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "cd $AUGER_BASE && make"
// End: