pmt_plot.cxx

Go to the documentation of this file.

//#define pmts_cxx
//#include "pmts.h"

#include <TROOT.h>
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include "Plex/PlexHandle.h"
#include <iostream>
#include <TGraphErrors.h>
#include <TF1.h>
#include <TSystem.h>
#include <TStyle.h>
#include <stdio.h>
//#include <fcntl.h> /* For STDIN_FILENO */
//#include <sys/select.h> /* For select() */
#include "mysql.h"
#include <TPostScript.h>
#include "TMySQLServer.h"
#include "TSQLResult.h"
#include "TSQLRow.h"
#include "Calibrator/Calibrator.h"


// N.B.: 
//loon [0] .include /usr/include/mysql
//loon [1] gSystem->Load("libMySQL")    

using namespace std;

//void pmts::Loop()
//{
//}

TCanvas* GetCanvas(const char* name)
{
  TCanvas* ctmp = (TCanvas*) gROOT->FindObject(name);
  if(!ctmp) ctmp = new TCanvas(name,name);
  return ctmp;
}

/*
void Pause()
{
  bool done = false;
  do{
    gSystem->ProcessEvents();
    fd_set read_fds;
    FD_ZERO(&read_fds);
    FD_SET( STDIN_FILENO, &read_fds );
    struct timeval timeout;
    timeout.tv_sec=0;
    timeout.tv_usec=1000;
    int res = select(STDIN_FILENO+1, &read_fds, NULL, NULL, &timeout);
    if(res>0) {
      getchar();
      done=true;
    }
  } while(!done);
}
*/

static std::vector<TObject*> cleanup;
static TGraph* gTemperature1= 0;
static TGraph* gTemperature2= 0;
TMySQLServer* server = 0;

void MakeTemperaturePlot(TMySQLServer* server,
                         Detector::Detector_t det,
                         UInt_t start, UInt_t stop)
{ 
  const char* query;
  if(det == Detector::kFar) {
    query = Form("select "
                 "UNIX_TIMESTAMP(DCS_ENV_FARVLD.TIMESTART), "
                 "DCS_ENV_FAR.FAR_TEMPERATURE1, "
                 "DCS_ENV_FAR.FAR_TEMPERATURE2 "
                 "from DCS_ENV_FARVLD, DCS_ENV_FAR "
                 "where (DCS_ENV_FAR.SEQNO=DCS_ENV_FARVLD.SEQNO) "
                 "and (UNIX_TIMESTAMP(DCS_ENV_FARVLD.TIMESTART)>%d) "
                 "and (UNIX_TIMESTAMP(DCS_ENV_FARVLD.TIMESTART)<%d) ",
                 start,stop
                 );
  } else {
    //return;
    query = Form("select "
                 "UNIX_TIMESTAMP(DCS_ENV_NEARVLD.TIMESTART), "
                 "DCS_ENV_NEAR.NEAR_THERM_V3b, "
                 "DCS_ENV_NEAR.NEAR_THERM_V3t "
                 "from DCS_ENV_NEARVLD, DCS_ENV_NEAR "
                 "where (DCS_ENV_NEAR.SEQNO=DCS_ENV_NEARVLD.SEQNO) "
                 "and (UNIX_TIMESTAMP(DCS_ENV_NEARVLD.TIMESTART)>%d) "
                 "and (UNIX_TIMESTAMP(DCS_ENV_NEARVLD.TIMESTART)<%d) ",
                 start,stop
                 );
  }
  
  cout << query << endl;

  TSQLResult* result = server->Query(query);
  if(!result) return;

  if(gTemperature1) delete gTemperature1;
  if(gTemperature2) delete gTemperature2;

  cout << "Building temperature plot with " 
       << "  rows=" << result->GetRowCount() << endl;

  gTemperature1 = new TGraph(result->GetRowCount());
  gTemperature2 = new TGraph(result->GetRowCount());

  double* x1 = gTemperature1->GetX();
  double* y1 = gTemperature1->GetY();
  double* x2 = gTemperature2->GetX();
  double* y2 = gTemperature2->GetY();
  
  int p = 0;

  TSQLRow* row = 0;
  while( ( row = result->Next() ) ) {

    UInt_t time;
    Float_t farenheit1;
    Float_t farenheit2;

    sscanf(row->GetField(0),"%u",&time);
    sscanf(row->GetField(1),"%f",&farenheit1);
    sscanf(row->GetField(2),"%f",&farenheit2);

    x1[p] = time;
    x2[p] = time;
    y1[p] = (farenheit1-32.0)/1.8;
    y2[p] = (farenheit2-32.0)/1.8;

    p++;   
  }
 
}


int PlotPmt(//pmts& tree,
             TMySQLServer* server,
             PlexPixelSpotId tube) 
{
  for(UInt_t i=0;i<cleanup.size();i++)  delete cleanup[i];
  cleanup.clear();

  UInt_t encoded = tube.GetUniquePmtEncodedValue();

  PlexHandle plex(VldContext(tube.GetDetector(), SimFlag::kData, VldTimeStamp()));
  PlexPixelSpotId psid = tube;
  psid.SetPixel(1);
  psid.SetSpot(1);
  RawChannelId rcid = plex.GetRawChannelId(psid);

  char title[1000];
  sprintf(title,"%s (%u) %s",tube.AsString("t"),encoded,rcid.AsString("e"));

  cout << title << endl;
 
  // Find number of points:
  Int_t point = 0;

  const char* query = Form("select "
                           "UNIX_TIMESTAMP(CALPMTDRIFTVLD.TIMESTART) "
                           ",CALPMTDRIFT.DRIFT "
                           ",CALPMTDRIFT.STATERROR "
                           ",CALPMTDRIFT.SYSERROR "
                           ",CALPMTDRIFT.SPOTSUSED "
                           ",CALPMTDRIFT.SPOTSAVAIL "
                           ",CALPMTDRIFT.CRUDEDRIFT "
                           ",CALPMTDRIFT.MINDRIFT "
                           ",CALPMTDRIFT.MAXDRIFT "
                           "from CALPMTDRIFTVLD, CALPMTDRIFT where "
                           "(CALPMTDRIFT.PMT = %u) "
                           "and (CALPMTDRIFT.SEQNO=CALPMTDRIFTVLD.SEQNO) "
                           "and (CALPMTDRIFTVLD.TASK=1) "
                           //"and (CALPMTDRIFTVLD.TIMESTART < '2005-07-31 19:00:04') " 
                           //"and (CALPMTDRIFTVLD.CREATIONDATE >= '2005-09-02 21:00:05') " 
                           //"and (CALPMTDRIFTVLD.TIMESTART > '2005-06-25 00:00:00') " 
                           "and (CALPMTDRIFTVLD.DETECTORMASK=%d) order by CALPMTDRIFTVLD.TIMESTART;",
                           encoded,
                           tube.GetDetector()                                   
                           );

  cout << query << endl;

  TSQLResult* result = server->Query(query);
  if(!result) return 0;

  cout << tube.AsString("t") 
       << "  " << encoded 
       << "  rows=" << result->GetRowCount() << endl;


  int npoints = result->GetRowCount();
  
  TGraphErrors* gFinal = new TGraphErrors(npoints);  cleanup.push_back(gFinal);
  TGraphErrors* gRaw   = new TGraphErrors(npoints);  cleanup.push_back(gRaw);
  TGraph* gHigh        = new TGraph(npoints);  cleanup.push_back(gHigh);
  TGraph* gLow         = new TGraph(npoints);  cleanup.push_back(gLow);
  TGraph* gNgood       = new TGraph(npoints);  cleanup.push_back(gNgood);
  TGraph* gNspot       = new TGraph(npoints);  cleanup.push_back(gNspot);
  

  TSQLRow* row = 0;
  while( ( row = result->Next() ) ) {

    if(point>npoints) break;
    UInt_t time;
    Float_t drift;
    Float_t staterr;
    Float_t syserr;
    Float_t ngood;
    Float_t nspot;
    Float_t rawdrift;
    Float_t mindrift = -999;
    Float_t maxdrift = -999;

    int irow=0;
    sscanf(row->GetField(irow++),"%u",&time);
    sscanf(row->GetField(irow++),"%f",&drift);
    sscanf(row->GetField(irow++),"%f",&staterr);
    sscanf(row->GetField(irow++),"%f",&syserr);
    sscanf(row->GetField(irow++),"%f",&ngood);
    sscanf(row->GetField(irow++),"%f",&nspot);
    sscanf(row->GetField(irow++),"%f",&rawdrift);
    sscanf(row->GetField(irow++),"%f",&mindrift);
    sscanf(row->GetField(irow++),"%f",&maxdrift);
    
    
    gFinal->SetPoint     (point,time, drift);
    gFinal->SetPointError(point,0, staterr);
    gRaw->SetPoint       (point,time, rawdrift);
    gRaw->SetPointError  (point,0, 0);
    gHigh->SetPoint      (point,time, mindrift);
    gLow->SetPoint       (point,time, maxdrift);
    gNgood->SetPoint     (point,time, ngood);
    gNspot->SetPoint     (point,time, nspot);
    
    //VldContext cx(det, SimFlag::kData, VldTimeStamp(time,0));
    //Calibrator::Instance().Reset(cx);
    //Float_t temp = Calibrator::Instance().GetTemperature(1);
    //gTemp->SetPoint      (point,time, temp);
    
    point++;   
  }

  if(point==0) return 0;
       
    // Find means of the gains over time.
  // Reject data points close to zero (likely HV failures)
  double finalmean = 0;
  double rawmean = 0;
  double finaln = 0;
  double rawn = 0;
  double* finalX = gFinal->GetX();
  double* finalY = gFinal->GetY();
  double* finalEY = gFinal->GetEY();
  double* rawY   = gRaw->GetY();
  double* rawEY =  gRaw->GetEY();
  for(int i=0;i<point;i++) {
    if(finalY[i]>10.) {
      finalmean += finalY[i];
      finaln+=1.;
    }
    if(rawY[i]>10.) {
      rawmean   += rawY[i];
      rawn += 1.0;
    }
  }
  finalmean /= rawn;
  rawmean /= finaln;


  TGraphErrors* gRelFinal = new TGraphErrors(npoints);  cleanup.push_back(gRelFinal);
  TGraphErrors* gRelRaw   = new TGraphErrors(npoints);  cleanup.push_back(gRelRaw);

  Double_t* relFinalX = gRelFinal->GetX();
  Double_t* relFinalY = gRelFinal->GetY();
  Double_t* relFinalEX= gRelFinal->GetEX();
  Double_t* relFinalEY= gRelFinal->GetEY();
  Double_t* relRawX = gRelRaw->GetX();
  Double_t* relRawY = gRelRaw->GetY();
  Double_t* relRawEX= gRelRaw->GetEX();
  Double_t* relRawEY= gRelRaw->GetEY();
  for(Int_t i=0;i<npoints;i++) {
    relFinalX[i] = finalX[i];
    relFinalY[i] = (finalY[i]-finalmean)/finalmean * 100.;
    relFinalEX[i] = 0;
    relFinalEY[i] = finalEY[i]/finalmean * 100.;

    relRawX[i] = finalX[i];
    relRawY[i] = (rawY[i]-rawmean)/rawmean * 100.;
    relRawEX[i] = 0;
    relRawEY[i] = rawEY[i]/rawmean * 100.;
  }


  TCanvas* cpmt = GetCanvas("cpmt");
  cpmt->cd();
  cpmt->Clear();
  cpmt->Divide(1,3,0.001,0.001,kWhite);
  cpmt->cd(1);

  
  double tlow = gFinal->GetX()[0];
  double thigh = gFinal->GetX()[point-1];

  //const double kMonth = 60.*60.*24.*365./12.;
  //const double tref = VldTimeStamp(2005,01,15,0,0,0,0).GetSec();

  //tlow  = ( (int)((tlow -tref)/kMonth)    )*kMonth + tref;
  //thigh = ( (int)((thigh-tref)/kMonth) + 1)*kMonth + tref;

  if(gTemperature1==0) MakeTemperaturePlot(server,
                                           tube.GetDetector(),
                                           (UInt_t)tlow,
                                           (UInt_t)thigh);
  
  double ymax = 150.;
  if(tube.GetDetector()==Detector::kNear) ymax=500.;
  TH2F* hFrame1 = new TH2F("hFrame1",
                           title,
                           100,tlow,thigh,
                           100,0,ymax);
  hFrame1->SetDirectory(0);
  cleanup.push_back(hFrame1);

  hFrame1->GetXaxis()->SetTimeDisplay(1);
  hFrame1->GetXaxis()->SetTimeFormat("%d/%b/%y %Hh");
  //hFrame1->SetLabelSize(0.04,"x");
  hFrame1->SetXTitle("Date");
  hFrame1->SetYTitle("Gain (ADC counts)");
  hFrame1->SetStats(0);
  hFrame1->Draw("");

 
  gRaw->SetLineColor(kRed);
  gRaw->SetMarkerColor(kRed);
  gRaw->SetMarkerStyle(20);
  gRaw->SetMarkerSize(0.45);
  gRaw->Draw("P");

  gFinal->SetMarkerStyle(20);
  gFinal->SetMarkerSize(0.5);
  gFinal->Draw("P");

  gLow->SetMarkerStyle(23);
  gLow->SetMarkerSize(0.6);
  gLow->SetMarkerColor(kBlue);
  gLow->Draw("P");
  gHigh->SetMarkerStyle(24);
  gHigh->SetMarkerSize(0.6);
  gHigh->SetMarkerColor(kBlue);
  gHigh->Draw("P");

  gNspot->SetLineColor(kBlack);
  gNspot->Draw("L");
  
  gNgood->SetLineColor(kGreen);
  gNgood->Draw("L");
  
  if(gTemperature1) {
    gTemperature1->SetLineColor(kBlue);
    gTemperature1->Draw("L");
  }
  

  cpmt->cd(2);
  TH2F* hFrame2 = new TH2F("hFrame2",
                           title,
                           100,tlow,thigh,
                           100,-4,4.);
  hFrame2->SetDirectory(0);
  cleanup.push_back(hFrame2);

  hFrame2->SetStats(0);
  hFrame2->GetXaxis()->SetTimeDisplay(1);
  hFrame2->GetXaxis()->SetTimeFormat("%d/%b/%y %Hh");
  hFrame2->SetLabelSize(0.04,"x");
  hFrame2->SetXTitle("Date");
  hFrame2->SetYTitle("Relative Gain (%)");
  hFrame2->Draw("");


  gRelRaw->SetLineColor(kRed);
  gRelRaw->SetMarkerColor(kRed);
  gRelRaw->SetMarkerStyle(21);
  gRelRaw->SetMarkerSize(0.5);
  gRelRaw->Draw("P");

  gRelFinal->SetMarkerStyle(20);
  gRelFinal->SetMarkerSize(0.5);
  gRelFinal->Draw("P");

  gStyle->SetOptFit(11111);

  TF1 *f1 = new TF1("f1","[0]",0,1e99);
  gRelFinal->Fit(f1,"RQ");
  f1->Draw("l same");


  cpmt->cd(3);
  TH2F* hFrame3 = new TH2F("hFrame3","Temperature",
                           100,tlow,thigh,
                           100,15,30.);
  hFrame3->SetDirectory(0);
  cleanup.push_back(hFrame3);

  hFrame3->SetStats(0);
  hFrame3->GetXaxis()->SetTimeDisplay(1);
  hFrame2->GetXaxis()->SetTimeFormat("%d/%b/%y %Hh");
  hFrame3->SetLabelSize(0.04,"x");
  hFrame3->SetXTitle("Date");
  hFrame3->SetYTitle("Temperature");
  hFrame3->Draw("");

  if(gTemperature2) {
    gTemperature2->SetLineColor(kCyan);
    gTemperature2->SetMarkerColor(kCyan);
    gTemperature2->SetMarkerStyle(20);
    gTemperature2->SetMarkerSize(0.01);
    gTemperature2->Draw("P");
  }

  if(gTemperature1) {
    gTemperature1->SetLineColor(kBlue);
    gTemperature1->SetMarkerColor(kBlue);
    gTemperature1->SetMarkerStyle(20);
    gTemperature1->SetMarkerSize(0.01);
    gTemperature1->Draw("P");
  }

  cout << "Done." << endl;

  cpmt->Update();
  return point;
}


void pmt_plot() 
{
  TCanvas* cpmt = (TCanvas*) gROOT->FindObject("cpmt");
  if(cpmt) delete cpmt;
  cpmt = new TCanvas("cpmt","cpmt",1000,750);
  
  TPostScript* ps = 0;
  ps = new TPostScript("pmt_drift.ps",112);
  gStyle->SetTimeOffset(0);


  // For some reason, I have to do this before I open the TMySQL interface.
  Calibrator::Instance();

  // Load db:
  TString host(getenv("ENV_TSQL_URL") );
  host.ReplaceAll("odbc:","");
  if(!server) {
    //server = new TMySQLServer( "mysql://minos-db1.fnal.gov/offline", "reader","minos_db");
    server = new TMySQLServer( "mysql://localhost/offline", "reader","numi96");
                             // host.Data(),
                             //  getenv("ENV_TSQL_USER"),
                             //  getenv("ENV_TSQL_PSWD")
                              // );
  }

  PlexPixelSpotId tube;
  
  /*
  tube = PlexPixelSpotId(2701264895);
  if(ps) ps->NewPage();
  PlotPmt(server, tube);
  Pause();
  
  tube = PlexPixelSpotId(2701266943);
  if(ps) ps->NewPage();
  PlotPmt(server, tube);
  Pause();
  */

  //Detector::Detector_t det = Detector::kNear;
  Detector::Detector_t det = Detector::kFar;

  
  VldContext cx(det,SimFlag::kData,VldTimeStamp());
  PlexHandle plex(cx);
  const std::vector<PlexPixelSpotId>& tubes = plex.GetAllTubes();

  //std::vector<PlexPixelSpotId> tubes;
  //ElecType::Elec_t elec = ElecType::kVA;
  //for(char ew = 'E'; ew<='W'; ew+='W'-'E')
  //  for(char ul = 'L'; ul<='U'; ul+='U'-'L') 
  //    for(int bay=1;bay<=16;bay++) 
  //    for(int inrack=0;inrack<=9;inrack++)
  //      for(int t=0;t<3;t++) {
  //        tubes.push_back(
  //                        PlexPixelSpotId(det,elec,ew,ul,bay,inrack,t,0,0)
  //                        );
  //      }

  if(ps) ps->NewPage();

  for(UInt_t i=0;i<tubes.size();i++) {
    tube = tubes[i];
    int res = PlotPmt(server, tube);
    if(res) {
      ps->NewPage();
      cpmt->Print(Form("pmt_drifts/%s.gif",tube.AsString("t")));
    }
  }
  
  if(ps) ps->Close();
}




void draw(UInt_t itube)
{
  // Load db:
  TString host(getenv("ENV_TSQL_URL") );
  host.ReplaceAll("odbc:","");
  if(!server) server = new TMySQLServer( host.Data(),
                                         getenv("ENV_TSQL_USER"),
                                         getenv("ENV_TSQL_PSWD")
                                         );

  TCanvas* cpmt = (TCanvas*) gROOT->FindObject("cpmt");
  if(!cpmt) cpmt = new TCanvas("cpmt","cpmt",1000,750);

  PlexPixelSpotId tube(itube);

  cout << "Trying tube " << itube << "  " << tube.AsString() << tube.GetUniquePmtEncodedValue() << endl;
  PlotPmt(server, tube); 
}

---