Skip to content
Snippets Groups Projects
Commit 16bd3582 authored by shengy3's avatar shengy3
Browse files

add batchWaveformFolder

parent e62b8894
No related branches found
No related tags found
1 merge request!1add batchWaveformFolder
Hide whitespace changes
Inline Side-by-side
Showing
with 525 additions and 514 deletions
outputDigitization.cc
+ 525
− 514
View file @ 16bd3582
#include "TFile.h"
#include "TH1F.h"
#include "TGraph.h"
#include "TF1.h"
#include "TRandom.h"
#include <TTree.h>
#include <vector>
#include <fstream>
using namespace std;
TGraph* getPMTdata( string fileName );
double evaluateSum( vector< TF1* >& funcVec, double x );
void generateWaveform(string fileName, bool DRAW = false);
void addDarkCurrent(vector< TF1* >& vec, float timeInterval, float timeBinWidth);
void addRandNoise(vector< float >* vec, float timeInterval);
void addTimingJitter(vector< float >* vec, float timeInterval);
Double_t PMTpulseFunction(Double_t *t, Double_t *par){
Float_t tt = t[0];
if(tt < par[1]) return 0.0; // cut out the time where the function behaves poorly
Double_t f = par[0]*pow((tt-par[1])/10.0, 3.4)*exp(-(tt-par[1])/10.0);
return f;
}
void batchWaveform(){
int fileStart = 0;
int fileStop = 100;
for(int fileNum = fileStart; fileNum < fileStop; fileNum++){
generateWaveform( Form("filename%d.root",fileNum) );
}
}
/*
* Calculate the number of Cherenkov photons likely to make it out
* of the ZDC starting from the output of the JZCaPA MC. Accounting
* for losses due to absorption in fused silica after heavy irradiation.
*/
void cutZDCoutput(string fileName, string rod, double xOffset, double zOffset){
int nBins = 3649;
double energy,X,Z,Px,Py,Pz;
TRandom rnd;
//Create input vectors
vector<double> *Xf=0, *Yf=0, *Zf=0, *Pxf=0, *Pyf=0, *Pzf=0, *Energyf=0, *NCherenkovs=0;
vector<int> *PIDf=0, *IDf=0, *EventNof=0;
//Create output vectors
vector<double> *x=0, *z=0,*x_0=0, *z_0=0, *px=0, *py=0, *pz=0, *Energy=0;
vector<int> *EventNo=0;
// Set the input and output files up
//
TFile* inFile = new TFile(fileName.c_str());
TTree* inTree = (TTree*)inFile->Get("ZDCtree");
inTree->SetBranchAddress("X",&Xf);
inTree->SetBranchAddress("Z",&Zf);
inTree->SetBranchAddress("Px",&Pxf);
inTree->SetBranchAddress("Py",&Pyf);
inTree->SetBranchAddress("Pz",&Pzf);
//inTree->SetBranchAddress("EventNo",&EventNf);
// intree->SetBranchAddress("NCherenkovs",&NCherenkovs);
//Output root file
TFile* outFile = new TFile( Form("output/%s_%s.root",fileName.substr(0,fileName.find_last_of(".")).c_str(), rod.c_str() ), "RECREATE");
TTree* outTree = new TTree( "tree", "tree" );
// outTree->Branch("X0",&x_0);
// outTree->Branch("Z0",&z_0);
outTree->Branch("X",&x);
outTree->Branch("Z",&z);
outTree->Branch("Px",&px);
outTree->Branch("Py",&py);
outTree->Branch("Pz",&pz);
outTree->Branch("Energy",&Energy);
outTree->Branch("EventNo",&EventNo);
// outTree->SetBranchAddress("NCherenkovs",&NCherenkovs);
// output txt file
ofstream outputFile( Form("output/%s_%s.txt",fileName.substr(0,fileName.find_last_of("." ) ).c_str(), rod.c_str() ) );
//// Get the tr.ansmission data
////
vector< double > transData;
string str;
ifstream file( Form("data/rod%s.txt", rod.c_str() ) );
if(!file.is_open()){
cout << Form("data/rod%s.txt didn't open", rod.c_str() ) << endl;
return;
}
int i = 0;
while(file){
getline(file,str);
transData.push_back( atof( str.c_str() ) );
i++;
}//End line loop
file.close();
// Load that data into a histogram so we can use the GetBin function
TH1F* hTrans = new TH1F( "TransData", "TransData", nBins, 197.2888, 1027.23596);
TH1F* hTrans2 = new TH1F( "TransData2", "TransData2", nBins, 197.2888, 1027.23596);
float trans, rand;
for(int bin = 0; bin < nBins; bin++){
hTrans->SetBinContent( bin, transData[bin] );
}
// Apply a modified baseline to the fused quartz spectrum
if(rod == "fusedQuartz"){
TF1* bLine = new TF1("baseline","0.0883*log(0.013087*(x-156.8) ) + 0.802",197,1028);
hTrans->Multiply(bLine);
}
//Get the number of events in the file
int nEvents = inTree->GetEntries();
//Find out how many events contain photons
float wl;
int realNevents = 0;
int totalCut = 0, totalPhotons = 0;
for (int ev = 0; ev < nEvents ; ev++){
inTree->GetEntry(ev);
if( Xf->size() > 0 ) realNevents++;
}
outputFile << Form("Total events = %d",realNevents) << endl;
for (int ev = 0; ev < nEvents ; ev++){
inTree->GetEntry(ev);
int nCut = 0;
int nPhotons = Xf->size();
if(nPhotons) outputFile << Form("Event %d",ev) << endl;
for (int k=0; k < nPhotons ; k++){ //begin loop over hits
Px = Pxf->at(k)*1e6;
Py = Pyf->at(k)*1e6;
Pz = Pzf->at(k)*1e6;
energy = sqrt( pow(Px,2) + pow(Py,2) + pow(Pz,2) );
// If the wavelength of the photon is < 197nm, set the transmission
// % equal to the value for the shortest wavelenght we have. Otherwise,
// get the data from the appropriate bin.
wl = 1240./energy;
if(wl < 197.2888){
trans = hTrans->GetBinContent(10);
}else{
trans = hTrans->GetBinContent( hTrans->FindBin( wl ) );
}
// Give the photon a random chance to make it weighted
// by transmission %
// if( rnd.Rndm() < trans ){
if( rnd.Rndm() < 1 ){
// If the photon is transmitted, add it to the vector
X = Xf->at(k) + xOffset;
Z = Zf->at(k) + zOffset;
x->push_back( X );
z->push_back( Z );
//We want momentum direction, not magnitude
Px/=energy;
Py/=energy;
Pz/=energy;
px->push_back( Px );
py->push_back( Py );
pz->push_back( Pz );
Energy->push_back( energy );
outputFile << Form("V,%17.11f,%17.11f,%17.11f",X,0.0,Z) << endl;
outputFile << Form("P,%17.14f,%17.14f,%17.14f,%17.13f",Px,Py,Pz,energy) << endl;
}else{
nCut++;
}//End cuts
}//end loop over fiber hits
if(nPhotons){
outputFile << Form("End event %d", ev) << endl;
cout << Form("%7d of %7d photons cut in event %3d",nCut,nPhotons,ev) << endl;
}
outTree->Fill();
totalCut += nCut;
totalPhotons += nPhotons;
//Clear vectors
x->clear();
z->clear();
px->clear();
py->clear();
pz->clear();
Energy->clear();
}//end loop over events
cout << Form("%5d of %5d photons cut %3.0f%%",totalCut,totalPhotons,100.*(float)totalCut/(float)totalPhotons ) << endl;
inFile->Close();
outFile->Close();
}
/*
* Apply quantum efficiency cuts to the output of the lightguide MC
*/
void PMTcuts(string fileName, int PMTmodel = 6091){
//Create input vectors
vector<double> *Xf=0, *Yf=0, *Zf=0, *Pxf=0, *Pyf=0, *Pzf=0, *Energyf=0, *NCherenkovs=0;
vector<int> *PIDf=0, *IDf=0, *EventNof=0;
// Set the input and output files up
//
TFile* inFile = new TFile(fileName.c_str());
if(!inFile->IsOpen()){
cout << "Input file didn't open" << endl;
return;
}
TTree* inTree = (TTree*)inFile->Get("lightGuide");
inTree->SetBranchAddress("hitX",&Xf);
inTree->SetBranchAddress("hitZ",&Zf);
inTree->SetBranchAddress("energy",&Energyf);
//// Get the quantum efficiency data
////
ifstream file( Form("data/model%dQE.txt",PMTmodel) );
if( !file.is_open() ){
cout << Form("data/model%dQE.csv didn't open... exiting",PMTmodel) << endl;
return;
}
vector<double> wavelength, efficiency;
string str;
size_t comma;
getline(file,str); //Burn the header
while( !file.eof() ){
getline(file,str);
if( str.length() ){
comma = str.find_first_of(",");
wavelength.push_back( atof( str.substr(0,comma).c_str() ) );
efficiency.push_back( atof( str.substr( comma+1, str.length()-comma ).c_str() ) );
}
}
TVectorD TvWL(wavelength.size(),&wavelength[0]);
TVectorD TvQE(efficiency.size(),&efficiency[0]);
TGraph* gQE = new TGraph(TvWL,TvQE);
//// begin loop over events
////
TRandom rnd;
int nCut, nPhotons, totalCut=0, totalPhotons=0;
double X,Z,wl,trans,rand;
int nEvents = inTree->GetEntries();
for (int ev = 0; ev < nEvents ; ev++){
inTree->GetEntry(ev);
nCut = 0;
nPhotons = Xf->size();
for (int k = 0; k < nPhotons; k++){ //begin loop over hits
// If the wavelength isn't defined by the QE data QE=0,
// so cut it and move on to the next photon
// Else give it a chance to make it under the curve
wl = 1240./(1e6*Energyf->at(k) );
if(wl < wavelength.front() || wl > wavelength.back()){
nCut++;
} else{
// Get transmission % from data
trans = gQE->Eval( wl );
// Give the photon a random chance to make it past the
// Quantum Efficiency check
rand = 100*rnd.Rndm();
if( 100*rnd.Rndm() > trans ){
nCut++;
}
}//end cuts
}//end nPhotons loop
totalCut += nCut;
totalPhotons += nPhotons;
cout << Form("%7d of %7d photons cut in event %3d",nCut,nPhotons,ev) << endl;
}//end event looop
cout << Form("%7d of %7d photons cut %3.0f%%",totalCut,totalPhotons,100.*(float)totalCut/(float)totalPhotons ) << endl;
inFile->Close();
delete inFile;
}//end PMTcuts
/* (experimental)
* Generate a waveform as the sum of single photon responses
* Pulse integral to amplitude conversion factor is a result of the following
* integral_21^∞ A*((t - t0)/10)^3.4 exp(-1/10 (t - t0)) dt = A*101.3581
* DRAW saves event display to .png
*/
void generateWaveform(string fileName, bool DRAW){
// Choose a path
bool theHardWay = true; // (superposition) Create a TF1 for every photon using each photon's time and energy/wavelength
bool theMiddleWay = false; // (per time bin) Create a TF1 for every time bin using the sum of that bin's photons energy/wavelength
bool theEasyWay = false; // (impulse) Create a single TF1 using the average photon energy (known value = 4.625eV) multiplied by nPhotons in that tile
float inputImpedance = 50.0; // ohms
float gain = 1e6; // unitless
float sampleFrequency = 5e9; // 5GHz
int nSamples = 1024;
float timeWindow = nSamples/sampleFrequency;
TH1D *h[16];
TCanvas *c = new TCanvas("simWF","Simulated Waveforms",800,600);
c->Divide(4,4);
TFile* inFile = new TFile(fileName.c_str());
if(inFile->IsZombie()){
cout << fileName << " didn't open" << endl;
return;
}
TTree* tree = (TTree*)inFile->Get("RPD1tree");
vector<double> *Px=0, *Py=0, *Pz=0, *time=0;
vector<int> *rodNo=0;
tree->SetBranchAddress("Px",&Px);
tree->SetBranchAddress("Py",&Py);
tree->SetBranchAddress("Pz",&Pz);
tree->SetBranchAddress("time",&time);
tree->SetBranchAddress("rodNo",&rodNo);
//Output root file
TFile* outFile = new TFile( "output.root", "RECREATE");
TTree* outTree = new TTree( "tree", "tree" );
vector< vector< float >* > waveforms(16,0);
for(int tile = 0; tile < 16; tile++){
waveforms[tile] = new vector< float >;
outTree->Branch( Form("Signal%d",tile), &waveforms[tile], Form("RPD tile %d simulated waveform",tile) );
if(DRAW) h[tile] = new TH1D( Form("tile%dwaveform",tile), Form("Tile %d Simulated Waveform;time (ns);Amplitude (mV)",tile), nSamples, 0, timeWindow);
}
if(theHardWay){
TGraph* g = getPMTdata( Form("data/model%dresponse.txt",2056) );
cout << "Average output = " << g->GetMean(2) << endl;
double wavelength, pulseIntegral, amplitude, photonEnergy;
vector< vector< TF1* > > pulses;
pulses.resize(16);
int nEntries = tree->GetEntries();
for(int ev = 0; ev < nEntries; ev++){
tree->GetEntry(ev);
int nHits = Px->size();
// if(nHits == 0) continue;
for(int hit = 0; hit < nHits; hit++){
if(hit%500 == 0) cout << "\r" << std::left << Form("Event %d, Hit %d",ev,hit) << flush;
photonEnergy = 1e6*sqrt( pow(Px->at(hit),2) + pow(Py->at(hit),2) + pow(Pz->at(hit),2) ); //eV
wavelength = 1240./photonEnergy;
pulseIntegral = g->Eval( wavelength )*gain*photonEnergy*1.6022e-19; // mC
amplitude = inputImpedance*pulseIntegral/((1./sampleFrequency)*101.358l); // mV
pulses[rodNo->at(hit)/16].push_back( new TF1( Form("ev%d",ev), PMTpulseFunction, -10, timeWindow + 10, 2) );
pulses[rodNo->at(hit)/16].back()->SetParameter( 0, amplitude);
pulses[rodNo->at(hit)/16].back()->SetParameter( 1, time->at(hit) + 21);
}//end hit loop
for(int tile = 0; tile < 16; tile++){
// Add dark current here to minimize calls to evaluate sum
addDarkCurrent(pulses[tile], 204.8, 1./sampleFrequency);
for(int bin = 0; bin < 1024; bin ++){
waveforms[tile]->push_back(evaluateSum( pulses[tile], 0.2*bin ) );
}//end waveform loop
addRandNoise(waveforms[tile], 204.8);
if(DRAW){
for(int bin = 0; bin < 1024; bin ++){
h[tile]->SetBinContent(bin, waveforms[tile]->at(bin));
}//end waveform loop
c->cd(tile+1);
h[tile]->Draw();
h[tile]->SetAxisRange(-10,750,"Y");
}
}//end tile loop
if(DRAW) c->Print( Form("event%d.png",ev) );
outTree->Fill();
//Delete all TF1
for(int tile = 0; tile < 16; tile++){
waveforms[tile]->clear();
int nHits = pulses[tile].size();
if(nHits == 0) continue;
for(int hit = 0; hit < nHits; hit++){
delete pulses[tile][hit];
}//end nHits loop
pulses[tile].clear();
}//end tile loop
}//end event loop
cout << endl;
}else if( theMiddleWay ){
}else if( theEasyWay ){
}//end the easy way
outFile->Close();
inFile->Close();
delete inFile;
delete outFile; //Deletes histograms
for(int tile = 0; tile < 16; tile++){
delete waveforms[tile];
}
if(!DRAW) delete c;
}
void addRandNoise(vector< float >* vec, float timeInterval){
float sigma = 5.0; //mV
float mean = 0.0; //mV
int size = vec->size();
TRandom2 rnd;
rnd.SetSeed(gSystem->Now());
for(int i = 0; i < size; i++){
vec->at(i) += rnd.Gaus(mean,sigma);
}
}
/*
* Adds dark current pulses to the given waveform vector
* Time interval is the duration of the waveform and must be given in ns
*/
void addDarkCurrent(vector< TF1* >& pulses, float timeInterval, float timeBinWidth){
// For R2496
// Average Cathode Radiant Sensitivity = 21.7475 mA/W (mC/J), gives pulse amplitude of 0.214561 for an average pulse
// Dark Current = 2nA typical, 50 nA max
float avePhotonEnergy = 4.625*1.6022e-19; // Joule
float aveRadSens = 21.7475; // mA/W aka mC/J
float gain = 1e6; // unitless
float aveChargePerPulse = avePhotonEnergy*aveRadSens*gain; // mC
float darkCurrent = 2e-6; // mA aka mC/s
float aveNpulsePerNanoSec = 1e-9*darkCurrent/aveChargePerPulse; // 1/ns
float inputImpedance = 50.0; // ohms
// The units we need to generate the dark current
float buffer = 80; // ns
timeInterval += buffer;
float aveNpulses = timeInterval*aveNpulsePerNanoSec; // unitless. Probably want to randomize this a bit
float avePulseAmp = aveChargePerPulse*inputImpedance/(timeBinWidth*101.3581); // mV
TRandom2 rnd;
rnd.SetSeed(gSystem->Now());
int nPulses = rnd.Poisson( aveNpulses );
for(int i = 0; i < nPulses; i++){
float delay = rnd.Uniform(-1*buffer,timeInterval);
pulses.push_back(new TF1( Form("darkCurrent%d",i), PMTpulseFunction, delay, timeInterval, 2) );
pulses.back()->SetParameter( 0, rnd.Gaus( avePulseAmp, 2*avePulseAmp) );
pulses.back()->SetParameter( 1, delay );
}
}
TGraph* getPMTdata(string fileName){
vector<double> x, y;
ifstream file( fileName.c_str() );
if( !file.is_open() ){
cout << Form( "%s didn't open... exiting", fileName.c_str() ) << endl;
return NULL;
}
string str;
size_t comma;
getline(file,str); //Burn the header
while( !file.eof() ){
getline(file,str);
if( str.length() ){
comma = str.find_first_of(",");
x.push_back( atof( str.substr(0,comma).c_str() ) );
y.push_back( atof( str.substr( comma+1, str.length()-comma ).c_str() ) );
}
}
TVectorD TvX(x.size(),&x[0]);
TVectorD TvY(y.size(),&y[0]);
TGraph* g = new TGraph(TvX,TvY);
return g;
}
double evaluateSum( vector< TF1* >& funcVec, double x ){
double value = 0.0;
for(int i = 0; i < funcVec.size(); i++){
value += funcVec.at(i)->Eval(x);
}
return value;
}
#include "TFile.h"
#include "TH1F.h"
#include "TGraph.h"
#include "TF1.h"
#include "TRandom.h"
#include <TTree.h>
#include <vector>
#include <fstream>
using namespace std;
TGraph* getPMTdata( string fileName );
double evaluateSum( vector< TF1* >& funcVec, double x );
void generateWaveform(string fileName, bool DRAW = false);
void addDarkCurrent(vector< TF1* >& vec, float timeInterval, float timeBinWidth);
void addRandNoise(vector< float >* vec, float timeInterval);
void addTimingJitter(vector< float >* vec, float timeInterval);
Double_t PMTpulseFunction(Double_t *t, Double_t *par){
Float_t tt = t[0];
if(tt < par[1]) return 0.0; // cut out the time where the function behaves poorly
Double_t f = par[0]*pow((tt-par[1])/10.0, 3.4)*exp(-(tt-par[1])/10.0);
return f;
}
void batchWaveformFolder(const char *dirname="./Output/", const char *ext=".root")
{ TSystemDirectory dir(dirname, dirname);
TList *files = dir.GetListOfFiles();
if (files) { TSystemFile *file;
TString fname;
TIter next(files);
while ((file=(TSystemFile*)next()))
{ fname = file->GetName();
if (!file->IsDirectory() && fname.EndsWith(ext))
{ generateWaveform(fname.Data()); } } } }
void batchWaveform(){
int fileStart = 0;
int fileStop = 100;
for(int fileNum = fileStart; fileNum < fileStop; fileNum++){
generateWaveform( Form("filename%d.root",fileNum) );
}
}
/*
* Calculate the number of Cherenkov photons likely to make it out
* of the ZDC starting from the output of the JZCaPA MC. Accounting
* for losses due to absorption in fused silica after heavy irradiation.
*/
void cutZDCoutput(string fileName, string rod, double xOffset, double zOffset){
int nBins = 3649;
double energy,X,Z,Px,Py,Pz;
TRandom rnd;
//Create input vectors
vector<double> *Xf=0, *Yf=0, *Zf=0, *Pxf=0, *Pyf=0, *Pzf=0, *Energyf=0, *NCherenkovs=0;
vector<int> *PIDf=0, *IDf=0, *EventNof=0;
//Create output vectors
vector<double> *x=0, *z=0,*x_0=0, *z_0=0, *px=0, *py=0, *pz=0, *Energy=0;
vector<int> *EventNo=0;
// Set the input and output files up
//
TFile* inFile = new TFile(fileName.c_str());
TTree* inTree = (TTree*)inFile->Get("ZDCtree");
inTree->SetBranchAddress("X",&Xf);
inTree->SetBranchAddress("Z",&Zf);
inTree->SetBranchAddress("Px",&Pxf);
inTree->SetBranchAddress("Py",&Pyf);
inTree->SetBranchAddress("Pz",&Pzf);
//inTree->SetBranchAddress("EventNo",&EventNf);
// intree->SetBranchAddress("NCherenkovs",&NCherenkovs);
//Output root file
TFile* outFile = new TFile( Form("output/%s_%s.root",fileName.substr(0,fileName.find_last_of(".")).c_str(), rod.c_str() ), "RECREATE");
TTree* outTree = new TTree( "tree", "tree" );
// outTree->Branch("X0",&x_0);
// outTree->Branch("Z0",&z_0);
outTree->Branch("X",&x);
outTree->Branch("Z",&z);
outTree->Branch("Px",&px);
outTree->Branch("Py",&py);
outTree->Branch("Pz",&pz);
outTree->Branch("Energy",&Energy);
outTree->Branch("EventNo",&EventNo);
// outTree->SetBranchAddress("NCherenkovs",&NCherenkovs);
// output txt file
ofstream outputFile( Form("output/%s_%s.txt",fileName.substr(0,fileName.find_last_of("." ) ).c_str(), rod.c_str() ) );
//// Get the tr.ansmission data
////
vector< double > transData;
string str;
ifstream file( Form("data/rod%s.txt", rod.c_str() ) );
if(!file.is_open()){
cout << Form("data/rod%s.txt didn't open", rod.c_str() ) << endl;
return;
}
int i = 0;
while(file){
getline(file,str);
transData.push_back( atof( str.c_str() ) );
i++;
}//End line loop
file.close();
// Load that data into a histogram so we can use the GetBin function
TH1F* hTrans = new TH1F( "TransData", "TransData", nBins, 197.2888, 1027.23596);
TH1F* hTrans2 = new TH1F( "TransData2", "TransData2", nBins, 197.2888, 1027.23596);
float trans, rand;
for(int bin = 0; bin < nBins; bin++){
hTrans->SetBinContent( bin, transData[bin] );
}
// Apply a modified baseline to the fused quartz spectrum
if(rod == "fusedQuartz"){
TF1* bLine = new TF1("baseline","0.0883*log(0.013087*(x-156.8) ) + 0.802",197,1028);
hTrans->Multiply(bLine);
}
//Get the number of events in the file
int nEvents = inTree->GetEntries();
//Find out how many events contain photons
float wl;
int realNevents = 0;
int totalCut = 0, totalPhotons = 0;
for (int ev = 0; ev < nEvents ; ev++){
inTree->GetEntry(ev);
if( Xf->size() > 0 ) realNevents++;
}
outputFile << Form("Total events = %d",realNevents) << endl;
for (int ev = 0; ev < nEvents ; ev++){
inTree->GetEntry(ev);
int nCut = 0;
int nPhotons = Xf->size();
if(nPhotons) outputFile << Form("Event %d",ev) << endl;
for (int k=0; k < nPhotons ; k++){ //begin loop over hits
Px = Pxf->at(k)*1e6;
Py = Pyf->at(k)*1e6;
Pz = Pzf->at(k)*1e6;
energy = sqrt( pow(Px,2) + pow(Py,2) + pow(Pz,2) );
// If the wavelength of the photon is < 197nm, set the transmission
// % equal to the value for the shortest wavelenght we have. Otherwise,
// get the data from the appropriate bin.
wl = 1240./energy;
if(wl < 197.2888){
trans = hTrans->GetBinContent(10);
}else{
trans = hTrans->GetBinContent( hTrans->FindBin( wl ) );
}
// Give the photon a random chance to make it weighted
// by transmission %
// if( rnd.Rndm() < trans ){
if( rnd.Rndm() < 1 ){
// If the photon is transmitted, add it to the vector
X = Xf->at(k) + xOffset;
Z = Zf->at(k) + zOffset;
x->push_back( X );
z->push_back( Z );
//We want momentum direction, not magnitude
Px/=energy;
Py/=energy;
Pz/=energy;
px->push_back( Px );
py->push_back( Py );
pz->push_back( Pz );
Energy->push_back( energy );
outputFile << Form("V,%17.11f,%17.11f,%17.11f",X,0.0,Z) << endl;
outputFile << Form("P,%17.14f,%17.14f,%17.14f,%17.13f",Px,Py,Pz,energy) << endl;
}else{
nCut++;
}//End cuts
}//end loop over fiber hits
if(nPhotons){
outputFile << Form("End event %d", ev) << endl;
cout << Form("%7d of %7d photons cut in event %3d",nCut,nPhotons,ev) << endl;
}
outTree->Fill();
totalCut += nCut;
totalPhotons += nPhotons;
//Clear vectors
x->clear();
z->clear();
px->clear();
py->clear();
pz->clear();
Energy->clear();
}//end loop over events
cout << Form("%5d of %5d photons cut %3.0f%%",totalCut,totalPhotons,100.*(float)totalCut/(float)totalPhotons ) << endl;
inFile->Close();
outFile->Close();
}
/*
* Apply quantum efficiency cuts to the output of the lightguide MC
*/
void PMTcuts(string fileName, int PMTmodel = 6091){
//Create input vectors
vector<double> *Xf=0, *Yf=0, *Zf=0, *Pxf=0, *Pyf=0, *Pzf=0, *Energyf=0, *NCherenkovs=0;
vector<int> *PIDf=0, *IDf=0, *EventNof=0;
// Set the input and output files up
//
TFile* inFile = new TFile(fileName.c_str());
if(!inFile->IsOpen()){
cout << "Input file didn't open" << endl;
return;
}
TTree* inTree = (TTree*)inFile->Get("lightGuide");
inTree->SetBranchAddress("hitX",&Xf);
inTree->SetBranchAddress("hitZ",&Zf);
inTree->SetBranchAddress("energy",&Energyf);
//// Get the quantum efficiency data
////
ifstream file( Form("data/model%dQE.txt",PMTmodel) );
if( !file.is_open() ){
cout << Form("data/model%dQE.csv didn't open... exiting",PMTmodel) << endl;
return;
}
vector<double> wavelength, efficiency;
string str;
size_t comma;
getline(file,str); //Burn the header
while( !file.eof() ){
getline(file,str);
if( str.length() ){
comma = str.find_first_of(",");
wavelength.push_back( atof( str.substr(0,comma).c_str() ) );
efficiency.push_back( atof( str.substr( comma+1, str.length()-comma ).c_str() ) );
}
}
TVectorD TvWL(wavelength.size(),&wavelength[0]);
TVectorD TvQE(efficiency.size(),&efficiency[0]);
TGraph* gQE = new TGraph(TvWL,TvQE);
//// begin loop over events
////
TRandom rnd;
int nCut, nPhotons, totalCut=0, totalPhotons=0;
double X,Z,wl,trans,rand;
int nEvents = inTree->GetEntries();
for (int ev = 0; ev < nEvents ; ev++){
inTree->GetEntry(ev);
nCut = 0;
nPhotons = Xf->size();
for (int k = 0; k < nPhotons; k++){ //begin loop over hits
// If the wavelength isn't defined by the QE data QE=0,
// so cut it and move on to the next photon
// Else give it a chance to make it under the curve
wl = 1240./(1e6*Energyf->at(k) );
if(wl < wavelength.front() || wl > wavelength.back()){
nCut++;
} else{
// Get transmission % from data
trans = gQE->Eval( wl );
// Give the photon a random chance to make it past the
// Quantum Efficiency check
rand = 100*rnd.Rndm();
if( 100*rnd.Rndm() > trans ){
nCut++;
}
}//end cuts
}//end nPhotons loop
totalCut += nCut;
totalPhotons += nPhotons;
cout << Form("%7d of %7d photons cut in event %3d",nCut,nPhotons,ev) << endl;
}//end event looop
cout << Form("%7d of %7d photons cut %3.0f%%",totalCut,totalPhotons,100.*(float)totalCut/(float)totalPhotons ) << endl;
inFile->Close();
delete inFile;
}//end PMTcuts
/* (experimental)
* Generate a waveform as the sum of single photon responses
* Pulse integral to amplitude conversion factor is a result of the following
* integral_21^∞ A*((t - t0)/10)^3.4 exp(-1/10 (t - t0)) dt = A*101.3581
* DRAW saves event display to .png
*/
void generateWaveform(string fileName, bool DRAW){
// Choose a path
bool theHardWay = true; // (superposition) Create a TF1 for every photon using each photon's time and energy/wavelength
bool theMiddleWay = false; // (per time bin) Create a TF1 for every time bin using the sum of that bin's photons energy/wavelength
bool theEasyWay = false; // (impulse) Create a single TF1 using the average photon energy (known value = 4.625eV) multiplied by nPhotons in that tile
float inputImpedance = 50.0; // ohms
float gain = 1e6; // unitless
float sampleFrequency = 5e9; // 5GHz
int nSamples = 1024;
float timeWindow = nSamples/sampleFrequency;
TH1D *h[16];
TCanvas *c = new TCanvas("simWF","Simulated Waveforms",800,600);
c->Divide(4,4);
TFile* inFile = new TFile(fileName.c_str());
if(inFile->IsZombie()){
cout << fileName << " didn't open" << endl;
return;
}
TTree* tree = (TTree*)inFile->Get("RPD1tree");
vector<double> *Px=0, *Py=0, *Pz=0, *time=0;
vector<int> *rodNo=0;
tree->SetBranchAddress("Px",&Px);
tree->SetBranchAddress("Py",&Py);
tree->SetBranchAddress("Pz",&Pz);
tree->SetBranchAddress("time",&time);
tree->SetBranchAddress("rodNo",&rodNo);
//Output root file
TFile* outFile = new TFile( "output.root", "RECREATE");
TTree* outTree = new TTree( "tree", "tree" );
vector< vector< float >* > waveforms(16,0);
for(int tile = 0; tile < 16; tile++){
waveforms[tile] = new vector< float >;
outTree->Branch( Form("Signal%d",tile), &waveforms[tile], Form("RPD tile %d simulated waveform",tile) );
if(DRAW) h[tile] = new TH1D( Form("tile%dwaveform",tile), Form("Tile %d Simulated Waveform;time (ns);Amplitude (mV)",tile), nSamples, 0, timeWindow);
}
if(theHardWay){
TGraph* g = getPMTdata( Form("data/model%dresponse.txt",2056) );
cout << "Average output = " << g->GetMean(2) << endl;
double wavelength, pulseIntegral, amplitude, photonEnergy;
vector< vector< TF1* > > pulses;
pulses.resize(16);
int nEntries = tree->GetEntries();
for(int ev = 0; ev < nEntries; ev++){
tree->GetEntry(ev);
int nHits = Px->size();
// if(nHits == 0) continue;
for(int hit = 0; hit < nHits; hit++){
if(hit%500 == 0) cout << "\r" << std::left << Form("Event %d, Hit %d",ev,hit) << flush;
photonEnergy = 1e6*sqrt( pow(Px->at(hit),2) + pow(Py->at(hit),2) + pow(Pz->at(hit),2) ); //eV
wavelength = 1240./photonEnergy;
pulseIntegral = g->Eval( wavelength )*gain*photonEnergy*1.6022e-19; // mC
amplitude = inputImpedance*pulseIntegral/((1./sampleFrequency)*101.358l); // mV
pulses[rodNo->at(hit)/16].push_back( new TF1( Form("ev%d",ev), PMTpulseFunction, -10, timeWindow + 10, 2) );
pulses[rodNo->at(hit)/16].back()->SetParameter( 0, amplitude);
pulses[rodNo->at(hit)/16].back()->SetParameter( 1, time->at(hit) + 21);
}//end hit loop
for(int tile = 0; tile < 16; tile++){
// Add dark current here to minimize calls to evaluate sum
addDarkCurrent(pulses[tile], 204.8, 1./sampleFrequency);
for(int bin = 0; bin < 1024; bin ++){
waveforms[tile]->push_back(evaluateSum( pulses[tile], 0.2*bin ) );
}//end waveform loop
addRandNoise(waveforms[tile], 204.8);
if(DRAW){
for(int bin = 0; bin < 1024; bin ++){
h[tile]->SetBinContent(bin, waveforms[tile]->at(bin));
}//end waveform loop
c->cd(tile+1);
h[tile]->Draw();
h[tile]->SetAxisRange(-10,750,"Y");
}
}//end tile loop
if(DRAW) c->Print( Form("event%d.png",ev) );
outTree->Fill();
//Delete all TF1
for(int tile = 0; tile < 16; tile++){
waveforms[tile]->clear();
int nHits = pulses[tile].size();
if(nHits == 0) continue;
for(int hit = 0; hit < nHits; hit++){
delete pulses[tile][hit];
}//end nHits loop
pulses[tile].clear();
}//end tile loop
}//end event loop
cout << endl;
}else if( theMiddleWay ){
}else if( theEasyWay ){
}//end the easy way
outFile->Close();
inFile->Close();
delete inFile;
delete outFile; //Deletes histograms
for(int tile = 0; tile < 16; tile++){
delete waveforms[tile];
}
if(!DRAW) delete c;
}
void addRandNoise(vector< float >* vec, float timeInterval){
float sigma = 5.0; //mV
float mean = 0.0; //mV
int size = vec->size();
TRandom2 rnd;
rnd.SetSeed(gSystem->Now());
for(int i = 0; i < size; i++){
vec->at(i) += rnd.Gaus(mean,sigma);
}
}
/*
* Adds dark current pulses to the given waveform vector
* Time interval is the duration of the waveform and must be given in ns
*/
void addDarkCurrent(vector< TF1* >& pulses, float timeInterval, float timeBinWidth){
// For R2496
// Average Cathode Radiant Sensitivity = 21.7475 mA/W (mC/J), gives pulse amplitude of 0.214561 for an average pulse
// Dark Current = 2nA typical, 50 nA max
float avePhotonEnergy = 4.625*1.6022e-19; // Joule
float aveRadSens = 21.7475; // mA/W aka mC/J
float gain = 1e6; // unitless
float aveChargePerPulse = avePhotonEnergy*aveRadSens*gain; // mC
float darkCurrent = 2e-6; // mA aka mC/s
float aveNpulsePerNanoSec = 1e-9*darkCurrent/aveChargePerPulse; // 1/ns
float inputImpedance = 50.0; // ohms
// The units we need to generate the dark current
float buffer = 80; // ns
timeInterval += buffer;
float aveNpulses = timeInterval*aveNpulsePerNanoSec; // unitless. Probably want to randomize this a bit
float avePulseAmp = aveChargePerPulse*inputImpedance/(timeBinWidth*101.3581); // mV
TRandom2 rnd;
rnd.SetSeed(gSystem->Now());
int nPulses = rnd.Poisson( aveNpulses );
for(int i = 0; i < nPulses; i++){
float delay = rnd.Uniform(-1*buffer,timeInterval);
pulses.push_back(new TF1( Form("darkCurrent%d",i), PMTpulseFunction, delay, timeInterval, 2) );
pulses.back()->SetParameter( 0, rnd.Gaus( avePulseAmp, 2*avePulseAmp) );
pulses.back()->SetParameter( 1, delay );
}
}
TGraph* getPMTdata(string fileName){
vector<double> x, y;
ifstream file( fileName.c_str() );
if( !file.is_open() ){
cout << Form( "%s didn't open... exiting", fileName.c_str() ) << endl;
return NULL;
}
string str;
size_t comma;
getline(file,str); //Burn the header
while( !file.eof() ){
getline(file,str);
if( str.length() ){
comma = str.find_first_of(",");
x.push_back( atof( str.substr(0,comma).c_str() ) );
y.push_back( atof( str.substr( comma+1, str.length()-comma ).c_str() ) );
}
}
TVectorD TvX(x.size(),&x[0]);
TVectorD TvY(y.size(),&y[0]);
TGraph* g = new TGraph(TvX,TvY);
return g;
}
double evaluateSum( vector< TF1* >& funcVec, double x ){
double value = 0.0;
for(int i = 0; i < funcVec.size(); i++){
value += funcVec.at(i)->Eval(x);
}
return value;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment