APDPlots.C

Go to the documentation of this file.

#include <TSQLServer.h>
#include <TSQLResult.h>
#include <iostream.h>
#include <fstream>
#include <iomanip.h>
#include <sstream>
#include <string>
#include <cstring>
#include <algorithm>

void APDPlots(int diblock, int dcm, int port) {

  gROOT->Reset();
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(1);
  gStyle->SetOptTitle(1);
  gStyle->SetPalette(1);
  gStyle->SetFillColor(0);
  gStyle->SetPadColor(0);
  gStyle->SetCanvasColor(0);
  gStyle->SetStatColor(0);
  gStyle->SetTitleColor(1);
  gStyle->SetTitleSize(0.06);
  gStyle->SetPadBorderMode(0);
  gStyle->SetPadBottomMargin(0.15);
  gStyle->SetPadLeftMargin(0.12);
  gStyle->SetPadRightMargin(0.09);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetCanvasBorderMode(0);

  TDatime da(2010,01,01,1,00,00);
  gStyle->SetTimeOffset(da.Convert());

  TDatime *TNow = new TDatime;//.TDatime();
  int XNow=TNow->Convert();
  unsigned int X2DaysAgo    =XNow - 2*60*60*24;
  unsigned int X4WeekAgo    =XNow - 60*60*24*7*4;
  TDatime *T2d = new TDatime;
  T2d->Set(X2DaysAgo,false);
  TDatime *T4w = new TDatime;
  T4w->Set(X4WeekAgo,false);

  TNow->Print();
  cout<<XNow<<"\t"<<X2DaysAgo<<"\t"<<60*60*24*2<<"\t"<<X2DaysAgo-XNow<<endl;
  cout<<"DATES: " <<TNow->AsSQLString()<< " " << T2d->AsSQLString()<< " " << T4w->AsSQLString() << endl;

  std::stringstream dbstream;

  dbstream << "select channel_info.pe,channel_info.pecorr,channel_info.rate_low,channel_info.rate_hi,"
           << " channel_info.nhits,channel_info.nn_rate,daq_info.channel,daq_info.pixel,"
           << " apd_header.time_last,temp_info.temp_ave,temp_info.temp_min,temp_info.temp_max,apd_header.num_temp,"
           << " channel_info.pesq,channel_info.pecorrsq,apd_header.time_elapsed"
           << " from channel_info,daq_info,apd_header,temp_info where "
           << " daq_info.dcm="<<dcm
           << " and temp_info.dcm="<<dcm
           << " and daq_info.diblock="<<diblock
           << " and temp_info.diblock="<<diblock
           << " and daq_info.port="<<port
           << " and temp_info.port="<<port
           << " and channel_info.channel=daq_info.channel"
           << " and apd_header.run=channel_info.run"
           << " and apd_header.subrun=channel_info.subrun"
           << " and apd_header.run=temp_info.run"
           << " and apd_header.subrun=temp_info.subrun"
           << " and apd_header.time_first> '"<<T2d->AsSQLString()<<"' "
           << " and apd_header.time_last<'"<<TNow->AsSQLString()<<" ' "
           << " order by daq_info.pixel";

  std::string dbstring = dbstream.str();

  TH1F *pe_pix_2d = new TH1F("pe_pix_2d", "pe_hist;pixel number;# p.e. per hit (last 48h)", 32, 0, 31);
  TH1F *pec_pix_2d = new TH1F("pec_pix_2d", "pecorr_hist;pixel number; # corrected p.e. per hit (if calibrated, last 48h)", 32, 0, 31);
  TH1F *nn_pix_2d = new TH1F("nn_pix_2d", "nn_hist;pixel number;nearest neighbor variable (last 48h)", 32, 0, 31);
  TH1F *rh_pix_2d = new TH1F("rh_pix_2d", "ratehi_hist;pixel number; number hits above 100 p.e. per second (last 48h)", 32, 0, 31);
  TH1F *rl_pix_2d = new TH1F("rl_pix_2d", "ratelo_hist;pixel number; number of hits below 100 p.e. per second (last 48h)", 32, 0, 31);

  TH1F *pe_apd_2d = new TH1F("pe_apd_2d", "pe_hist;(48 - #) hours ago; # p.e. per hit (whole APD)", 48, 1, 48);
  TH1F *pec_apd_2d = new TH1F("pec_apd_2d", "pecorr_hist;(48 - #) hours ago; # corrected p.e. per hit (if calibrated, whole APD)", 48, 1, 48);
  TH1F *nn_apd_2d = new TH1F("nn_apd_2d", "nn_hist;(48 - #) hours ago;nearest neighbor variable (whole APD)", 48, 1, 48);
  TH1F *rh_apd_2d = new TH1F("rh_apd_2d", "ratehi_hist;(48 - #) hours ago; # hits above 100 p.e. per second (whole APD)", 48, 1, 48);
  TH1F *rl_apd_2d = new TH1F("rl_apd_2d", "ratelo_hist;(48 - #) hours ago; # hits below 100 p.e. per second (whole APD)", 48, 1, 48);
  TH1F *apd_temp = new TH1F("apd_temp", "apd_temp;# hours ago; average temperature of APD (degrees C)", 48, 1, 48);

  char * dbcstr;
  dbcstr = new char [dbstring.size()+1];
  strcpy(dbcstr, dbstring.c_str());

  //password removed
  TSQLServer *db = TSQLServer::Connect("pgsql://131.215.116.171", "", "");

  std::cout << "Connected? " << db->IsConnected() << std::endl;

  TSQLStatement *stmt = db->Statement(dbcstr,1000);

  std::cout << " about to send statement: " << dbcstr << std::endl;

  if(stmt->Process()){
    // store result of statement in buffer
    stmt->StoreResult();

    //    std::cout << " we got some results I think! " << std::endl;

    double pe_pix[34] = {0.};
    double pecorr_pix[34] = {0.};
    double pesq_pix[34] = {0.};
    double pecorrsq_pix[34] = {0.};
    double nn_pix[34] = {0.};
    double rate_hi_pix[34] = {0.};
    double rate_low_pix[34] = {0.};
    int nhit_pix[34] = {0.};
    int count_pix[34] = {0.};
    int telapse_pix[34] = {0.};

    double pe_2d[50] = {0.};
    double pecorr_2d[50] = {0.};
    double pesq_2d[50] = {0.};
    double pecorrsq_2d[50] = {0.};
    double nn_2d[50] = {0.};
    double rate_hi_2d[50] = {0.};
    double rate_low_2d[50] = {0.};
    int nhit_2d[50] = {0.};
    int count_2d[50] = {0.};
    double tave[50] = {0.};
    double tmin[50] = {0.};
    double tmax[50] = {0.};
    int telapse_2d[50] = {0.};

    int tpix, tnh, tchan, ttn, telapse;
    double tpe, tpec, tnn, trh, trl, tpesq, tpecorrsq;
    double ttave, ttmin, ttmax;

    int counter = 0;

    while (stmt->NextResultRow())
      {
        counter++;
        //        std::cout << " reading in row number: " << counter << std::endl;
        // Extract rows one after another

        tpe = stmt->GetDouble(0);
        tpec = stmt->GetDouble(1);
        trl = stmt->GetDouble(2);
        trh = stmt->GetDouble(3);
        tnh = stmt->GetInt(4);
        tnn = stmt->GetDouble(5);
        tchan = stmt->GetInt(6);
        tpix = stmt->GetInt(7);
        const char* time_last = stmt->GetString(8);
        ttave = stmt->GetDouble(9);
        ttmin = stmt->GetDouble(10);
        ttmax = stmt->GetDouble(11);
        ttn = stmt->GetInt(12);
        tpesq = stmt->GetDouble(13);
        tpecorrsq = stmt->GetDouble(14);
        telapse = stmt->GetInt(15);

        pe_pix[tpix+1] += tnh*tpe;
        pecorr_pix[tpix+1] += tnh*tpec;
        pesq_pix[tpix+1] += tnh*tpesq;
        pecorrsq_pix[tpix+1] += tnh*tpecorrsq;
        nn_pix[tpix+1] += tnn*telapse;
        rate_hi_pix[tpix+1] += trh*telapse;
        rate_low_pix[tpix+1] += trl*telapse;
        nhit_pix[tpix+1] += tnh;
        count_pix[tpix+1] ++;
        telapse_pix[tpix+1] += telapse;

        TDatime *ltime = new TDatime;
        ltime->Set(time_last);
        unsigned int tlint = ltime->Convert();

        int tbin = (tlint-X2DaysAgo) / 3600;
        if(tbin == 48) tbin = 47;
        tbin = tbin + 1;

        //        std::cout << " " << tpe << " " << tpix << " " << tnh << " " << ttmin << " " << ttmax << " " << tbin << std::endl;

        pe_2d[tbin] += tpe*tnh;
        pecorr_2d[tbin] += tpec*tnh;
        pesq_2d[tbin] += tpesq*tnh;
        pecorrsq_2d[tbin] += tpecorrsq*tnh;
        nn_2d[tbin] += tnn*telapse;
        rate_hi_2d[tbin] += trh*telapse;
        rate_low_2d[tbin] += trl*telapse;
        nhit_2d[tbin] += tnh;
        count_2d[tbin] += ttn;
        telapse_2d[tbin] += telapse;

        tave[tbin] += ttn*ttave;
        if(ttmin==0) ttmin = ttave - (ttmax-ttave);  // sometimes we get no temperature data, not the same as bad data, cheat a reasonable value
        if(tmin[tbin]==0) tmin[tbin]=ttmin;
        else tmin[tbin] = min(ttmin,tmin[tbin]);
        if(tmax[tbin]==0) tmax[tbin]=ttmax;
        else tmax[tbin] = max(ttmax,tmax[tbin]);

      }

    double error1[34] = {0.};
    double error2[34] = {0.};
    double error3[34] = {0.};
    double error4[34] = {0.};
    double error5[34] = {0.};

    for(unsigned int n = 1; n < 33; ++n){

      if(nhit_pix[n]!=0) {

        error4[n] = sqrt(rate_hi_pix[n]) / telapse_pix[n];
        error5[n] = sqrt(rate_low_pix[n]) / telapse_pix[n];

        pe_pix[n] = pe_pix[n]/nhit_pix[n];
        pecorr_pix[n] = pecorr_pix[n]/nhit_pix[n];
        pesq_pix[n] = pesq_pix[n]/nhit_pix[n];
        pecorrsq_pix[n] = pecorrsq_pix[n]/nhit_pix[n];
        nn_pix[n] = nn_pix[n]/telapse_pix[n];
        rate_hi_pix[n] = rate_hi_pix[n]/telapse_pix[n];
        rate_low_pix[n] = rate_low_pix[n]/telapse_pix[n];

        error3[n] = ( sqrt(nhit_pix[n]) / nhit_pix[n] ) * nn_pix[n];

        if(pesq_pix[n]==0) {
          error1[n] = ( 1 / sqrt(nhit_pix[n]) ) * pe_pix[n];
          error2[n] = ( 1 / sqrt(nhit_pix[n]) ) * pecorr_pix[n];
        }
        else if(pesq_pix[n]>0) {
          error1[n] = sqrt( fabs(pesq_pix[n] - (pe_pix[n]*pe_pix[n])) ) * 1.5 / sqrt(nhit_pix[n]);
          error2[n] = sqrt( fabs(pecorrsq_pix[n] - (pecorr_pix[n]*pecorr_pix[n])) * 1.5 / sqrt(nhit_pix[n]) );
        }
        else {
          std::cout << "Negative squares???  Something is wrong!" << std::endl;
        }

      }
    }

    pe_pix_2d->SetContent(pe_pix);
    pec_pix_2d->SetContent(pecorr_pix);
    nn_pix_2d->SetContent(nn_pix);
    rh_pix_2d->SetContent(rate_hi_pix);
    rl_pix_2d->SetContent(rate_low_pix);

    pe_pix_2d->SetError(error1);
    pec_pix_2d->SetError(error2);
    nn_pix_2d->SetError(error3);
    rh_pix_2d->SetError(error4);
    rl_pix_2d->SetError(error5);

    double err1[50] = {0.};
    double err2[50] = {0.};
    double err3[50] = {0.};
    double err4[50] = {0.};
    double err5[50] = {0.};
    double err6[50] = {0.};

    for(unsigned int n = 1; n < 49; ++n){

      if(nhit_2d[n]!=0) {

        err4[n] = sqrt(rate_hi_2d[n]) / telapse_2d[n];
        err5[n] = sqrt(rate_low_2d[n]) / telapse_2d[n];

        pe_2d[n] = pe_2d[n]/nhit_2d[n];
        pecorr_2d[n] = pecorr_2d[n]/nhit_2d[n];
        pesq_2d[n] = pesq_2d[n]/nhit_2d[n];
        pecorrsq_2d[n] = pecorrsq_2d[n]/nhit_2d[n];
        nn_2d[n] = nn_2d[n]/nhit_2d[n];
        rate_hi_2d[n] = rate_hi_2d[n]/telapse_2d[n];
        rate_low_2d[n] = rate_low_2d[n]/telapse_2d[n];

        err3[n] = ( sqrt(nhit_2d[n]) / nhit_2d[n] ) * nn_2d[n];

        if(pesq_2d[n]==0) {
          err1[n] = ( 1 / sqrt(nhit_2d[n]) ) * pe_2d[n];
          err2[n] = ( 1 / sqrt(nhit_2d[n]) ) * pecorr_2d[n];
        }
        else if(pesq_2d[n]>0) {
          err1[n] = sqrt( fabs(pesq_2d[n] - (pe_2d[n]*pe_2d[n])) ) * 1.5 / sqrt(nhit_2d[n]);
          err2[n] = sqrt( fabs(pecorrsq_2d[n] - (pecorr_2d[n]*pecorr_2d[n])) * 1.5 / sqrt(nhit_2d[n]) );
        }
        else {
          std::cout << "Negative squares???  Something is wrong!" << std::endl;
        }

      }
      if(count_2d[n]!=0) {
        tave[n] = tave[n]/count_2d[n];
        err6[n] = (tmax[n]-tmin[n])/2.0;
      }
    }

    pe_apd_2d->SetContent(pe_2d);
    pec_apd_2d->SetContent(pecorr_2d);
    nn_apd_2d->SetContent(nn_2d);
    rh_apd_2d->SetContent(rate_hi_2d);
    rl_apd_2d->SetContent(rate_low_2d);
    apd_temp->SetContent(tave);

    pe_apd_2d->SetError(err1);
    pec_apd_2d->SetError(err2);
    nn_apd_2d->SetError(err3);
    rh_apd_2d->SetError(err4);
    rl_apd_2d->SetError(err5);
    apd_temp->SetError(err6);

    delete stmt;

  }

  std::stringstream dbstream2;

  dbstream2 << "select pixel,channel from daq_info where "
           << " dcm="<<dcm
           << " and diblock="<<diblock
           << " and port="<<port
           << " order by pixel";

  std::string dbstring2 = dbstream2.str();

  char * dbcstr2;
  dbcstr2 = new char [dbstring2.size()+1];
  strcpy(dbcstr2, dbstring2.c_str());

  TSQLStatement *stmt2 = db->Statement(dbcstr2,1000);

  if(stmt2->Process()){
    // store result of statement in buffer
    stmt2->StoreResult();

    while (stmt2->NextResultRow()) {

        const char* pix = stmt2->GetString(0);
        const char* chann = stmt2->GetString(1);

        cout << "pixel: " << pix << " channel: " << chann << endl;

    }
  }

  delete db;

}