Finally got units right after adding 1/ns to amplitude calculation

outputDigitization.cc

@@ -12,7 +12,7 @@ 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);
+void addDarkCurrent(vector< TF1* >& vec, float timeInterval, float timeBinWidth);
 void addRandNoise(vector< float >*  vec, float timeInterval);
 void addTimingJitter(vector< float >* vec, float timeInterval);
 
@@ -310,6 +310,9 @@ void generateWaveform(string fileName, bool DRAW){
 
   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);
@@ -340,7 +343,7 @@ void generateWaveform(string fileName, bool DRAW){
 	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), 1024, 0, 204.8);
+		if(DRAW) h[tile] = new TH1D( Form("tile%dwaveform",tile), Form("Tile %d Simulated Waveform;time (ns);Amplitude (mV)",tile), nSamples, 0, timeWindow);
 	}
 
 
@@ -357,17 +360,17 @@ void generateWaveform(string fileName, bool DRAW){
 		for(int ev = 0; ev < nEntries; ev++){
 			tree->GetEntry(ev);
 			int nHits = Px->size();
-			if(nHits == 0) continue;
+			// 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) );
+        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 )*photonEnergy; // mC
-				amplitude = 1e-3*inputImpedance*pulseIntegral/101.358l; // mV
+				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, 0, 204.8, 1) );
+				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);
 
@@ -375,7 +378,7 @@ void generateWaveform(string fileName, bool DRAW){
 
 			for(int tile = 0; tile < 16; tile++){
 				// Add dark current here to minimize calls to evaluate sum
-				addDarkCurrent(pulses[tile], 204.8);
+				addDarkCurrent(pulses[tile], 204.8, 1./sampleFrequency);
 
 				for(int bin = 0; bin < 1024; bin ++){
 					waveforms[tile]->push_back(evaluateSum( pulses[tile], 0.2*bin ) );
@@ -384,12 +387,11 @@ void generateWaveform(string fileName, bool DRAW){
 				addRandNoise(waveforms[tile], 204.8);
 
 
-				if(DRAW && tile==10){
+				if(DRAW){
 					for(int bin = 0; bin < 1024; bin ++){
 						h[tile]->SetBinContent(bin, waveforms[tile]->at(bin));
 					}//end waveform loop
-					// c->cd(tile+1);
-					c->cd(0);
+					c->cd(tile+1);
 					h[tile]->Draw();
 					h[tile]->SetAxisRange(-10,750,"Y");
 				}
@@ -444,15 +446,15 @@ void addRandNoise(vector< float >* vec, float timeInterval){
  * 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){
+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 = 1e-3*21.7475; // A/W aka mC/J
+	float aveRadSens = 21.7475; // mA/W aka mC/J
   float gain = 1e6; // unitless
-	float aveChargePerPulse = 1e-3*avePhotonEnergy*aveRadSens*gain; // Coulomb
-	float darkCurrent = 2e-9; // Amps aka C/s
+	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
 
@@ -460,7 +462,7 @@ void addDarkCurrent(vector< TF1* >& pulses, float timeInterval){
 	float buffer = 80; // ns
 	timeInterval += buffer;
 	float aveNpulses = timeInterval*aveNpulsePerNanoSec; // unitless. Probably want to randomize this a bit
-	float avePulseAmp = 1e3*aveChargePerPulse*gain*inputImpedance/101.3581; // mV
+	float avePulseAmp = aveChargePerPulse*inputImpedance/(timeBinWidth*101.3581); // mV
 
   TRandom2 rnd;
   rnd.SetSeed(gSystem->Now());