#include <RAT/DU/DSReader.hh>
#include <RAT/DU/Utility.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/MC.hh>
#include <RAT/DS/EV.hh>
#include <RAT/DS/PMT.hh>
#include <RAT/DU/PMTInfo.hh>
#include <vector>
#include "TH2D.h"
#include "TH1D.h"
#include <TVector3.h>
#include <TMath.h>
#include <TROOT.h>
#include "TFile.h"
using namespace std ;
void WL()
{
RAT::DU::DSReader dsReader("WLS_beta_turnoffWater_100evts_track.root");
const RAT::DU::PMTInfo& pmtInfo = RAT::DU::Utility::Get()->GetPMTInfo();
TH2D* PMTIDHitTime = new TH2D("PMTIDHitTime", "PMTID vs Hit Time", 10000,1,10000,600,0,600);
PMTIDHitTime->GetXaxis()->SetTitle("PMTID");
PMTIDHitTime->GetYaxis()->SetTitle("hitTime [nsecs]");
TH1D *WavelengthEmission = new TH1D("WavelengthEmission","Wavelength which re-emitted after absorption", 600,0,600);
WavelengthEmission->GetXaxis()->SetTitle("Wavelength [nm]");
TH1D *WavelengthAbsorption = new TH1D("WavelengthAbsorption","Wavelength which absorbed", 600,0,600);
WavelengthAbsorption->GetXaxis()->SetTitle("Wavelength [nm]");
TH2D *hYZxAbsPosbin2 = new TH2D("Abs_yzXpos_uniXrange","x",1000, 0, 6000, 1000, 0, 6000);
TH2D *hYZxReEmitPosbin2 = new TH2D("ReEmit_yzXpos_uniXrange","x",1000, 0, 6000, 1000, 0, 6000);
TH1D *hAbsDistance = new TH1D("Abs_Distance","",1000,0,6000);
TH1D *hReEmitDistance = new TH1D("ReEmit_Distance","",1000,0,6000);
TH1D *hAbsCosTheta = new TH1D("Abs_cosTheta","",1000,0.0,1.0);
TH1D *hReEmitCosTheta = new TH1D("ReEmit_cosTheta","",1000,0.0,1.0);
TH2D *DistanceCosTheta = new TH2D("ReEmitDistance-cosTheta","CosTheta",1000,0,6000,1000,0,1.0);
TVector3 PMTVector, calpmtVector, MomentumVec, eventPosition, momentumCerenkov, posCherenkov ;
eventPosition.SetXYZ(0,0,0);
double grVelocity = 2.17554021555098529e+02 ;
double eventTime = 0.0 ;
double KECerenkov = 0.0 ;
double KEWLS = 0.0 ;
double cerenkovangle = 0.0 ;
double CerenkovWavelength = 0.0 ;
double WLSShiftedWavelength =0.0 ;
string process1 = "Scintillation" ;
string process2 = "Cerenkov" ;
string process3 = "OpAbsorption" ;
string process4 = "Reemission_from_comp1" ;
double energy, wavelength ;
UInt_t UTrackID = 999999 ;
Int_t nPhotons_r =0;
Int_t nPhotons_c =0;
for( size_t iEntry = 0; iEntry<100; iEntry++) //dsReader.GetEntryCount(); iEntry++ )
{
//std:cout << " event ID "<< iEntry <<std::endl ;
const RAT::DS::Entry& rDS = dsReader.GetEntry( iEntry );
const RAT::DS::MC& rmc= rDS.GetMC();
nPhotons_r =nPhotons_r+ rmc.GetNReemittedPhotons() ;
nPhotons_c = nPhotons_c + rmc.GetNCherPhotons() ;
if(rmc.GetMCParticle(0).GetPDGCode()==11)
MomentumVec.SetXYZ(rmc.GetMCParticle(0).GetMomentum().x(),rmc.GetMCParticle(0).GetMomentum().y(),rmc.GetMCParticle(0).GetMomentum().z()) ;
TVector3 startPos_track2, endPos_track2, startPos_track1, endPos_track1;
double startTime_track2 =0.0, endTime_track2=0.0, startTime_track1= 0.0, endTime_track1=0.0;
for(size_t iTrack=0; iTrack< rmc.GetMCTrackCount(); iTrack++)
{
if(rmc.GetMCTrack(iTrack).GetLastMCTrackStep().GetProcess()==process3 && rmc.GetMCTrack(iTrack).GetFirstMCTrackStep().GetProcess()==process2 && rmc.GetMCTrack(iTrack).GetLastMCTrackStep().GetEndVolume()=="inner_av")
{
startPos_track1.SetXYZ(rmc.GetMCTrack(iTrack).GetFirstMCTrackStep().GetPosition().x(),rmc.GetMCTrack(iTrack).GetFirstMCTrackStep().GetPosition().y(),rmc.GetMCTrack(iTrack).GetFirstMCTrackStep().GetPosition().z()) ;
endPos_track1.SetXYZ(rmc.GetMCTrack(iTrack).GetLastMCTrackStep().GetPosition().x(),rmc.GetMCTrack(iTrack).GetLastMCTrackStep().GetPosition().y(),rmc.GetMCTrack(iTrack).GetLastMCTrackStep().GetPosition().z());
KECerenkov = rmc.GetMCTrack(iTrack).GetFirstMCTrackStep().GetKineticEnergy();
CerenkovWavelength = ((6.625e-34)*(3e+8)/(KECerenkov*1e+6*1.6e-19))*1e+9;
WavelengthAbsorption->Fill(CerenkovWavelength);
hYZxAbsPosbin2->Fill(TMath::Sqrt(pow(endPos_track1.Z(),2)+pow(endPos_track1.Y(),2)),endPos_track1.X());
hAbsDistance->Fill((endPos_track1-startPos_track1).Mag());
hAbsCosTheta->Fill(endPos_track1.X()/endPos_track1.Mag());
}
}
for(size_t iTr=0; iTr< rmc.GetMCTrackCount(); iTr++)
{
if(rmc.GetMCTrack(iTr).GetFirstMCTrackStep().GetProcess()==process4)
{
startPos_track2.SetXYZ(rmc.GetMCTrack(iTr).GetFirstMCTrackStep().GetPosition().x(),rmc.GetMCTrack(iTr).GetFirstMCTrackStep().GetPosition().y(),rmc.GetMCTrack(iTr).GetFirstMCTrackStep().GetPosition().z());
endPos_track2.SetXYZ(rmc.GetMCTrack(iTr).GetLastMCTrackStep().GetPosition().x(),rmc.GetMCTrack(iTr).GetLastMCTrackStep().GetPosition().y(),rmc.GetMCTrack(iTr).GetLastMCTrackStep().GetPosition().z());
KEWLS = rmc.GetMCTrack(iTr).GetFirstMCTrackStep().GetKineticEnergy() ;
WLSShiftedWavelength = ((6.625e-34)*(3e+8)/(KEWLS*1e+6*1.6e-19))*1e+9 ;
WavelengthEmission->Fill(WLSShiftedWavelength);
hYZxReEmitPosbin2->Fill(TMath::Sqrt(pow(startPos_track2.Z(),2)+pow(startPos_track2.Y(),2)),startPos_track2.X());
hReEmitDistance->Fill((startPos_track1-startPos_track2).Mag());
hReEmitCosTheta->Fill(startPos_track2.X()/startPos_track2.Mag());
DistanceCosTheta->Fill((startPos_track1-startPos_track2).Mag(),startPos_track2.X()/startPos_track2.Mag());
}
}
}
double ReEmitArea=WavelengthEmission->Integral();
double AbsArea=WavelengthAbsorption->Integral();
cout<<"ratio:"<<ReEmitArea/AbsArea<<endl;
TFile *f=new TFile("chgWater_OPTICS_5MeV_processCheck.root","RECREATE");
f->cd();
WavelengthEmission->Write();
WavelengthAbsorption->Write();
hYZxAbsPosbin2->Write();
hYZxReEmitPosbin2->Write();
hAbsDistance->Write();
hAbsCosTheta->Write();
hReEmitDistance->Write();
hReEmitCosTheta->Write();
DistanceCosTheta->Write();
f->Close();
std::cout<<"Cerenkov photons "<< nPhotons_c << " Remitted photons " << nPhotons_r <<std::endl ;
}