Skip to content
Snippets Groups Projects
Commit 7b8b5e12 authored by shengy3's avatar shengy3
Browse files

Remove applying DRS4 response

parent f1dd0498
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 590 additions and 589 deletions
outputDigitization.cc
+ 590
− 589
View file @ 7b8b5e12
#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* >& vec, double x );
void addTimingJitter ( vector< float >* vec, float timeInterval);
void generateWaveform ( string fileName, bool DRAW = false);
void getFile (string filename, /*out*/ ifstream& file);
/*
* This PMT waveform function came from the ZDC_PileUpTool in Athena
* The function was made piecewise to avoid extreme values before the
* delay time
* Parameter 0: Amplitude
* Parameter 1: Delay
* Time units are in ns
*/
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="./data/", 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(Form("%s%s",dirname, 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 transmission 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 wavelength 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) Not implemented yet. Create a TF1 for every time bin using the sum of that bin's photons energy/wavelength
bool theEasyWay = false; // (impulse) Not implemented yet. Create a single TF1 using the average photon energy (known value = 4.625eV) multiplied by nPhotons in that tile
////////////////////// Readout constants //////////////////////
float sampleFrequency = 5.; // GHz
int nSamples = 1024;
float timeBinWidth = 1./sampleFrequency; // ns
float timeWindow = nSamples*timeBinWidth; // ns
//Get the intgral of a pulse with unit amplitude
float triggerDelay = 21.0;
TF1 *testPulse = new TF1("testPulse", PMTpulseFunction, 0.0, timeWindow, 2);
testPulse->SetParameters(1.0,triggerDelay);
float unitPulseIntegral = testPulse->Integral( 0.0, timeWindow, 1e-4 );
TF1* DRS4_response = new TF1("DRS4Response", "(x<700)*(0.27698+0.965037*x) + \(x>700 && x<780)*( -9.63227511e+6 + 5.24719827e4*x -107.112621*pow(x,2) + 9.71101289e-2*pow(x,3) - (3.29905732e-5)*pow(x,4) ) + \(x>780)*( -6.29891342e+4 + 2.31093738e+2*x -2.78617584e-1*pow(x,2) + 1.11952166e-4*pow(x,3) )", 0, 1000);
////////////////////// PMT constants //////////////////////
// Will probably make a PMT object when appropriate
// For R2496
// Dark Current = 2nA typical, 50 nA max
TGraph* g = getPMTdata( Form("model%dresponse.txt",2056) );
if(g == NULL) return;
float avePhotonEnergy = (1240./g->GetMean(1))*1.6022e-19; // Joule
float aveRadSens = g->GetMean(2); // mA/W aka mC/J
float gain = 1e6; // unitless
float darkCurrent = 2e-6; // mA aka mC/s. Mfr claims this is 2nA typical, but up to 50nA
float outputImpedance = 50.0; // ohms
float aveChargePerPulse = avePhotonEnergy*aveRadSens*gain; // mC
float aveNpulsePerNanoSec = 1e-9*darkCurrent/aveChargePerPulse; // 1/ns
float aveNpulses = timeWindow*aveNpulsePerNanoSec; // unitless
float avePulseAmp = aveChargePerPulse*outputImpedance/(1e-9*timeBinWidth*unitPulseIntegral); // mV
////////////////////// Noise constants //////////////////////
float sigma = 5.0; //mV
float mean = 0.0; //mV
TH1D *h[16];
TCanvas *c = new TCanvas("simWF","Simulated Waveforms",800,600);
c->Divide(4,4);
// Input root file
TFile* inFile = new TFile(fileName.c_str());
if(inFile->IsZombie()){
cout << fileName << " didn't open" << endl;
return;
}
TTree* eventDatatree = (TTree*)inFile->Get("EventData");
TTree* tree = (TTree*)inFile->Get("RPD1tree");
vector<double> *Px=0, *Py=0, *Pz=0, *time=0, *LastStepZf = 0;
vector<int> *rodNo=0;
double Xf=0, Yf=0, Zf=0;
tree->SetBranchAddress("Px",&Px);
tree->SetBranchAddress("Py",&Py);
tree->SetBranchAddress("Pz",&Pz);
tree->SetBranchAddress("time",&time);
tree->SetBranchAddress("rodNo",&rodNo);
eventDatatree->SetBranchAddress("gunPosX",&Xf);
eventDatatree->SetBranchAddress("gunPosY",&Yf);
eventDatatree->SetBranchAddress("gunPosZ",&Zf);
eventDatatree->SetBranchAddress("lastStepZ",&LastStepZf);
// Output root file
string outputName = fileName;
if(outputName.find("/") != string::npos) outputName.erase( 0, outputName.find_last_of("/") + 1 );
if(outputName.find(".") != string::npos) outputName.erase( outputName.find_last_of(".") );
TFile* outFile = new TFile( Form("%sWF.root", outputName.c_str() ), "RECREATE");
TTree* outTree = new TTree( "tree", "tree" );
outTree->Branch("gunPosX",&Xf);
outTree->Branch("gunPosY",&Yf);
outTree->Branch("gunPosZ",&Zf);
outTree->Branch("lastStepZ",&LastStepZf);
vector< vector< float >* > waveforms(16,0);
for(int tile = 0; tile < 16; tile++){
waveforms[tile] = new vector< float >;
outTree->Branch( Form("RawSignal%d",tile), &waveforms[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){
double wavelength, pulseIntegral, amplitude, photonEnergy;
vector< vector< TF1* > > pulses;
pulses.resize(16);
TRandom2 rnd;
rnd.SetSeed(gSystem->Now());
int nEntries = tree->GetEntries();
for(int ev = 0; ev < nEntries; ev++){
eventDatatree->GetEntry(ev);
tree->GetEntry(ev);
int nHits = Px->size();
for(int hit = 0; hit < nHits; hit++){
if(hit%500 == 0) cout << "\r" << std::left << Form("Event %d, Hit %d",ev,hit) << flush;
// Get the photon energy and determine the pulse amplitude based on that energy
photonEnergy = 1e6*sqrt( pow(Px->at(hit),2) + pow(Py->at(hit),2) + pow(Pz->at(hit),2) ); //eV
wavelength = 1240./photonEnergy; // nm
pulseIntegral = g->Eval( wavelength )*gain*photonEnergy*1.6022e-19; // mC
amplitude = outputImpedance*pulseIntegral/(1e-9*timeBinWidth*unitPulseIntegral); // mV
// Add the pulse to the pulses vector for each tile
pulses[rodNo->at(hit)/16].push_back( new TF1( Form("ev%d",ev), PMTpulseFunction, 0., timeWindow, 2) );
pulses[rodNo->at(hit)/16].back()->SetParameters( amplitude, time->at(hit) + triggerDelay);
}//end hit loop
for(int tile = 0; tile < 16; tile++){
// Add dark current pulses
int nPulses = rnd.Poisson( aveNpulses );
for(int i = 0; i < nPulses; i++){
pulses[tile].push_back(new TF1( Form("darkCurrent%d",i), PMTpulseFunction, 0., timeWindow, 2) );
pulses[tile].back()->SetParameters( rnd.Gaus( avePulseAmp, 2*avePulseAmp), rnd.Uniform(-80,timeWindow) );
}
// Evaluate the TF1s for each time bin + gaussian noise
for(int bin = 0; bin < 1024; bin ++){
waveforms[tile]->push_back(DRS4_response->Eval(evaluateSum( pulses[tile], timeBinWidth*bin ) + rnd.Gaus(mean,sigma) ));
}//end waveform loop
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 objects, clear vectors, reset histograms
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();
if(DRAW)h[tile]->Reset();
}//end tile loop
}//end event loop
cout << endl;
}else if( theMiddleWay ){
// Not implemented yet
}else if( theEasyWay ){
// Not implemented yet
}//end the easy way
outFile->Write();
outFile->Close();
delete outFile; //Deletes histograms
inFile->Close();
delete inFile;
for(int tile = 0; tile < 16; tile++){
delete waveforms[tile];
}
if(!DRAW) delete c;
}
TGraph* getPMTdata(string fileName){
vector<double> x, y;
ifstream file( Form("data/%s", fileName.c_str() ) );
if( !file.is_open() ){
file.open( 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;
}
/*
* Generate timing data to be used by JZCaPA
* Argument is digitization speed in GHz
*/
void generateTimingData( float digiSpeed = 5.0 ){
ofstream outputFile( Form("output/MC_%3.1fGHz.txt", digiSpeed ) );
outputFile << Form("%3.1fGHz simulated timing data", digiSpeed) << endl;
outputFile << " Module1 Module2 Module3 Module4 Module5" << endl;
float binWidth = 1./digiSpeed;
for(int i = 0; i < 1024; i++){
outputFile << Form("%8.3f, %8.3f, %8.3f, %8.3f, %8.3f", i*binWidth, i*binWidth, i*binWidth, i*binWidth, i*binWidth) << endl;
}
}
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;
}
bool checkFileExistencein(const string& filename){
ifstream f(filename.c_str());
return f.is_open();
}
void getFile(string filename, /*out*/ ifstream& file){
const bool file_exists = checkFileExistencein(Form("data/%s", filename.c_str() ) );
if (!file_exists) {
cout << "File in " << Form("data/%s", filename.c_str() ) << " not found." << endl;
cout << "Search " << filename << endl;
bool file_exists = checkFileExistencein(filename);
if (file_exists) {
file.open(filename.c_str());
cout << filename << " is found in the same level" << endl;
}
}
}
#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* >& vec, double x );
void addTimingJitter ( vector< float >* vec, float timeInterval);
void generateWaveform ( string fileName, bool DRAW = false);
void getFile (string filename, /*out*/ ifstream& file);
/*
* This PMT waveform function came from the ZDC_PileUpTool in Athena
* The function was made piecewise to avoid extreme values before the
* delay time
* Parameter 0: Amplitude
* Parameter 1: Delay
* Time units are in ns
*/
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="./data/", 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(Form("%s%s",dirname, 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 transmission 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 wavelength 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) Not implemented yet. Create a TF1 for every time bin using the sum of that bin's photons energy/wavelength
bool theEasyWay = false; // (impulse) Not implemented yet. Create a single TF1 using the average photon energy (known value = 4.625eV) multiplied by nPhotons in that tile
////////////////////// Readout constants //////////////////////
float sampleFrequency = 5.; // GHz
int nSamples = 1024;
float timeBinWidth = 1./sampleFrequency; // ns
float timeWindow = nSamples*timeBinWidth; // ns
//Get the intgral of a pulse with unit amplitude
float triggerDelay = 21.0;
TF1 *testPulse = new TF1("testPulse", PMTpulseFunction, 0.0, timeWindow, 2);
testPulse->SetParameters(1.0,triggerDelay);
float unitPulseIntegral = testPulse->Integral( 0.0, timeWindow, 1e-4 );
TF1* DRS4_response = new TF1("DRS4Response", "(x<700)*(0.27698+0.965037*x) + \(x>700 && x<780)*( -9.63227511e+6 + 5.24719827e4*x -107.112621*pow(x,2) + 9.71101289e-2*pow(x,3) - (3.29905732e-5)*pow(x,4) ) + \(x>780)*( -6.29891342e+4 + 2.31093738e+2*x -2.78617584e-1*pow(x,2) + 1.11952166e-4*pow(x,3) )", 0, 1000);
////////////////////// PMT constants //////////////////////
// Will probably make a PMT object when appropriate
// For R2496
// Dark Current = 2nA typical, 50 nA max
TGraph* g = getPMTdata( Form("model%dresponse.txt",2056) );
if(g == NULL) return;
float avePhotonEnergy = (1240./g->GetMean(1))*1.6022e-19; // Joule
float aveRadSens = g->GetMean(2); // mA/W aka mC/J
float gain = 1e6; // unitless
float darkCurrent = 2e-6; // mA aka mC/s. Mfr claims this is 2nA typical, but up to 50nA
float outputImpedance = 50.0; // ohms
float aveChargePerPulse = avePhotonEnergy*aveRadSens*gain; // mC
float aveNpulsePerNanoSec = 1e-9*darkCurrent/aveChargePerPulse; // 1/ns
float aveNpulses = timeWindow*aveNpulsePerNanoSec; // unitless
float avePulseAmp = aveChargePerPulse*outputImpedance/(1e-9*timeBinWidth*unitPulseIntegral); // mV
////////////////////// Noise constants //////////////////////
float sigma = 5.0; //mV
float mean = 0.0; //mV
TH1D *h[16];
TCanvas *c = new TCanvas("simWF","Simulated Waveforms",800,600);
c->Divide(4,4);
// Input root file
TFile* inFile = new TFile(fileName.c_str());
if(inFile->IsZombie()){
cout << fileName << " didn't open" << endl;
return;
}
TTree* eventDatatree = (TTree*)inFile->Get("EventData");
TTree* tree = (TTree*)inFile->Get("RPD1tree");
vector<double> *Px=0, *Py=0, *Pz=0, *time=0, *LastStepZf = 0;
vector<int> *rodNo=0;
double Xf=0, Yf=0, Zf=0;
tree->SetBranchAddress("Px",&Px);
tree->SetBranchAddress("Py",&Py);
tree->SetBranchAddress("Pz",&Pz);
tree->SetBranchAddress("time",&time);
tree->SetBranchAddress("rodNo",&rodNo);
eventDatatree->SetBranchAddress("gunPosX",&Xf);
eventDatatree->SetBranchAddress("gunPosY",&Yf);
eventDatatree->SetBranchAddress("gunPosZ",&Zf);
eventDatatree->SetBranchAddress("lastStepZ",&LastStepZf);
// Output root file
string outputName = fileName;
if(outputName.find("/") != string::npos) outputName.erase( 0, outputName.find_last_of("/") + 1 );
if(outputName.find(".") != string::npos) outputName.erase( outputName.find_last_of(".") );
TFile* outFile = new TFile( Form("%sWF.root", outputName.c_str() ), "RECREATE");
TTree* outTree = new TTree( "tree", "tree" );
outTree->Branch("gunPosX",&Xf);
outTree->Branch("gunPosY",&Yf);
outTree->Branch("gunPosZ",&Zf);
outTree->Branch("lastStepZ",&LastStepZf);
vector< vector< float >* > waveforms(16,0);
for(int tile = 0; tile < 16; tile++){
waveforms[tile] = new vector< float >;
outTree->Branch( Form("RawSignal%d",tile), &waveforms[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){
double wavelength, pulseIntegral, amplitude, photonEnergy;
vector< vector< TF1* > > pulses;
pulses.resize(16);
TRandom2 rnd;
rnd.SetSeed(gSystem->Now());
int nEntries = tree->GetEntries();
for(int ev = 0; ev < nEntries; ev++){
eventDatatree->GetEntry(ev);
tree->GetEntry(ev);
int nHits = Px->size();
for(int hit = 0; hit < nHits; hit++){
if(hit%500 == 0) cout << "\r" << std::left << Form("Event %d, Hit %d",ev,hit) << flush;
// Get the photon energy and determine the pulse amplitude based on that energy
photonEnergy = 1e6*sqrt( pow(Px->at(hit),2) + pow(Py->at(hit),2) + pow(Pz->at(hit),2) ); //eV
wavelength = 1240./photonEnergy; // nm
pulseIntegral = g->Eval( wavelength )*gain*photonEnergy*1.6022e-19; // mC
amplitude = outputImpedance*pulseIntegral/(1e-9*timeBinWidth*unitPulseIntegral); // mV
// Add the pulse to the pulses vector for each tile
pulses[rodNo->at(hit)/16].push_back( new TF1( Form("ev%d",ev), PMTpulseFunction, 0., timeWindow, 2) );
pulses[rodNo->at(hit)/16].back()->SetParameters( amplitude, time->at(hit) + triggerDelay);
}//end hit loop
for(int tile = 0; tile < 16; tile++){
// Add dark current pulses
int nPulses = rnd.Poisson( aveNpulses );
for(int i = 0; i < nPulses; i++){
pulses[tile].push_back(new TF1( Form("darkCurrent%d",i), PMTpulseFunction, 0., timeWindow, 2) );
pulses[tile].back()->SetParameters( rnd.Gaus( avePulseAmp, 2*avePulseAmp), rnd.Uniform(-80,timeWindow) );
}
// Evaluate the TF1s for each time bin + gaussian noise
for(int bin = 0; bin < 1024; bin ++){
//waveforms[tile]->push_back(DRS4_response->Eval(evaluateSum( pulses[tile], timeBinWidth*bin ) + rnd.Gaus(mean,sigma) ));
waveforms[tile]->push_back(evaluateSum( pulses[tile], timeBinWidth*bin ) + rnd.Gaus(mean,sigma) );
}//end waveform loop
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 objects, clear vectors, reset histograms
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();
if(DRAW)h[tile]->Reset();
}//end tile loop
}//end event loop
cout << endl;
}else if( theMiddleWay ){
// Not implemented yet
}else if( theEasyWay ){
// Not implemented yet
}//end the easy way
outFile->Write();
outFile->Close();
delete outFile; //Deletes histograms
inFile->Close();
delete inFile;
for(int tile = 0; tile < 16; tile++){
delete waveforms[tile];
}
if(!DRAW) delete c;
}
TGraph* getPMTdata(string fileName){
vector<double> x, y;
ifstream file( Form("data/%s", fileName.c_str() ) );
if( !file.is_open() ){
file.open( 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;
}
/*
* Generate timing data to be used by JZCaPA
* Argument is digitization speed in GHz
*/
void generateTimingData( float digiSpeed = 5.0 ){
ofstream outputFile( Form("output/MC_%3.1fGHz.txt", digiSpeed ) );
outputFile << Form("%3.1fGHz simulated timing data", digiSpeed) << endl;
outputFile << " Module1 Module2 Module3 Module4 Module5" << endl;
float binWidth = 1./digiSpeed;
for(int i = 0; i < 1024; i++){
outputFile << Form("%8.3f, %8.3f, %8.3f, %8.3f, %8.3f", i*binWidth, i*binWidth, i*binWidth, i*binWidth, i*binWidth) << endl;
}
}
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;
}
bool checkFileExistencein(const string& filename){
ifstream f(filename.c_str());
return f.is_open();
}
void getFile(string filename, /*out*/ ifstream& file){
const bool file_exists = checkFileExistencein(Form("data/%s", filename.c_str() ) );
if (!file_exists) {
cout << "File in " << Form("data/%s", filename.c_str() ) << " not found." << endl;
cout << "Search " << filename << endl;
bool file_exists = checkFileExistencein(filename);
if (file_exists) {
file.open(filename.c_str());
cout << filename << " is found in the same level" << endl;
}
}
}
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