Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ExtractTimeOffsetsAndCEff.C
Go to the documentation of this file.
1 // Script to extract time-walk constants for the BCAL
2 
3 #include "TString.h"
4 #include "TFile.h"
5 #include "TH1.h"
6 #include "TH2.h"
7 #include "TSystem.h"
8 #include "TF1.h"
9 #include "TProfile.h"
10 #include "TFitResult.h"
11 #include "TCanvas.h"
12 #include "TPad.h"
13 #include "TStyle.h"
14 #include "TString.h"
15 // #include "T.h"
16 
17 #include <iostream>
18 #include <fstream>
19 #include <sstream>
20 
21 namespace ExtractTimeOffsetsAndCEffNS {
22  //Leave this global so the accesors don't need the pointer as an argument
23  TFile *thisFile;
24 
25  // Accessor functions to grab histograms from our file
26  // (Makes things easy with the HistogramTools fills)
27  TH1I * Get1DHistogram(const char * plugin, const char * directoryName, const char * name, bool print = true){
28  TH1I * histogram;
29  TString fullName = TString(plugin) + "/" + TString(directoryName) + "/" + TString(name);
30  thisFile->GetObject(fullName, histogram);
31  if (histogram == 0){
32  if (print) cout << "Unable to find histogram " << fullName.Data() << endl;
33  return NULL;
34  }
35  return histogram;
36  }
37 
38  TH1I * Get1DHistogramBasedir(const char * name, bool print = true){
39  TH1I * histogram;
40  TString fullName = TString(name);
41  thisFile->GetObject(fullName, histogram);
42  if (histogram == 0){
43  if (print) cout << "Unable to find histogram " << fullName.Data() << endl;
44  return NULL;
45  }
46  return histogram;
47  }
48 
49  TH2I * Get2DHistogram(const char * plugin, const char * directoryName, const char * name){
50  TH2I * histogram;
51  TString fullName = TString(plugin) + "/" + TString(directoryName) + "/" + TString(name);
52  thisFile->GetObject(fullName, histogram);
53  if (histogram == 0){
54  cout << "Unable to find histogram " << fullName.Data() << endl;
55  return NULL;
56  }
57  return histogram;
58  }
59 }
60 
61 TH1D* GetGausFitMean(TH2I *hist2d, TString tag, bool print);
62 
63 
64 
65 // Do the extraction of the actual constants
66 void ExtractTimeOffsetsAndCEff(int run = 2931, TString filename = "hd_root.root", TString tag="default"){
67  printf("ExtractTimeOffsetsAndCEff\n");
68 
69  // Open our input and output file
71  char outfilename[255];
72  sprintf(outfilename,"BCALTimeOffsetAndCEff_Results_%i_%s.root",run,tag.Data());
73  TFile *outputFile = TFile::Open(outfilename, "RECREATE");
74  outputFile->mkdir("Fits");
75  outputFile->mkdir("ResultOverview");
76 
77  // Check to make sure it is open
79  cout << "Unable to open file " << filename.Data() << "...Exiting" << endl;
80  return;
81  }
82 
83  // We need to have the value of C_effective used for the position determination to get the time offsets
84  double C_eff = 16.75; // cm/ns
85 
86  // Since we are updating existing constants we will need their current values. We can pipe them in from the CCDB...
87  // We need a place to store them...
88  // double tdiff_u_d[768];
89  //double channel_global_offset[768];
90 
91  //Pipe the current constants into this macro
92  //NOTE: This dumps the "LATEST" values. If you need something else, modify this script.
93  // char command[100];
94  // sprintf(command, "ccdb dump /BCAL/tdiff_u_d:%i:default", run);
95  // FILE* locInputFile = gSystem->OpenPipe(command, "r");
96  // if(locInputFile == NULL)
97  // return;
98  // //get the first (comment) line
99  // char buff[1024];
100  // if(fgets(buff, sizeof(buff), locInputFile) == NULL)
101  // return;
102  // //get the remaining lines
103  // double time;
104  // int counter = 0;
105  // while(fgets(buff, sizeof(buff), locInputFile) != NULL){
106  // std::istringstream locConstantsStream(buff);
107  // locConstantsStream >> time;
108  // cout << "t_up - t_down = " << time << endl;
109  // tdiff_u_d[counter] = time;
110  // counter++;
111  // }
112  // if (counter != 768) cout << "Wrong number of tdiff_u_d entries (" << counter << " != 768)?" << endl;
113  // //Close the pipe
114  // gSystem->ClosePipe(locInputFile);
115 
116  /*
117  sprintf(command, "ccdb dump /BCAL/channel_global_offset:%i:default", run);
118  locInputFile = gSystem->OpenPipe(command, "r");
119  if(locInputFile == NULL)
120  return;
121  //get the first (comment) line
122  //char buff[1024];
123  if(fgets(buff, sizeof(buff), locInputFile) == NULL)
124  return;
125  //get the remaining lines
126  counter = 0;
127  while(fgets(buff, sizeof(buff), locInputFile) != NULL){
128  istringstream locConstantsStream(buff);
129  locConstantsStream >> time;
130  cout << "Channel Global Offset = " << time << endl;
131  channel_global_offset[counter] = time;
132  counter++;
133  }
134  if (counter != 768) cout << "Wrong number of channel global offsets (" << counter << " != 768)?" << endl;
135  //Close the pipe
136  gSystem->ClosePipe(locInputFile);
137  */
138  printf("Loaded CCDB parameters\n");
139 
140  // This stream will be for outputting the results in a format suitable for the CCDB
141  // Will wait to open until needed
142  ofstream channel_global_offset_File;
143  //, tdiff_u_d_File;
144  sprintf(outfilename,"channel_global_offset_BCAL_%i_%s.txt",run,tag.Data());
145  channel_global_offset_File.open(outfilename);
146  //sprintf(outfilename,"tdiff_u_d_BCAL_%i.txt",run);
147  //tdiff_u_d_File.open(outfilename);
148 
149  // Declaration of the fit funtion
150  TF1 *f1 = new TF1("f1", "[0]+[1]*x", -200, 200);
151  f1->SetParLimits(0, -20.0, 20.0);
152  f1->SetParLimits(1, 0.85, 1.15);
153  outputFile->cd("ResultOverview");
154  // Make some histograms to get the distributions of the fit parameters
155  // TH1I *h1_c0 = new TH1I("h1_c0", "Distribution of parameter c_{0}", 100, -15, 15);
156  // TH1I *h1_c1 = new TH1I("h1_c1", "Distribution of parameter c_{1}", 100, 0.85, 1.15);
157  // TH1F *h1_c0_all = new TH1F ("h1_c0_all", "Value of c0; CCDB Index; c0 [cm]", 768, 0.5, 768.5);
158  // TH1F *h1_c1_all = new TH1F ("h1_c1_all", "Value of c1; CCDB Index; c1 [cm]", 768, 0.5, 768.5);
159  // TH2I *h2_c0_c1 = new TH2I("h2_c0_c1", "c_{1} Vs. c_{0}; c_{0}; c_{1}", 100, -15, 15, 100, 0.85, 1.15);
160 
161  TH1I *OldTimeOffsets = new TH1I("OldTimeOffsets", "Total time offset (pre)", 2000, -4, 4);
162  TH1F *OldTimeOffsetsVsChannel = new TH1F("OldTimeOffsetsVsChannel", "Total time offset (pre) vs channel;CCDB channel", 768, 0.5, 768.5);
163  TH1I *NewTimeOffsets = new TH1I("NewTimeOffsets", "Total time offset (post)", 2000, -4, 4);
164  TH1F *NewTimeOffsetsVsChannel = new TH1F("NewTimeOffsetsVsChannel", "Total time offset (post) vs channel;CCDB channel", 768, 0.5, 768.5);
165  TH1I *NewTimeOffsetsFit = new TH1I("NewTimeOffsetsFit", "Total time offset fit (post)", 2000, -4, 4);
166  TH1F *NewTimeOffsetsFitVsChannel = new TH1F("NewTimeOffsetsFitVsChannel", "Total time offset fit (post) vs channel;CCDB channel", 768, 0.5, 768.5);
167  TH1I *RelativeTimeOffsetsFit = new TH1I("RelativeTimeOffsetsFit", "Change in time offset", 2000, -4, 4);
168  TH1F *RelativeTimeOffsetsFitVsChannel = new TH1F("RelativeTimeOffsetsFitVsChannel", "Change in time offset vs channel;CCDB channel", 768, 0.5, 768.5);
169 
170  TH1D *priorOffsetHist = NULL;
171  priorOffsetHist = (TH1D*)ExtractTimeOffsetsAndCEffNS::thisFile->Get("BCAL_Global_Offsets/Target Time/CCDB_raw_channel_global_offset");
172  //auto priorOffsetHist = ExtractTimeOffsetsAndCEffNS::Get1DHistogram("BCAL_Global_Offsets", "Target Time", "CCDB_raw_channel_global_offset");
173 
174  //auto residualHist = ExtractTimeOffsetsAndCEffNS::Get2DHistogram("BCAL_Global_Offsets", "Target Time", "deltaTVsCell_q-");
175  auto residualHist = ExtractTimeOffsetsAndCEffNS::Get2DHistogram("BCAL_Global_Offsets", "Target Time", "deltaTVsCell_q0");
176 
177  if (priorOffsetHist != NULL) {
178  for (int i = 1 ; i <= priorOffsetHist->GetNbinsX(); i++){
179  printf("%4i %-.3f\n", i,priorOffsetHist->GetBinContent(i));
180  }
181  } else {
182  printf("Failed to load %s\n","CCDB_raw_channel_global_offset");
183  }
184 
185 
186  if(residualHist != NULL && priorOffsetHist != NULL){
187  int nBinsX = residualHist->GetNbinsX();
188  int nBinsY = residualHist->GetNbinsY();
189 
190  TString newtag = tag + "_channel";
191  bool print=1;
192  TH1D *meanhist = GetGausFitMean(residualHist,tag.Data(),print);
193 
194  for (int i = 1 ; i <= nBinsX; i++){
195  TH1D *projY = residualHist->ProjectionY("temp", i, i);
196  // Scan over the histogram
197  float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
198  float timeWindow = 0.5; //ns (Full Width)
199  int binWindow = int(timeWindow / nsPerBin);
200  double maxEntries = 0;
201  double maxMean = 0;
202  for (int j = 1 ; j <= projY->GetNbinsX();j++){
203  int minBin = j;
204  int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
205  double sum = 0, nEntries = 0;
206  for (int bin = minBin; bin <= maxBin; bin++){
207  sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
208  nEntries += projY->GetBinContent(bin);
209  if (bin == maxBin){
210  if (nEntries > maxEntries) {
211  maxMean = sum / nEntries;
212  maxEntries = nEntries;
213  }
214  }
215  }
216  }
217  float priorOffset = priorOffsetHist->GetBinContent(i);
218  float residualOffset = maxMean;
219  float residualOffsetFit = meanhist->GetBinContent(i);
220  float newOffset = priorOffset + residualOffset;
221  float newOffsetFit = priorOffset + residualOffsetFit;
222 
223  OldTimeOffsets->Fill(priorOffset);
224  OldTimeOffsetsVsChannel->SetBinContent(i,priorOffset);
225  NewTimeOffsets->Fill(newOffset);
226  NewTimeOffsetsVsChannel->SetBinContent(i,newOffset);
227  NewTimeOffsetsFit->Fill(newOffsetFit);
228  NewTimeOffsetsFitVsChannel->SetBinContent(i,newOffsetFit);
229  RelativeTimeOffsetsFit->Fill(residualOffsetFit);
230  RelativeTimeOffsetsFitVsChannel->SetBinContent(i,residualOffsetFit);
231 
232  printf("%4i %6.3f %+.3f = %6.3f %+.3f = %6.3f\n",
233  i,priorOffset,residualOffset,newOffset,residualOffsetFit,newOffsetFit);
234  channel_global_offset_File << newOffsetFit << endl;;
235 
236  //selectedBCALOffset->SetBinContent(i, maxMean);
237  //BCALOffsetDistribution->Fill(maxMean);
238  }
239  }
240 
241  // outputFile->cd("Fits");
242  // // Now we want to loop through all available module/layer/sector and try to make a fit of each one
243  // for (unsigned int iModule = 1; iModule <=48; iModule++){
244  // for (unsigned int iLayer = 1; iLayer <= 4; iLayer++){ // Only 3 layers with TDCs
245  // for (unsigned int iSector = 1; iSector <= 4; iSector++){
246  // int the_cell = (iModule - 1) * 16 + (iLayer - 1) * 4 + iSector;
247  // int the_tdc_cell = (iModule - 1) * 12 + (iLayer - 1) * 4 + iSector; // One less layer of TDCs
248  // // Format the string to lookup the histogram by name
249  // char name[200];
250  // sprintf(name, "Module%.2iLayer%.2iSector%.2i", iModule, iLayer, iSector);
251 
252  // // These histograms are created on the fly in the plugin, so there is a chance that they do not exist, in which case the pointer will be NULL
253 
254  // TH2I *h_offsets = ExtractTimeOffsetsAndCEffNS::Get2DHistogram ("BCAL_TDC_Offsets", "Z Position", name);
255 
256  // // Use FitSlicesY routine to extract the mean of each x bin
257  // TObjArray ySlices;
258 
259  // if (h_offsets != NULL) {
260  // h_offsets->RebinX(5);
261  // TProfile *profile = h_offsets->ProfileX();
262  // f1->SetParameters(0, 1); // Just out initial guess
263  // TFitResultPtr fr = profile->Fit(f1, "SQR");
264  // Int_t fitStatus = fr;
265  // if (fitStatus == 0){
266  // double c0 = fr->Parameter(0);
267  // double c0_err = fr->ParError(0);
268  // double c1 = fr->Parameter(1);
269  // double c1_err = fr->ParError(1);
270  // if (c0 == 10.0 || c0 == -10.0 || c1 == 0.9 || c1 == 1.1){
271  // cout << "WARNING: Parameter hit limit " << name << endl;
272  // }
273  // h1_c0->Fill(c0); h1_c1->Fill(c1); h2_c0_c1->Fill(c0,c1);
274  // h1_c0_all->SetBinContent(the_cell, c0); h1_c0_all->SetBinError(the_cell, c0_err);
275  // h1_c1_all->SetBinContent(the_cell, c1); h1_c1_all->SetBinError(the_cell, c1_err);
276  // tdiff_u_d_File << tdiff_u_d[the_cell - 1] + c0 / C_eff << endl;
277  // }
278  // else {
279  // cout << "WARNING: Fit Status "<< fitStatus << " for Upstream " << name << endl;
280  // tdiff_u_d_File << tdiff_u_d[the_cell - 1] << endl;
281  // }
282  // }
283  // else{
284  // tdiff_u_d_File << tdiff_u_d[the_cell - 1] << endl;
285  // }
286  // }
287  // }
288  // }
289  channel_global_offset_File.close();
290  // tdiff_u_d_File.close();
291  outputFile->Write();
293 }
294 
295 
296 TH1D* GetGausFitMean(TH2I *hist2d, TString tag, bool print=0) {
297  // Take a 2D histogram and return the Gaussian mean of each X bin (somewhat like a profile)
298  double outermintime=-2, outermaxtime=2;
299 
300  TAxis * xaxis = hist2d->GetXaxis();
301  TAxis * yaxis = hist2d->GetYaxis();
302  int nBinsX = hist2d->GetNbinsX();
303  int nBinsY = hist2d->GetNbinsY();
304 
305  printf("bins (%i,%i)\n",nBinsX,nBinsY);
306 
307  TH1D *xmean = hist2d->ProjectionX("xmean");
308  xmean->Reset();
309 
310  gStyle->SetOptStat(2210);
311  gStyle->SetOptFit(1);
312 
313  TCanvas *canvas_fit;
314  canvas_fit = new TCanvas("canvas_fit","canvas_fit");
315  gPad->SetLogy();
316 
317  for (int i = 1 ; i <= nBinsX; i++){
318  TH1D *singlebin_projY = hist2d->ProjectionY("temp", i, i);
319  float center = yaxis->GetBinCenter(singlebin_projY->GetMaximumBin());
320  float width=1;
321  float mintime = center-width;
322  float maxtime = center+width;
323 
324  TF1 *f_gaus = new TF1("f_gaus","gaus",mintime,maxtime);
325  //f_gaus->SetLineColor(kRed);
326  //f_gaus->SetLineWidth(1);
327  singlebin_projY->SetAxisRange(mintime,maxtime);
328  double mean=0, meanerr=0;
329 
330  singlebin_projY->Fit(f_gaus,"RQ");
331  mean = f_gaus->GetParameter(1);
332  meanerr = f_gaus->GetParError(1);
333  xmean->SetBinContent(i,mean);
334  xmean->SetBinError(i,meanerr);
335  if (print || meanerr>0.5) {
336  char plotname[255];
337  sprintf(plotname,"plots/gausfitformean/fit_%s_%i.png",tag.Data(),i);
338  canvas_fit->Print(plotname);
339  }
340  singlebin_projY->Reset();
341  }
342  return xmean;
343 }
TH1F * hist2d[Idx+1]
Definition: readhist.C:16
sprintf(text,"Post KinFit Cut")
TString filename
TH1I * Get1DHistogramBasedir(const char *name, bool print=true)
TFile f1(fnam)
TH2I * Get2DHistogram(const char *plugin, const char *directoryName, const char *name)
TH1I * Get1DHistogram(const char *plugin, const char *directoryName, const char *name, bool print=true)
void ExtractTimeOffsetsAndCEff(int run=2931, TString filename="hd_root.root", TString tag="default")
printf("string=%s", string)
TH1D * GetGausFitMean(TH2I *hist2d, TString tag, bool makeplots)