FDsimG4Write2ROOT.cc

Go to the documentation of this file.

#include "FDsimG4Defs.hh"
#ifdef __TELESCOPE_SIMULATOR_LX_USE_ROOT__

#include "FDsimG4Write2ROOT.hh"
#include "FDsimG4StoreOpticalHit.hh"
#include <TFile.h>
#include <TTree.h>
#include <TDirectory.h>
#include <utl/Photon.h>
#include <utl/AugerUnits.h>
#include <iostream>
#include <string>

namespace TelescopeSimulatorLX{double RCam = 1.743*utl::m;}

using namespace std;
using namespace TelescopeSimulatorLX;

FDsimG4Write2ROOT::FDsimG4Write2ROOT() : 
  fRootFile(NULL),
  fTreePhotons(NULL)
{
}


FDsimG4Write2ROOT::FDsimG4Write2ROOT(const string& RootFilename)
{
  Open(RootFilename);
  CreateTree();
}


FDsimG4Write2ROOT::~FDsimG4Write2ROOT()
{
  SaveAndClose();
}


void FDsimG4Write2ROOT::Open(const string& filename)
{
  TDirectory *tmpSaveDir = gDirectory;
  fRootFile = new TFile(filename.c_str(), "RECREATE");
  fSaveDir=gDirectory;
  tmpSaveDir->cd();
}


void FDsimG4Write2ROOT::SaveAndClose()
{
  if (fRootFile) {
    if (fRootFile->IsOpen()) {
      TDirectory *tmpSaveDir = gDirectory;
      fSaveDir->cd();
      fTreePhotons->Write("",TObject::kOverwrite);
      fRootFile->Close();
      tmpSaveDir->cd(); 
      fTreePhotons  = NULL;
    }
  }

  if (fRootFile) {
    delete fRootFile; 
    fRootFile = NULL;
  }
}


void FDsimG4Write2ROOT::CreateTree()
{
  TDirectory *tmpSaveDir = gDirectory;
  fSaveDir->cd();
  if (!fRootFile->Get("Photons")){
    fTreePhotons = new TTree("Photons","Photons Tree");
    
    fTreePhotons->Branch("G4",
                         &fG4_info,
                         "x/D:y/D:z/D:alfa/D:beta/D:E/D:w/D:t/D:LocX/D:LocY/D:LocZ/D:PixX/D:PixY/D:PixZ/D:HitID/D:PMTid/D");
    
    fTreePhotons->Branch("G4Filter",&fG4Filter_info,"x/D:y/D:z/D");
    fTreePhotons->Branch("G4Lens",&fG4Lens_info,"x/D:y/D:z/D");
    fTreePhotons->Branch("G4Mirror",&fG4Mirror_info,"x/D:y/D:z/D");
    fTreePhotons->Branch("G4Mercedes",&fG4Mercedes_info,"x/D:y/D:z/D");
    fTreePhotons->Branch("G4FS",&fG4FS_info,"x/D:y/D:z/D:alfa/D:beta/D:th/D:phi/D");
    fTreePhotons->Branch("nHit",
                         &fn_info,
                         "G4Filter/I:G4Lens/I:G4Mirror/I:G4Mercedes/I:G4FS/I:G4/I");
    
    fTreePhotons->Branch("IN",
                         &fIN_info,
                         "x/D:y/D:z/D:DX/D:DY/D:DZ/D:th/D:phi/D:wantedPix/D:wantedTh/D:wantedPhi/D");
    fTreePhotons->Branch("INExp",
                         &fINExp_info,
                         "x/D:y/D:z/D:alfa/D:beta/D:DX/D:DY/D:DZ/D:th/D:phi/D");
  }
  else 
    {
      fTreePhotons = (TTree*)fRootFile->Get("Photons");
      
      fTreePhotons->SetBranchAddress("G4",&fG4_info);
      fTreePhotons->SetBranchAddress("G4Filter",&fG4Filter_info);
      fTreePhotons->SetBranchAddress("G4Lens",&fG4Lens_info);
      fTreePhotons->SetBranchAddress("G4Mirror",&fG4Mirror_info);
      fTreePhotons->SetBranchAddress("G4Mercedes",&fG4Mercedes_info);
      fTreePhotons->SetBranchAddress("G4FS",&fG4FS_info);
      fTreePhotons->SetBranchAddress("nHit",&fn_info);
      fTreePhotons->SetBranchAddress("IN",&fIN_info);
      fTreePhotons->SetBranchAddress("INExp",&fINExp_info);
    }
  tmpSaveDir->cd();
}


void FDsimG4Write2ROOT::SetG4FilterPhoton(FDsimG4StoreOpticalHit* G4_photon)
{
  fG4Filter_info.x = G4_photon->GetPosX()*utl::mm;
  fG4Filter_info.y = G4_photon->GetPosY()*utl::mm;
  fG4Filter_info.z = G4_photon->GetPosZ()*utl::mm;
  fn_info.G4Filter = 1;
}


void FDsimG4Write2ROOT::SetG4LensPhoton(FDsimG4StoreOpticalHit* G4_photon)
{
  fG4Lens_info.x = G4_photon->GetPosX()*utl::mm;
  fG4Lens_info.y = G4_photon->GetPosY()*utl::mm;
  fG4Lens_info.z = G4_photon->GetPosZ()*utl::mm;
  fn_info.G4Lens = 1;
}


void FDsimG4Write2ROOT::SetG4MirrorPhoton(FDsimG4StoreOpticalHit* G4_photon)
{
  fG4Mirror_info.x = G4_photon->GetPosX()*utl::mm;
  fG4Mirror_info.y = G4_photon->GetPosY()*utl::mm;
  fG4Mirror_info.z = G4_photon->GetPosZ()*utl::mm;
  fn_info.G4Mirror = 1;
}


void FDsimG4Write2ROOT::SetG4MercedesPhoton(FDsimG4StoreOpticalHit* G4_photon)
{
  fG4Mercedes_info.x = G4_photon->GetPosX()*utl::mm;
  fG4Mercedes_info.y = G4_photon->GetPosY()*utl::mm;
  fG4Mercedes_info.z = G4_photon->GetPosZ()*utl::mm;
  fn_info.G4Mercedes = 1;
}


void FDsimG4Write2ROOT::SetG4FSPhoton(FDsimG4StoreOpticalHit* G4_photon)
{
  fG4FS_info.x = G4_photon->GetPosX()*utl::mm;
  fG4FS_info.y = G4_photon->GetPosY()*utl::mm;
  fG4FS_info.z = G4_photon->GetPosZ()*utl::mm;
//calcula  

  fG4FS_info.alfa = asin(fG4FS_info.x/RCam);
  fG4FS_info.beta = asin(fG4FS_info.y/-RCam);
 
  utl::Vector tempVec= utl::Vector(fG4FS_info.x,fG4FS_info.y,fG4FS_info.z,ftelCS);
  tempVec*=-1;
   //debugpj//cout << tempVec.GetTheta(ftelCS)/3.1416*180. <<" telG4 "<<tempVec.GetPhi(ftelCS)/3.1416*180.<<" tel "<<tempVec.GetMag() <<endl;
   //debugpj//cout << tempVec.GetTheta(feyeCS)/3.1416*180. <<"  G$ "<<tempVec.GetPhi(feyeCS)/3.1416*180.<<endl;

//Changed for Valerio//  
fG4FS_info.th = tempVec.GetTheta(feyeCS);
//Changed for Valerio//  
fG4FS_info.phi = tempVec.GetPhi(feyeCS);
//Changed for Valerio
//  fG4FS_info.th = tempVec.GetTheta(ftelCS);
//  fG4FS_info.phi = tempVec.GetPhi(ftelCS);


  fn_info.G4FS = 1;
}


void FDsimG4Write2ROOT::SetG4Photon(FDsimG4StoreOpticalHit* G4_photon)
{
  fG4_info.x = G4_photon->GetPosX()*utl::mm;
  fG4_info.y = G4_photon->GetPosY()*utl::mm;
  fG4_info.z = G4_photon->GetPosZ()*utl::mm;
  fn_info.G4 = 1;
//calcula
  fG4_info.alfa = asin(fG4_info.x/RCam);
  fG4_info.beta = asin(fG4_info.y/-RCam);
// Other info

  fG4_info.E = G4_photon->GetEnergy()*utl::eV;        
  fG4_info.w= G4_photon->GetWavelength()*utl::nanometer;
  fG4_info.t = G4_photon->GetTime()*utl::ns;
  
  fG4_info.LocX = G4_photon->GetLocalPosX()*utl::mm;   
  fG4_info.LocY = G4_photon->GetLocalPosY()*utl::mm;     
  fG4_info.LocZ = G4_photon->GetLocalPosZ()*utl::mm;     
  
  fG4_info.PixX = G4_photon->GetPixelPosX()*utl::mm;     
  fG4_info.PixY = G4_photon->GetPixelPosY()*utl::mm;     
  fG4_info.PixZ = G4_photon->GetPixelPosZ()*utl::mm;     
  
  fG4_info.HitID = G4_photon->GetHitID();      
  fG4_info.PMTid = G4_photon->GetPMTid();        
}

void FDsimG4Write2ROOT::SetINPhoton(FDsimG4StoreOpticalHit* G4_photon)
{  

  double X = G4_photon->GetPosX()*utl::mm;
  double Y = G4_photon->GetPosY()*utl::mm;
  double Z = G4_photon->GetPosZ()*utl::mm;

  double DX = G4_photon->GetDirX();
  double DY = G4_photon->GetDirY();
  double DZ = G4_photon->GetDirZ();

  double lambda = G4_photon->GetWavelength()*utl::nanometer;
  double weight = G4_photon->GetWeight();

  utl::Point Pos = utl::Point(X,Y,Z,ftelCS);
  utl::Vector Dir = utl::Vector(DX,DY,DZ,ftelCS);

  utl::Photon IN_photon(Pos,Dir,lambda,weight);

  SetINPhoton(IN_photon,ftelCS,feyeCS);
}


void 
FDsimG4Write2ROOT::SetINPhoton(utl::Photon IN_photon , 
                               utl::CoordinateSystemPtr telCS,
                               utl::CoordinateSystemPtr eyeCS)
{
  utl::CoordinateSystemPtr CS=eyeCS;
  static double RCam=1.743*utl::m;
  fIN_info.x   = IN_photon.GetPosition().GetX(telCS);
  fIN_info.y   = IN_photon.GetPosition().GetY(telCS);
  fIN_info.z   = IN_photon.GetPosition().GetZ(telCS);
  fIN_info.DX  = IN_photon.GetDirection().GetX(telCS);
  fIN_info.DY = IN_photon.GetDirection().GetY(telCS); 
  fIN_info.DZ  = IN_photon.GetDirection().GetZ(telCS);

//atencao que isto aponta para dentro do Tel
  fIN_info.th  = IN_photon.GetDirection().GetTheta(telCS);
  fIN_info.phi = IN_photon.GetDirection().GetPhi(telCS);

  //temporary container for thphi
  fINExp_info.th=IN_photon.GetDirection().GetTheta(telCS);
  fINExp_info.phi=IN_photon.GetDirection().GetPhi(telCS);

  fINExp_info.x = RCam*sin(fINExp_info.th)*cos(fINExp_info.phi);
  fINExp_info.y = RCam*sin(fINExp_info.th)*sin(fINExp_info.phi);
  fINExp_info.z = -RCam*cos(fINExp_info.th);
  fINExp_info.alfa = asin(fINExp_info.x/RCam);
  fINExp_info.beta = asin(fINExp_info.y/-RCam);

//   utl::Vector tempVec= utl::Vector(fINExp_info.x,fINExp_info.y,-fINExp_info.z,ftelCS);
//   tempVec*=-1;
//   fINExp_info.th = tempVec.GetTheta(feyeCS);
//   fINExp_info.phi = tempVec.GetPhi(feyeCS);
}


void FDsimG4Write2ROOT::FillTree()
{
  WriteIt();
}


void FDsimG4Write2ROOT::WriteIt()
{
    fTreePhotons->Fill();
}


void FDsimG4Write2ROOT::Clean()
{
  fG4Filter_info.x = nan("");
  fG4Filter_info.y = nan("");
  fG4Filter_info.z = nan("");

  fG4Lens_info.x = nan("");
  fG4Lens_info.y = nan("");
  fG4Lens_info.z = nan("");

  fG4Mirror_info.x = nan("");
  fG4Mirror_info.y = nan("");
  fG4Mirror_info.z = nan("");

  fG4Mercedes_info.x = nan("");
  fG4Mercedes_info.y = nan("");
  fG4Mercedes_info.z = nan("");

  fG4FS_info.x = nan("");
  fG4FS_info.y = nan("");
  fG4FS_info.z = nan("");
  fG4FS_info.alfa = nan("");
  fG4FS_info.beta = nan("");
  fG4FS_info.th=nan("");
  fG4FS_info.phi=nan("");

  fG4_info.x = nan("");
  fG4_info.y = nan("");
  fG4_info.z = nan("");
  fG4_info.E = nan("");
  fG4_info.w = nan("");
  fG4_info.t = nan("");
  fG4_info.LocX = nan("");
  fG4_info.LocY = nan("");
  fG4_info.LocZ = nan("");
  fG4_info.PixX = nan("");
  fG4_info.PixY = nan("");
  fG4_info.PixZ = nan("");
  fG4_info.HitID = nan("");
  fG4_info.PMTid = nan("");

  fIN_info.x= nan("");
  fIN_info.y= nan("");
  fIN_info.z= nan("");
  fIN_info.DX= nan("");
  fIN_info.DY= nan("");
  fIN_info.DZ= nan("");
  fIN_info.th= nan("");
  fIN_info.phi= nan("");

  fINExp_info.x= nan("");
  fINExp_info.y= nan("");
  fINExp_info.z= nan("");
  fINExp_info.alfa= nan("");
  fINExp_info.beta= nan("");
  fINExp_info.DX= nan("");
  fINExp_info.DY= nan("");
  fINExp_info.DZ= nan("");
  fINExp_info.th= nan("");
  fINExp_info.phi= nan("");

  fn_info.G4Filter = 0;
  fn_info.G4Lens = 0;
  fn_info.G4Mirror = 0;
  fn_info.G4Mercedes = 0;
  fn_info.G4FS =0;
  fn_info.G4 = 0;
}

#endif