DSOPlotMaker.C

Go to the documentation of this file.

//C++ INCLUDES
#include <iostream>
#include <fstream>
#include <sstream>
#include <stdlib.h>
#include <vector>
#include <stdint.h>

//ROOT INCLUDES
#include "TROOT.h"
#include "TSystem.h"
#include "TApplication.h"
#include "TFile.h"
#include "TTree.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "TString.h"
#include "TRegexp.h"

//..............................................................................
int GetDetID(std::string detector)
{
  int detID;
  if (detector == "NearDet") {
    detID = 1;
  } else if (detector == "FarDet") {
    detID = 2;
  } else if (detector == "TestBeam") {
    detID = 5;
  } else {
    detID = 0;
    std::cerr << "Invalid detector chosen" << std::endl;
  }
  
  return detID;
}

//..............................................................................
int GetNumberOfDiblocks(std::string detector)
{
  int nDB;
  if (detector == "NearDet") {
    nDB = 4;
  } else if (detector == "FarDet") {
    nDB = 14;
  } else if (detector == "TestBeam") {
    nDB = 1;
  } else {
    nDB = 0;
    std::cerr << "Invalid detector chosen" << std::endl;
  }
  
  return nDB;
}

//..............................................................................
/// Open the FEB file and parse the contents
std::vector<std::vector<int>> ParseFEBList(std::string FEBListPath, std::string dsoScanDir, std::string detector)
{
  
  std::cout << FEBListPath << std::endl;
  std::vector<std::vector<int>> FEBList;
  
  std::ifstream FEBFile(FEBListPath);
  
  // Check that the file opens
  if (!FEBFile)
  {
    std::cout << "The file could not be opened\n";
    exit(1);
  }
   
  std::stringstream FEBFileStream;
  FEBFileStream << FEBFile.rdbuf();

  FEBFile.close();
   
  // Loop over FEBs in file
  while (FEBFileStream){    
    int aDiblock;
    int aDCM;
    int anFEB;
    
    std::cout << "Found " << aDiblock << " " << aDCM << " " << anFEB << std::endl;
    
    // Reads the text file into these int variables 
    FEBFileStream >> aDiblock >> aDCM >> anFEB;
    
    std::vector<int> thisFEB = {aDiblock, aDCM, anFEB};
    FEBList.push_back(thisFEB);
  }
  return FEBList;
}

//..............................................................................
void MakePlots(std::string detector, std::string dsoFile, std::string plotOutDir, int aDB, int aDCM, int aFEB)
{
  int detID = GetDetID(detector);
  
  // Specify ROOT file where we'll save this stuff
  char rootFile[1000];
  sprintf(rootFile, "%s/DSO_apd-%d-%02d-%02d-%02d.root", plotOutDir.c_str(), detID, aDB, aDCM, aFEB);
  
  int thefeb;                         //The FEB number stored in the branch
  int pix;                            //The pixel number stored in the branch
  int ph;                             //The raw pulse height
  uint32_t time;                      //The time word
  int thenum;                         //The number of samples of a given pix and feb
  int num[32]={0};                    //Count of samples for a given pixel
  std::vector<int32_t> rawph[32];     //Raw pulse height vectora
  std::vector<uint32_t> rawticks[32]; //Time sorting vector

  int ret_status; //File status
  
  // Set up color scheme to match that of HW Watch List
  // Array of 32 reasonable colors for pixel level histos
  int hahacolor[33];
  hahacolor[0]  = 632;      //red
  hahacolor[1]  = 920 + 2;  //gray
  hahacolor[2]  = 632 - 7;  //red
  hahacolor[3]  = 900 + 5;  //pink
  hahacolor[4]  = 400 - 6;  //yellow
  hahacolor[5]  = 616;      //mgnt
  hahacolor[6]  = 616 + 2;  //mgnt
  hahacolor[7]  = 616 - 6;  //mgnt
  hahacolor[8]  = 800 + 10; //orange
  hahacolor[9]  = 880 - 7;  //violet
  hahacolor[10] = 600;      //blue
  hahacolor[11] = 600 + 2;  //blue
  hahacolor[12] = 600 - 3;  //blue
  hahacolor[13] = 600 + 3;  //blue
  hahacolor[14] = 600 - 6;  //blue
  hahacolor[15] = 860 + 7;  //azure
  hahacolor[16] = 432 + 1;  //cyan
  hahacolor[17] = 880 - 5;  //violet
  hahacolor[18] = 432 - 6;  //cyan
  hahacolor[19] = 416;      //green
  hahacolor[20] = 416 + 3;  //green
  hahacolor[21] = 416 - 6;  //green
  hahacolor[22] = 416 + 2;  //green
  hahacolor[23] = 820 - 5;  //spring
  hahacolor[24] = 820 - 7;  //spring
  hahacolor[25] = 820 + 4;  //spring
  hahacolor[26] = 400 + 2;  //yellow
  hahacolor[27] = 400 + 3;  //yellow
  hahacolor[28] = 800 + 7;  //orange
  hahacolor[29] = 632 + 2;  //red
  hahacolor[30] = 860 + 7;  //azure
  hahacolor[31] = 900 + 6;  //pink
  hahacolor[32] = 800 + 4;  //orange
  
  // Set up histograms
  TH1F* hadcpix[32];
  TH1F* hadctpix[32];
  TH1F* hdsotpix[32];
  TH1F* hfftpix[32];
  TCanvas* cpix[32];
  TH2F* hthresholdmap = new TH2F("hthreshold", "", 8, -0.5, 7.5, 4, -0.5, 3.5);
  TH2F* hthreshold2map = new TH2F("hthreshold2", "", 8, -0.5, 7.5, 4, -0.5, 3.5); 
  
  // Open ROOT file for writing
  TFile file(rootFile, "RECREATE");
  
  // Open DSO file
  FILE* infile;
  infile = fopen(dsoFile.c_str(), "r");
  for (;;) {
      ret_status = fscanf(infile, "%d, %d, %u, %u, %*x, %*x, %*x\n", &thefeb, &pix, &ph, &time);
      if (ret_status<0) break;
      rawph[pix].push_back(ph); 
      rawticks[pix].push_back(time);
      num[pix]++;
  }
  
  //pix has to exist if the file existed so let's use it to find the step
  int step = 0;
  step = rawticks[pix].at(1) - rawticks[pix].at(0);
  if (rawticks[pix].at(2) == (rawticks[pix].at(1)+step))
    std::cout << "Consistent step found to be " << step << " ticks" << std::endl;
  else step = 1;
  
  float width = 15.625e-9*step; 
  float fftmax = 1.0/width;
  for (int p=0;p<32;p++) {
    thenum = num[p]; 
    std::ostringstream sstrd;
    sstrd <<  aDB << "-" << aDCM << "-" << aFEB << ":" << p;
    hadcpix[p]  = new TH1F(("hadcpix" + (sstrd.str())).c_str(), ";DCS", 300, -150, 150);
    hadctpix[p] = new TH1F(("hadctpix" + (sstrd.str())).c_str(), ";sample;DCS", thenum, -0.5 * width, (thenum - 0.5) * width);
    hdsotpix[p] = new TH1F(("hdsotpix" + (sstrd.str())).c_str(), ";sample;DSO", thenum, -0.5 * width, (thenum - 0.5) * width);
    hfftpix[p]  = new TH1F(("hfftpix" + (sstrd.str())).c_str(), ";frequency [Hz]", thenum, 0.0, fftmax);
    //hfftpix[p] = 0;
    cpix[p] = new TCanvas(("cpix" + (sstrd.str())).c_str(), "", 400, 500);

    for (int j=0; j<thenum-3; j++) {
      hadcpix[p]->Fill(rawph[p][j + 3] - rawph[p][j]);
      hadctpix[p]->Fill(j * width,(rawph[p][j + 3] - rawph[p][j]));
      hdsotpix[p]->Fill(j * width,(rawph[p][j]));
    }
    
    TF1 *mygaus = new TF1("mygaus", "gaus");
    
    int status1 = int(hadcpix[p]->Fit(mygaus));
    int fthreshold = 0.0;
    
    if( status1 && hadcpix[p]->GetRMS() < mygaus->GetParameter(2))
      fthreshold = (int)4.0 * mygaus->GetParameter(2);
    else
      fthreshold = (int)4.0 * hadcpix[p]->GetRMS();
      
    delete mygaus;
    
    int mapi = p%8;
    int mapj = (int)((p - mapi) / 8);
    
    hthresholdmap->Fill(mapi, mapj, fthreshold);
    hadcpix[p]->SetLineColor(hahacolor[p]);
    hadctpix[p]->SetLineColor(hahacolor[p]);
    hdsotpix[p]->SetLineColor(hahacolor[p]);
    //TVirtualFFT::SetTransform(0);
    hdsotpix[p]->FFT(hfftpix[p], "MAG");
    hfftpix[p]->Scale(1.0 / sqrt(thenum));
    //hfftpix[p]->SetAxisRange(2000.0,0.5*fftmax); 
    hfftpix[p]->SetAxisRange(fftmax / 1000.0, 0.5 * fftmax); // low freq cutoff for 1000 samples
    hfftpix[p]->SetLineColor(hahacolor[p]);

    hadcpix[p]->GetXaxis()->SetLabelSize(0.06);
    hadcpix[p]->GetYaxis()->SetLabelSize(0.06);
    hadctpix[p]->GetXaxis()->SetLabelSize(0.06);
    hadctpix[p]->GetYaxis()->SetLabelSize(0.06);
    hdsotpix[p]->GetXaxis()->SetLabelSize(0.06);
    hdsotpix[p]->GetYaxis()->SetLabelSize(0.06);
    hfftpix[p]->GetXaxis()->SetLabelSize(0.06);
    hfftpix[p]->GetYaxis()->SetLabelSize(0.06);
  }
  
  file.cd();
  hthresholdmap->SetMinimum(0);
  hthresholdmap->SetMaximum(100);
  hthresholdmap->SetMarkerSize(3);
  hthresholdmap->SetStats(0);
  hthresholdmap->Write(); 
  hthreshold2map->SetMinimum(0);
  hthreshold2map->SetMaximum(100);
  hthreshold2map->SetMarkerSize(3);
  hthreshold2map->SetStats(0);
  hthreshold2map->Write(); 
  for(int j = 0; j < 32; j++){
    hadctpix[j]->Write();
    hdsotpix[j]->Write();
    hadcpix[j]->Write();
    hfftpix[j]->Write();
  }
  
  // int nDB = GetNumberOfDiblocks(detector);
  
  // Print png files
  for(int pix = 0; pix < 32; ++pix){
    TCanvas* cadcpix  = new TCanvas;
    hadcpix[pix]->Draw();
    cadcpix->SetLogy();
    char outfilename[1000];
    sprintf(outfilename, "%s/DIBLOCK%d/noisyAPD_dcm-%d-%d-%d-%d_hadcpix_%d.png", plotOutDir.c_str(), aDB, detID, aDB, aDCM, aFEB, pix);
    cadcpix->Print(outfilename);
    
    TCanvas* cadctpix = new TCanvas;
    hadctpix[pix]->Draw("hist");
    sprintf(outfilename, "%s/DIBLOCK%d/noisyAPD_dcm-%d-%d-%d-%d_hadctpix_%d.png", plotOutDir.c_str(), aDB, detID, aDB, aDCM, aFEB, pix);
    cadctpix->Print(outfilename);
    
    TCanvas* cdsotpix = new TCanvas;
    hdsotpix[pix]->Draw("hist");
    sprintf(outfilename, "%s/DIBLOCK%d/noisyAPD_dcm-%d-%d-%d-%d_hdsotpix_%d.png", plotOutDir.c_str(), aDB, detID, aDB, aDCM, aFEB, pix);
    cdsotpix->Print(outfilename);
    
    TCanvas* cfftpix  = new TCanvas;
    hfftpix[pix]->Draw("hist");
    sprintf(outfilename, "%s/DIBLOCK%d/noisyAPD_dcm-%d-%d-%d-%d_hfftpix_%d.png", plotOutDir.c_str(), aDB, detID, aDB, aDCM, aFEB, pix);
    cfftpix->Print(outfilename);
  }
}

//..............................................................................
void FEBLoop(std::string detector, std::string dsoScanDir, std::vector<std::vector<int>> FEBlist, std::string plotOutDir)
{
  int detID = GetDetID(detector);
  
  // Loop over files
  for(unsigned int FEBIdx = 0; FEBIdx < FEBlist.size(); ++FEBIdx){
    int aDB  = FEBlist.at(FEBIdx).at(0);
    int aDCM = FEBlist.at(FEBIdx).at(1);
    int aFEB = FEBlist.at(FEBIdx).at(2);
    
    // Check to see if the data exists
    char dsoFile[1000];
    sprintf(dsoFile, "%s/dcm-%d-%02d-%02d/standardVoltage/standardVoltage-FEB%02d.dat", dsoScanDir.c_str(), detID, aDB, aDCM, aFEB);
    std::cout << dsoFile << std::endl;
    
    ifstream File(dsoFile);
    
    if(File.good()){
      std::cout << "We are now using Diblock: " << aDB << " DCM: " << aDCM << " FEB: " << aFEB << std::endl;
      MakePlots(detector, dsoFile, plotOutDir, aDB, aDCM, aFEB);
    }
    
    if(!File){
      std::cout << "We cannot use Diblock: " << aDB << " DCM: " << aDCM << " FEB: " << aFEB << std::endl;
    }
  }
}

//..............................................................................
void DSOPlotMaker(std::string detector, std::string dsoScanDir, std::string FEBListPath, std::string plotOutDir)
{
  // Determine which FEBs are in the list
  std::vector<std::vector<int>> FEBlist = ParseFEBList(FEBListPath, dsoScanDir, detector);
  std::cout << "FEB list size = " << FEBlist.size() << std::endl;
  
  // Loop over files from FEB list and generate plots
  FEBLoop(detector, dsoScanDir, FEBlist, plotOutDir);
}