Started the ZDC digitizer project

.gitignore

+target1.root
+target2.root
+target3.root
+target4.root
+target5.root
+test.txt
+testlight.root
+Out_testlight.root
+testtest.txt
+trial.txt
+generateSample.cc

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