DSOPlotMaker.C
Go to the documentation of this file.
1 //C++ INCLUDES
2 #include <iostream>
3 #include <fstream>
4 #include <sstream>
5 #include <stdlib.h>
6 #include <vector>
7 #include <stdint.h>
8 
9 //ROOT INCLUDES
10 #include "TROOT.h"
11 #include "TSystem.h"
12 #include "TApplication.h"
13 #include "TFile.h"
14 #include "TTree.h"
15 #include "TF1.h"
16 #include "TH1F.h"
17 #include "TH2F.h"
18 #include "TCanvas.h"
19 #include "TString.h"
20 #include "TRegexp.h"
21 
22 //..............................................................................
24 {
25  int detID;
26  if (detector == "NearDet") {
27  detID = 1;
28  } else if (detector == "FarDet") {
29  detID = 2;
30  } else if (detector == "TestBeam") {
31  detID = 5;
32  } else {
33  detID = 0;
34  std::cerr << "Invalid detector chosen" << std::endl;
35  }
36 
37  return detID;
38 }
39 
40 //..............................................................................
42 {
43  int nDB;
44  if (detector == "NearDet") {
45  nDB = 4;
46  } else if (detector == "FarDet") {
47  nDB = 14;
48  } else if (detector == "TestBeam") {
49  nDB = 1;
50  } else {
51  nDB = 0;
52  std::cerr << "Invalid detector chosen" << std::endl;
53  }
54 
55  return nDB;
56 }
57 
58 //..............................................................................
59 /// Open the FEB file and parse the contents
60 std::vector<std::vector<int>> ParseFEBList(std::string FEBListPath, std::string dsoScanDir, std::string detector)
61 {
62 
63  std::cout << FEBListPath << std::endl;
64  std::vector<std::vector<int>> FEBList;
65 
66  std::ifstream FEBFile(FEBListPath);
67 
68  // Check that the file opens
69  if (!FEBFile)
70  {
71  std::cout << "The file could not be opened\n";
72  exit(1);
73  }
74 
75  std::stringstream FEBFileStream;
76  FEBFileStream << FEBFile.rdbuf();
77 
78  FEBFile.close();
79 
80  // Loop over FEBs in file
81  while (FEBFileStream){
82  int aDiblock;
83  int aDCM;
84  int anFEB;
85 
86  std::cout << "Found " << aDiblock << " " << aDCM << " " << anFEB << std::endl;
87 
88  // Reads the text file into these int variables
89  FEBFileStream >> aDiblock >> aDCM >> anFEB;
90 
91  std::vector<int> thisFEB = {aDiblock, aDCM, anFEB};
92  FEBList.push_back(thisFEB);
93  }
94  return FEBList;
95 }
96 
97 //..............................................................................
98 void MakePlots(std::string detector, std::string dsoFile, std::string plotOutDir, int aDB, int aDCM, int aFEB)
99 {
100  int detID = GetDetID(detector);
101 
102  // Specify ROOT file where we'll save this stuff
103  char rootFile[1000];
104  sprintf(rootFile, "%s/DSO_apd-%d-%02d-%02d-%02d.root", plotOutDir.c_str(), detID, aDB, aDCM, aFEB);
105 
106  int thefeb; //The FEB number stored in the branch
107  int pix; //The pixel number stored in the branch
108  int ph; //The raw pulse height
109  uint32_t time; //The time word
110  int thenum; //The number of samples of a given pix and feb
111  int num[32]={0}; //Count of samples for a given pixel
112  std::vector<int32_t> rawph[32]; //Raw pulse height vectora
113  std::vector<uint32_t> rawticks[32]; //Time sorting vector
114 
115  int ret_status; //File status
116 
117  // Set up color scheme to match that of HW Watch List
118  // Array of 32 reasonable colors for pixel level histos
119  int hahacolor[33];
120  hahacolor[0] = 632; //red
121  hahacolor[1] = 920 + 2; //gray
122  hahacolor[2] = 632 - 7; //red
123  hahacolor[3] = 900 + 5; //pink
124  hahacolor[4] = 400 - 6; //yellow
125  hahacolor[5] = 616; //mgnt
126  hahacolor[6] = 616 + 2; //mgnt
127  hahacolor[7] = 616 - 6; //mgnt
128  hahacolor[8] = 800 + 10; //orange
129  hahacolor[9] = 880 - 7; //violet
130  hahacolor[10] = 600; //blue
131  hahacolor[11] = 600 + 2; //blue
132  hahacolor[12] = 600 - 3; //blue
133  hahacolor[13] = 600 + 3; //blue
134  hahacolor[14] = 600 - 6; //blue
135  hahacolor[15] = 860 + 7; //azure
136  hahacolor[16] = 432 + 1; //cyan
137  hahacolor[17] = 880 - 5; //violet
138  hahacolor[18] = 432 - 6; //cyan
139  hahacolor[19] = 416; //green
140  hahacolor[20] = 416 + 3; //green
141  hahacolor[21] = 416 - 6; //green
142  hahacolor[22] = 416 + 2; //green
143  hahacolor[23] = 820 - 5; //spring
144  hahacolor[24] = 820 - 7; //spring
145  hahacolor[25] = 820 + 4; //spring
146  hahacolor[26] = 400 + 2; //yellow
147  hahacolor[27] = 400 + 3; //yellow
148  hahacolor[28] = 800 + 7; //orange
149  hahacolor[29] = 632 + 2; //red
150  hahacolor[30] = 860 + 7; //azure
151  hahacolor[31] = 900 + 6; //pink
152  hahacolor[32] = 800 + 4; //orange
153 
154  // Set up histograms
155  TH1F* hadcpix[32];
156  TH1F* hadctpix[32];
157  TH1F* hdsotpix[32];
158  TH1F* hfftpix[32];
159  TCanvas* cpix[32];
160  TH2F* hthresholdmap = new TH2F("hthreshold", "", 8, -0.5, 7.5, 4, -0.5, 3.5);
161  TH2F* hthreshold2map = new TH2F("hthreshold2", "", 8, -0.5, 7.5, 4, -0.5, 3.5);
162 
163  // Open ROOT file for writing
164  TFile file(rootFile, "RECREATE");
165 
166  // Open DSO file
167  FILE* infile;
168  infile = fopen(dsoFile.c_str(), "r");
169  for (;;) {
170  ret_status = fscanf(infile, "%d, %d, %u, %u, %*x, %*x, %*x\n", &thefeb, &pix, &ph, &time);
171  if (ret_status<0) break;
172  rawph[pix].push_back(ph);
173  rawticks[pix].push_back(time);
174  num[pix]++;
175  }
176 
177  //pix has to exist if the file existed so let's use it to find the step
178  int step = 0;
179  step = rawticks[pix].at(1) - rawticks[pix].at(0);
180  if (rawticks[pix].at(2) == (rawticks[pix].at(1)+step))
181  std::cout << "Consistent step found to be " << step << " ticks" << std::endl;
182  else step = 1;
183 
184  float width = 15.625e-9*step;
185  float fftmax = 1.0/width;
186  for (int p=0;p<32;p++) {
187  thenum = num[p];
188  std::ostringstream sstrd;
189  sstrd << aDB << "-" << aDCM << "-" << aFEB << ":" << p;
190  hadcpix[p] = new TH1F(("hadcpix" + (sstrd.str())).c_str(), ";DCS", 300, -150, 150);
191  hadctpix[p] = new TH1F(("hadctpix" + (sstrd.str())).c_str(), ";sample;DCS", thenum, -0.5 * width, (thenum - 0.5) * width);
192  hdsotpix[p] = new TH1F(("hdsotpix" + (sstrd.str())).c_str(), ";sample;DSO", thenum, -0.5 * width, (thenum - 0.5) * width);
193  hfftpix[p] = new TH1F(("hfftpix" + (sstrd.str())).c_str(), ";frequency [Hz]", thenum, 0.0, fftmax);
194  //hfftpix[p] = 0;
195  cpix[p] = new TCanvas(("cpix" + (sstrd.str())).c_str(), "", 400, 500);
196 
197  for (int j=0; j<thenum-3; j++) {
198  hadcpix[p]->Fill(rawph[p][j + 3] - rawph[p][j]);
199  hadctpix[p]->Fill(j * width,(rawph[p][j + 3] - rawph[p][j]));
200  hdsotpix[p]->Fill(j * width,(rawph[p][j]));
201  }
202 
203  TF1 *mygaus = new TF1("mygaus", "gaus");
204 
205  int status1 = int(hadcpix[p]->Fit(mygaus));
206  int fthreshold = 0.0;
207 
208  if( status1 && hadcpix[p]->GetRMS() < mygaus->GetParameter(2))
209  fthreshold = (int)4.0 * mygaus->GetParameter(2);
210  else
211  fthreshold = (int)4.0 * hadcpix[p]->GetRMS();
212 
213  delete mygaus;
214 
215  int mapi = p%8;
216  int mapj = (int)((p - mapi) / 8);
217 
218  hthresholdmap->Fill(mapi, mapj, fthreshold);
219  hadcpix[p]->SetLineColor(hahacolor[p]);
220  hadctpix[p]->SetLineColor(hahacolor[p]);
221  hdsotpix[p]->SetLineColor(hahacolor[p]);
222  //TVirtualFFT::SetTransform(0);
223  hdsotpix[p]->FFT(hfftpix[p], "MAG");
224  hfftpix[p]->Scale(1.0 / sqrt(thenum));
225  //hfftpix[p]->SetAxisRange(2000.0,0.5*fftmax);
226  hfftpix[p]->SetAxisRange(fftmax / 1000.0, 0.5 * fftmax); // low freq cutoff for 1000 samples
227  hfftpix[p]->SetLineColor(hahacolor[p]);
228 
229  hadcpix[p]->GetXaxis()->SetLabelSize(0.06);
230  hadcpix[p]->GetYaxis()->SetLabelSize(0.06);
231  hadctpix[p]->GetXaxis()->SetLabelSize(0.06);
232  hadctpix[p]->GetYaxis()->SetLabelSize(0.06);
233  hdsotpix[p]->GetXaxis()->SetLabelSize(0.06);
234  hdsotpix[p]->GetYaxis()->SetLabelSize(0.06);
235  hfftpix[p]->GetXaxis()->SetLabelSize(0.06);
236  hfftpix[p]->GetYaxis()->SetLabelSize(0.06);
237  }
238 
239  file.cd();
240  hthresholdmap->SetMinimum(0);
241  hthresholdmap->SetMaximum(100);
242  hthresholdmap->SetMarkerSize(3);
243  hthresholdmap->SetStats(0);
244  hthresholdmap->Write();
245  hthreshold2map->SetMinimum(0);
246  hthreshold2map->SetMaximum(100);
247  hthreshold2map->SetMarkerSize(3);
248  hthreshold2map->SetStats(0);
249  hthreshold2map->Write();
250  for(int j = 0; j < 32; j++){
251  hadctpix[j]->Write();
252  hdsotpix[j]->Write();
253  hadcpix[j]->Write();
254  hfftpix[j]->Write();
255  }
256 
257  // int nDB = GetNumberOfDiblocks(detector);
258 
259  // Print png files
260  for(int pix = 0; pix < 32; ++pix){
261  TCanvas* cadcpix = new TCanvas;
262  hadcpix[pix]->Draw();
263  cadcpix->SetLogy();
264  char outfilename[1000];
265  sprintf(outfilename, "%s/DIBLOCK%d/noisyAPD_dcm-%d-%d-%d-%d_hadcpix_%d.png", plotOutDir.c_str(), aDB, detID, aDB, aDCM, aFEB, pix);
266  cadcpix->Print(outfilename);
267 
268  TCanvas* cadctpix = new TCanvas;
269  hadctpix[pix]->Draw("hist");
270  sprintf(outfilename, "%s/DIBLOCK%d/noisyAPD_dcm-%d-%d-%d-%d_hadctpix_%d.png", plotOutDir.c_str(), aDB, detID, aDB, aDCM, aFEB, pix);
271  cadctpix->Print(outfilename);
272 
273  TCanvas* cdsotpix = new TCanvas;
274  hdsotpix[pix]->Draw("hist");
275  sprintf(outfilename, "%s/DIBLOCK%d/noisyAPD_dcm-%d-%d-%d-%d_hdsotpix_%d.png", plotOutDir.c_str(), aDB, detID, aDB, aDCM, aFEB, pix);
276  cdsotpix->Print(outfilename);
277 
278  TCanvas* cfftpix = new TCanvas;
279  hfftpix[pix]->Draw("hist");
280  sprintf(outfilename, "%s/DIBLOCK%d/noisyAPD_dcm-%d-%d-%d-%d_hfftpix_%d.png", plotOutDir.c_str(), aDB, detID, aDB, aDCM, aFEB, pix);
281  cfftpix->Print(outfilename);
282  }
283 }
284 
285 //..............................................................................
286 void FEBLoop(std::string detector, std::string dsoScanDir, std::vector<std::vector<int>> FEBlist, std::string plotOutDir)
287 {
288  int detID = GetDetID(detector);
289 
290  // Loop over files
291  for(unsigned int FEBIdx = 0; FEBIdx < FEBlist.size(); ++FEBIdx){
292  int aDB = FEBlist.at(FEBIdx).at(0);
293  int aDCM = FEBlist.at(FEBIdx).at(1);
294  int aFEB = FEBlist.at(FEBIdx).at(2);
295 
296  // Check to see if the data exists
297  char dsoFile[1000];
298  sprintf(dsoFile, "%s/dcm-%d-%02d-%02d/standardVoltage/standardVoltage-FEB%02d.dat", dsoScanDir.c_str(), detID, aDB, aDCM, aFEB);
299  std::cout << dsoFile << std::endl;
300 
301  ifstream File(dsoFile);
302 
303  if(File.good()){
304  std::cout << "We are now using Diblock: " << aDB << " DCM: " << aDCM << " FEB: " << aFEB << std::endl;
305  MakePlots(detector, dsoFile, plotOutDir, aDB, aDCM, aFEB);
306  }
307 
308  if(!File){
309  std::cout << "We cannot use Diblock: " << aDB << " DCM: " << aDCM << " FEB: " << aFEB << std::endl;
310  }
311  }
312 }
313 
314 //..............................................................................
315 void DSOPlotMaker(std::string detector, std::string dsoScanDir, std::string FEBListPath, std::string plotOutDir)
316 {
317  // Determine which FEBs are in the list
318  std::vector<std::vector<int>> FEBlist = ParseFEBList(FEBListPath, dsoScanDir, detector);
319  std::cout << "FEB list size = " << FEBlist.size() << std::endl;
320 
321  // Loop over files from FEB list and generate plots
322  FEBLoop(detector, dsoScanDir, FEBlist, plotOutDir);
323 }
void DSOPlotMaker(std::string detector, std::string dsoScanDir, std::string FEBListPath, std::string plotOutDir)
Definition: DSOPlotMaker.C:315
fVtxDx Fit("f")
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
OStream cerr
Definition: OStream.cxx:7
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
string outfilename
knobs that need extra care
string infile
int GetNumberOfDiblocks(std::string detector)
Definition: DSOPlotMaker.C:41
const double j
Definition: BetheBloch.cxx:29
std::vector< std::vector< int > > ParseFEBList(std::string FEBListPath, std::string dsoScanDir, std::string detector)
Open the FEB file and parse the contents.
Definition: DSOPlotMaker.C:60
OStream cout
Definition: OStream.cxx:6
void FEBLoop(std::string detector, std::string dsoScanDir, std::vector< std::vector< int >> FEBlist, std::string plotOutDir)
Definition: DSOPlotMaker.C:286
void MakePlots(std::string detector, std::string dsoFile, std::string plotOutDir, int aDB, int aDCM, int aFEB)
Definition: DSOPlotMaker.C:98
int num
Definition: f2_nu.C:119
exit(0)
const int File
TFile * file
Definition: cellShifts.C:17
int GetDetID(std::string detector)
Definition: DSOPlotMaker.C:23
enum BeamMode string