// global variables
TFile *f = 0;
TTree *T = 0;

int run(const char *fname = "Feb24-12-57-28_16bit.root")
{
if (!(fname && *fname)) return; // just a precaution

if (T) {delete T; T = 0;} // make sure you delete it, if it already exists
if (f) {delete f; f = 0;} // make sure you close the file, if it is opened

//f = new TFile("Feb15-12-16-20_16bit.root", "READ");
f = new TFile(fname, "READ");
if (!f) return; // just a precaution

T = (TTree *)f->FindObjectAny("T");
if (!T) return; // just a precaution

TBranch *b_voltagesI = T->GetBranch("voltages_chan1");
TBranch *b_voltagesQ = T->GetBranch("voltages_chan2");
TBranch *b_timestamps = T->GetBranch("timestamps");

// To check receiver at room temperature with terminator at input:
MakeRoomTemperaturePlot(b_voltagesI, b_voltagesQ);
MakeCryoTemperaturePlot(b_voltagesI, b_voltagesQ);

#if 0 /* 0 or 1 */
if (T) {delete T; T = 0;} // make sure you delete it, if it already exists
if (f) {delete f; f = 0;} // make sure you close the file, if it is opened
#endif /* 0 or 1 */

return 0;

}

void MakeCryoTemperaturePlot(TBranch *b_voltagesI, TBranch *b_voltagesQ)
{
TCanvas *c = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cc");
if (c) delete c; // make sure you delete it, if it already exists
c = new TCanvas("cc","Cryo Temperature",800,1000);
c->Divide(2,3);

c->cd(3);
Double_t video_offset1 = GetVideoOffset(b_voltagesI);
c->cd(4);
Double_t video_offset2 = GetVideoOffset(b_voltagesQ);

c->cd(5);
voltages_sqhisto(b_voltagesI, video_offset1);
c->cd(6);
voltages_sqhisto(b_voltagesQ, video_offset2);

c->cd(1);
Double_t voltage_RMSI = GetVoltageRMS(b_voltagesI, video_offset1, "Real Noise Power");
c->cd(2);
Double_t voltage_RMSQ = GetVoltageRMS(b_voltagesQ, video_offset2, "Imaginary Noise Power");

Double_t real_noise_power = GetNoisePower(voltage_RMSI);
Double_t imaginary_noise_power = GetNoisePower(voltage_RMSQ);

Double_t CryoTemperatureI = GetCryoTemperature(real_noise_power);
c->cd(1); 
TText l;
l.SetTextSize(0.04);
char buffer[50];
int n = sprintf(buffer, "cryo temperature = %.1f K", CryoTemperatureI);
l.DrawText(0.15,700., buffer);

Double_t CryoTemperatureQ = GetCryoTemperature(imaginary_noise_power);
c->cd(2); 
TText l;
l.SetTextSize(0.04);
char buffer[50];
int n = sprintf(buffer, "cryo temperature = %.1f K", CryoTemperatureQ);
l.DrawText(0.15,700., buffer);

}





void MakeRoomTemperaturePlot(TBranch *b_voltagesI, TBranch *b_voltagesQ)
{
TCanvas *c = (TCanvas *)gROOT->GetListOfCanvases()->FindObject("cr");
if (c) delete c; // make sure you delete it, if it already exists
c = new TCanvas("cr","Room Temperature",800,1000);
c->Divide(2,3);

c->cd(3);
Double_t video_offset1 = GetVideoOffset(b_voltagesI);
c->cd(4);
Double_t video_offset2 = GetVideoOffset(b_voltagesQ);

c->cd(5);
voltages_sqhisto(b_voltagesI, video_offset1);
c->cd(6);
voltages_sqhisto(b_voltagesQ, video_offset2);

c->cd(1);
Double_t voltage_RMSI = GetVoltageRMS(b_voltagesI, video_offset1, "Real Noise Power");
c->cd(2);
Double_t voltage_RMSQ = GetVoltageRMS(b_voltagesQ, video_offset2, "Imaginary Noise Power");

Double_t real_noise_power = GetNoisePower(voltage_RMSI);
Double_t imaginary_noise_power = GetNoisePower(voltage_RMSQ);

Double_t RoomTemperatureI = GetRoomTemperature(real_noise_power);
c->cd(1); 
TText l;
l.SetTextSize(0.04);
char buffer[50];
int n = sprintf(buffer, "room temperature = %.0f K", RoomTemperatureI);
l.DrawText(0.15,700., buffer);

Double_t RoomTemperatureQ = GetRoomTemperature(imaginary_noise_power);
c->cd(2); 
TText l;
l.SetTextSize(0.04);
char buffer[50];
int n = sprintf(buffer, "room temperature = %.0f K", RoomTemperatureQ);
l.DrawText(0.15,700., buffer);

}




double GetNoisePower(double voltage_RMS)
{
// Noise power is RMS of histogram.  In dBm/Hz:  
Double_t noise_power=0.;
// With 4 MHz LPF at output, voltage gain of 5X:
noise_power = 10.*log10(voltage_RMS*voltage_RMS*1.e3/25./50./4.e6);
printf("Output noise power is %f dBm/Hz \n", noise_power);  
return noise_power;
}



double GetCryoTemperature(double noise_power)
{

static const int ncomps = 12;

Double_t HEMT_noise_temp = 22.0; // Kelvin

//                                 HEMT      mx1     fif1    ampif1  fif2    mx2     fif2   ampif2   ampif2  fif2    mx3    lpf1    
Double_t gains_db[ncomps]=        { 30.0,    -10.,    0.0,    60.0,   0.0,   -6.0,   -1.0,   34.0,    0.0,   -1.0,   -8.5,   0.0};
Double_t noise_figures_db[ncomps]={ 0.01,      10.5,   0.0,    1.5,    0.0,    8.5,    1.0,    2.0,    0.0,    1.0,    10.0,  1.0};
Double_t gains_lin[ncomps] = {};
Double_t noise_figures_lin[ncomps] = {};
Double_t total_gains_db[ncomps] = {};
Double_t total_noise_figures_lin[ncomps] = {};
Double_t total_noise_figures_db[ncomps] = {};
Double_t gain_cascade[ncomps] = {};


for (Int_t i=0; i<ncomps; i++)
  {
  gains_lin[i]=pow(10.,gains_db[i]/10.);
  noise_figures_lin[i]=pow(10.,noise_figures_db[i]/10.);  
  }

for (Int_t i=0; i<ncomps; i++)
  {
  for (Int_t j=0; j<i+1; j++)
    {
    total_gains_db[i] += gains_db[j];
    }
  }


total_noise_figures_lin[0] = noise_figures_lin[0];
total_noise_figures_db[0] = 10.*log10(total_noise_figures_lin[0]);

for (Int_t i=1; i<ncomps; i++)
  {
  gain_cascade[i] = 1.0;
  for (Int_t j=0; j<i; j++)
    gain_cascade[i] *= gains_lin[j];
  total_noise_figures_lin[i] = total_noise_figures_lin[i-1] + (noise_figures_lin[i]-1.) / gain_cascade[i];  
  total_noise_figures_db[i] = 10.*log10(total_noise_figures_lin[i]);
      
  }  


/*
for (Int_t i=0; i<ncomps; i++)
  {
  printf("\n gains_lin %d is %f\n", i, gains_lin[i]);
  printf("noise_figures_lin %d is %f\n", i, noise_figures_lin[i]);
  printf("gain cascade %d is %f\n", i, gain_cascade[i]);
  printf("total noise figures lin %d is %f\n", i, total_noise_figures_lin[i]);
  printf("total noise figures db %d is %f\n", i, total_noise_figures_db[i]);
  printf("total gains db %d is %f\n\n\n", i, total_gains_db[i]);
  }
*/


printf("input power (dBm/Hz) is %f\n", noise_power-total_gains_db[ncomps-1]-total_noise_figures_db[ncomps-1]);
Double_t CryoTemperature = pow(10.,(noise_power-total_gains_db[ncomps-1]-total_noise_figures_db[ncomps-1])/10.)/1.38e-20;
printf("Cryo temperature is %f\n", CryoTemperature);

return CryoTemperature;


}








double GetRoomTemperature(double noise_power)
{

static const int ncomps = 11;


//                                mx1     fif1    ampif1  fif2    mx2     fif2   ampif2   ampif2  fif2    mx3    lpf1    
Double_t gains_db[ncomps]=       {-10.,    0.0,    60.0,   0.0,   -6.0,   -1.0,   34.0,    0.0,   -1.0,   -8.5,   0.0};
Double_t noise_figures_db[ncomps]={10.5,   0.0,    1.5,    0.0,    8.5,    1.0,    2.0,    0.0,    1.0,    10.0,  1.0};
Double_t gains_lin[ncomps] = {};
Double_t noise_figures_lin[ncomps] = {};
Double_t total_gains_db[ncomps] = {};
Double_t total_noise_figures_lin[ncomps] = {};
Double_t total_noise_figures_db[ncomps] = {};
Double_t gain_cascade[ncomps] = {};


for (Int_t i=0; i<ncomps; i++)
  {
  gains_lin[i]=pow(10.,gains_db[i]/10.);
  noise_figures_lin[i]=pow(10.,noise_figures_db[i]/10.);  
  }

for (Int_t i=0; i<ncomps; i++)
  {
  for (Int_t j=0; j<i+1; j++)
    {
    total_gains_db[i] += gains_db[j];
    }
  }

total_noise_figures_lin[0] = noise_figures_lin[0];
total_noise_figures_db[0] = 10.*log10(total_noise_figures_lin[0]);

for (Int_t i=1; i<ncomps; i++)
  {
  gain_cascade[i] = 1.0;
  for (Int_t j=0; j<i; j++)
    gain_cascade[i] *= gains_lin[j];
  total_noise_figures_lin[i] = total_noise_figures_lin[i-1] + (noise_figures_lin[i]-1.) / gain_cascade[i];  
  total_noise_figures_db[i] = 10.*log10(total_noise_figures_lin[i]);
      
  }  

/*
for (Int_t i=0; i<ncomps; i++)
  {
  printf("\ngains_lin %d is %f\n", i, gains_lin[i]);
  printf("noise_figures_lin %d is %f\n", i, noise_figures_lin[i]);
  printf("gain cascade %d is %f\n", i, gain_cascade[i]);
  printf("total noise figures lin %d is %f\n", i, total_noise_figures_lin[i]);
  printf("total noise figures db %d is %f\n", i, total_noise_figures_db[i]);
  printf("total gains db %d is %f\n\n", i, total_gains_db[i]);
  }
*/


printf("input power (dBm/Hz) is %f\n", noise_power-total_gains_db[ncomps-1]-total_noise_figures_db[ncomps-1]);
Double_t RoomTemperature = pow(10.,(noise_power-total_gains_db[ncomps-1]-total_noise_figures_db[ncomps-1])/10.)/1.38e-20;
printf("Room temperature is %f\n", RoomTemperature);

return RoomTemperature;


}







void voltages_gaussian_sq()
{

TH1F *hgaussian_sq = new TH1F("hgaussian_sq", "Fake Gaussian, squared; Volts^2; Counts", 100, -3., 10.);
TRandom3 r;
Double_t x=0.;
for (Int_t i=0; i<10000; i++)
  {
  x=r.Gaus(0.,0.5);
  hgaussian_sq->Fill(x*x);
  }
hgaussian_sq->DrawCopy(); 
delete hgaussian_sq;
}




void voltages_gaussian()
{

TH1F *hgaussian = new TH1F("hgaussian", "Fake Gaussian; Volts; Counts", 100, -3., 3.);
TRandom3 r;
for (Int_t i=0; i<10000; i++)
  hgaussian->Fill(r.Gaus(0.,0.5));
hgaussian->DrawCopy(); 
delete hgaussian;
}

void voltages_sqhisto(TBranch *b_voltages, double video_offset)
{
Int_t nentries = (Int_t)b_voltages->GetEntries();
Double_t voltages[1024] = {};
b_voltages->SetAddress(voltages);

TH1F *hvoltages_sq = new TH1F("hvoltages_sq", "O-scope data minus offset, squared; Volts^2; Counts", 100, -0.01, 0.05);

for (Int_t i=0;i<nentries;i++) 
  {
  b_voltages->GetEntry(i);
  // fill the histogram with the voltages
  for (Int_t j=0; j<1024; j++)
    hvoltages_sq->Fill((voltages[j]-video_offset)*(voltages[j]-video_offset));  
  }
  hvoltages_sq->DrawCopy();

delete hvoltages_sq;
}



double GetVoltageRMS(TBranch *b_voltages, double video_offset, char *htitle)
{
Int_t nentries = (Int_t)b_voltages->GetEntries();
Double_t voltages[1024] = {};
b_voltages->SetAddress(voltages);

TH1F *hvoltages_m_offset = new TH1F("hvoltages_m_offset", (TString(htitle) += ";Volts;Counts").Data(), 100, -0.5, 0.5);

for (Int_t i=0;i<nentries;i++) 
  {
  b_voltages->GetEntry(i);
  // fill the histogram with the voltages
  for (Int_t j=0; j<1024; j++)
    hvoltages_m_offset->Fill(voltages[j]-video_offset);  
  }


hvoltages_m_offset->DrawCopy();


Double_t RMS = hvoltages_m_offset->GetRMS();
delete hvoltages_m_offset;

return RMS;



}





double GetVideoOffset(TBranch *b_voltages)
{
Int_t nentries = (Int_t)b_voltages->GetEntries();
Double_t voltages[1024] = {};
b_voltages->SetAddress(voltages);
Double_t histomean=0.;
Double_t noise_power;

TH1F *hvoltages = new TH1F("hvoltages", "Raw O-scope data; Volts; Counts", 100, -0.5, 0.5);

for (Int_t i=0;i<nentries;i++) 
  {
  b_voltages->GetEntry(i);
  // fill the histogram with the voltages
  for (Int_t j=0; j<1024; j++)
    hvoltages->Fill(voltages[j]);  
  }
  hvoltages->DrawCopy();


Double_t offset = hvoltages->GetMean();
delete hvoltages;


return offset;


}