outputDigitization.cc

/////////////////////////////////////////
//	I've only tested cutZDCoutput so far
//	Run with:
//		root[0] .L outputDigitizer.cc
//		root[1] cutZDCoutput(testlight.root,5)
//
//	It works with the transmission spectra from
//	all 7 rods 1,2,3,4,5,6,9 (old numbering scheme, I'll change soon)
//
//	Todo:
//		1) Get quantum efficiency curves from Hamamatsu data sheets
//		2) Work those into PMTcuts
//		3) Get rid of readCSV since QE curves will have x and y points
//
//
//	\author YaBoi
/////////////////////////////////////////
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TRandom.h"
#include <TTree.h>
#include <vector>
#include <fstream>


using namespace std;

TFile* inFile;
TTree* inTree;
TFile* outFile;
TTree* outTree;

//Create fiber vectors
vector<double> *Xf=0, *Yf=0, *Zf=0, *X0=0, *Y0=0, *Z0=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, *y=0, *z=0,*x_0=0, *y_0=0, *z_0=0, *px=0, *py=0, *pz=0, *Energy=0;
vector<int> *EventNo=0;

void SetBranches(string fileName, bool input = false);
vector< double > readCSV(string fileName);

/*
 * 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, int rodNum){
	SetBranches(fileName);
	int nBins = 3649;
	int EventNo;
	double energy,X,Y,Z,Px,Py,Pz;
	TRandom rnd;
	ofstream outputFile( Form("%s_Out.txt",fileName.c_str() ) );

	//// Get the transmission data
	////
	vector< double > transData;
	string str;
	ifstream file( Form("rod%dlongitudinal.txt", rodNum ) );
	if(!file.is_open()){
	    cout << fileName << " didn't open" << 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);
	float trans, rand;

	for(int bin = 0; bin < nBins; bin++){
		hTrans->SetBinContent( bin, transData[bin] );
	}


	//Get the number of events in the file
	int nEvents = inTree->GetEntries();

	//Find out how many events contain photons
	int realNevents = 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) );
			// Get transmission % from data
			trans = hTrans->GetBinContent( hTrans->GetBin( energy ) );

			// Give the photon a random chance to make it weighted
			// by transmission %
			if( rnd.Rndm() < trans ){

				// If the photon is transmitted, add it to the vector
				X = Xf->at(k);
				Z = Zf->at(k);
				x->push_back( X );
				z->push_back( Z );

				//We only 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("%5d of %5d photons cut in event %4d",nCut,nPhotons,ev) << endl;
		}
		outTree->Fill();

	}//end loop over events
	inFile->Close();
	outFile->Close();
}


/*
 * Apply quantum efficiency cuts to the output of the lightguide MC
 */
void PMTcuts(string fileName, int PMTmodel = 6091){
	SetBranches(true);

	//// Get the quantum efficiency data
	////
	ifstream file("data/model%dQE.csv");
	if( !file.is_open() ){
		cout << "Quantum efficiency data file didn't open... exiting" << endl;
		return;
	}
	vector<double> wavelength, efficiency;
	string str;
	size_t comma;

	while( !file.eof() ){
		getline(file,str);
		if( str.length() ){
			comma = str.find_first_of(",");

			wavelength.push_back( atof( str.substr(0,comma) ) );
			efficiency.push_back( atof( str.substr( comma, str.length()-comma ) ) );
		}
	}

	// Load that data into a histogram so we can use the GetBin function
	TH1F* hQE = new TH1F( "TransData", "TransData", wavelength.size(), wavelength.front(),  wavelength.back());
	for(int bin = 0; bin < wavelength.size(); bin++0){
		hQE->SetBinContent(bin,efficiency[bin]);
	}


	//// begin loop over events
	////
	TRandom rnd;
	float energy, trans;
	int nEvents = inTree->GetEntries();
	for (int ev = 0; ev < nEvents ; ev++){
		inTree->GetEntry(ev);

		int Kf = Xf->size();
		for (int k=0; k<Kf; 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) );
			// Get transmission % from data
			trans = hQE->GetBinContent( hQE->GetBin( energy ) );

			// Give the photon a random chance to make it past the
			// Quantum Efficiency check
			if( rnd.Rndm() < trans ){
					// If the PMT detects the photon, add it to the vector
					X = Xf->at(k);
					Z = Zf->at(k);
					x->push_back( X );
					z->push_back( Z );

					px->push_back( Px );
					py->push_back( Py );
					pz->push_back( Pz );
					Energy->push_back( energy );
			}//end cuts
		}//end vector loop
	}//end event looop
}//end PMTcuts


/*
 * Set branch addresses of the global TTrees
 * input = true means save the origin position data for the particle
 */
void SetBranches(string fileName, bool input = false){
		inFile = new TFile(fileName.c_str());
	 	inTree = (TTree*)inFile->Get("ZDCtree");

		inTree->SetBranchAddress("X",&Xf);
		inTree->SetBranchAddress("Z",&Zf);
		//inTree->SetBranchAddress("X0",&X0);
		//inTree->SetBranchAddress("Z0",&Z0);
		inTree->SetBranchAddress("Px",&Pxf);
		inTree->SetBranchAddress("Py",&Pyf);
		inTree->SetBranchAddress("Pz",&Pzf);
		//inTree->SetBranchAddress("Pid",&PIDf);
		//inTree->SetBranchAddress("Energy",&Energyf);
		//inTree->SetBranchAddress("ID",&IDf);
		//inTree->SetBranchAddress("EventNo",&EventNof);
		//intree->SetBranchAddress("NCherenkovs",&NCherenkovs);

		outFile = new TFile( Form("Out_%s",fileName.c_str() ), "RECREATE");
		outTree = new TTree( "tree", "tree" );

		if(input){
			outTree->Branch("X0",&x_0);
			outTree->Branch("Y0",&y_0);
			outTree->Branch("Z0",&z_0);
		}
		outTree->Branch("X",&x);
		outTree->Branch("Y",&y);
		outTree->Branch("Z",&z);
		outTree->Branch("Px",&px);
		outTree->Branch("Py",&py);
		outTree->Branch("Pz",&pz);
		outTree->Branch("Energy",&Energy);
		outTree->Branch("EventNo",&EventNo);
}//end SetBranches