Based pmt pulses on quantum efficiency and single PE pulses instead of cathode radiant sensitivity

outputDigitization.cc

@@ -344,19 +344,18 @@ void generateWaveform(string fileName, bool DRAW){
   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 );
+  float unitPulseIntegral = testPulse->Integral( 0.0, timeWindow, 1e-4 ); //charge for 1mV pulse
 
   ////////////////////// 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) );
+  TGraph* g = getPMTdata( Form("model%dQE.txt",6091) );
   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 electronPerCoulomb = 6.24e18;
   float aveChargePerPulse = GetAverageCurrent ( g, gain ); // mC
   float aveNpulsePerNanoSec = 1e-9*darkCurrent/aveChargePerPulse; // 1/ns
   float aveNpulses = timeWindow*aveNpulsePerNanoSec; // unitless
@@ -421,13 +420,9 @@ void generateWaveform(string fileName, bool DRAW){
 	outTree->Branch("lastStepZ",&LastStepZf);
 
 	vector< vector< float >* > waveforms(16,0);
-	vector< vector< float >* > heights(16,0);
-	vector< vector< float >* > charges(16,0);
 	for(int tile = 0; tile < 16; tile++){
 		waveforms[tile] = new vector< float >;
 		outTree->Branch( Form("RawSignal%d",tile), &waveforms[tile] );
-		outTree->Branch( Form("pulseHeight%d",tile), &heights[tile] );
-		outTree->Branch( Form("pulseCharge%d",tile), &charges[tile] );
 		if(DRAW) h[tile] = new TH1D( Form("tile%dwaveform",tile), Form("Tile %d Simulated Waveform;time (ns);Amplitude (mV)",tile), nSamples, 0, timeWindow);
 	}
 
@@ -438,7 +433,10 @@ void generateWaveform(string fileName, bool DRAW){
   // //////////////////// The hard way //////////////////////
 	if(theHardWay){
 
-		double wavelength, pulseIntegral, amplitude, photonEnergy;
+		double wavelength, pulseIntegral, photonEnergy;
+    // double onePEamp = outputImpedance/(1e-9*timeBinWidth*unitPulseIntegral); // The amplitude of a single photon pulse in mV
+    double onePEamp = (gain/electronPerCoulomb)/unitPulseIntegral; // V=IR The amplitude of a single photon pulse in mV
+    cout << "onePEamp = " << onePEamp << endl;
 		vector< vector< TF1* > > pulses;
 		pulses.resize(16);
 
@@ -455,23 +453,28 @@ void generateWaveform(string fileName, bool DRAW){
 
 				if(hit%500 == 0) cout << "\r" << std::left << Form("Event %d, Hit %d",eventNo,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
-
-        charges[rodNo->at(hit)/16]->push_back(pulseIntegral);
-        heights[rodNo->at(hit)/16]->push_back(amplitude);
-
-        // Add the pulse to the pulses vector for each tile
-				pulses[rodNo->at(hit)/16].push_back( new TF1( Form("ev%d",eventNo), PMTpulseFunction, 0., timeWindow, 2) );
-				pulses[rodNo->at(hit)/16].back()->SetParameters( amplitude, time->at(hit) + triggerDelay);
+        //The photon has a chance to be detected based on the quantum efficiency of the PMT
+        //roll the dice here
+        if(rnd.Rndm() < g->Eval( wavelength ) ){
+
+          // Add the pulse to the pulses vector for each tile
+  				pulses[rodNo->at(hit)/16].push_back( new TF1( Form("ev%d",eventNo), PMTpulseFunction, 0., timeWindow, 2) );
+  				pulses[rodNo->at(hit)/16].back()->SetParameters( onePEamp, time->at(hit) + triggerDelay);
+  				// pulses[rodNo->at(hit)/16].push_back( new TF1( Form("ev%d",eventNo), "[A]*ROOT::Math::lognormal_pdf(x,[m],[s],[x0])", 0., timeWindow, 2) );
+  				// pulses[rodNo->at(hit)/16].back()->SetParameters( onePEamp, time->at(hit) + triggerDelay);
+        }
 
 			}//end hit loop
 
+      cout << endl << endl;
+
 			for(int tile = 0; tile < 16; tile++){
 
+        if( tile%4 == 0 ) cout << endl;
+        cout << pulses[tile].size() << ", ";
+
         // Add dark current pulses
         int nPulses = rnd.Poisson( aveNpulses );
       	for(int i = 0; i < nPulses; i++){
@@ -605,7 +608,7 @@ void generateWaveform(string fileName, bool DRAW){
 }
 
 
-double GetAverageCurrent( TGraph* PMTcathodeRadSens, double gain ){
+double GetAverageCurrent( TGraph* PMTquantumEff, double gain ){
   TFile* f = new TFile("data/Energy_dist.root","read");
   TH1D* energyHist = (TH1D*)f->Get("Energy dist");
 
@@ -619,7 +622,7 @@ double GetAverageCurrent( TGraph* PMTcathodeRadSens, double gain ){
     photonEnergy = 1e6*energyHist->GetBinCenter(bin); //eV
 
     // current is fraction of photons in this bin * the single photon current at this energy/wavelength
-    current += energyHist->GetBinContent(bin) * PMTcathodeRadSens->Eval( 1240./photonEnergy )*gain*photonEnergy*1.6022e-19; // mC
+    current += energyHist->GetBinContent(bin) * PMTquantumEff->Eval( 1240./photonEnergy )*gain; // mC
   }
 
   f->Close();