Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Calibration/HLDetectorTiming/FitScripts/ExtractTrackBasedTiming.C
Go to the documentation of this file.
1 namespace ExtractTrackBasedTimingNS {
2  TFile * thisFile;
3 
4  TH1I * Get1DHistogram(const char * plugin, const char * directoryName, const char * name){
5  TH1I * histogram;
6  TString fullName = TString(plugin) + "/" + TString(directoryName) + "/" + TString(name);
7  thisFile->GetObject(fullName, histogram);
8  if (histogram == 0){
9  cout << "Unable to find histogram " << fullName.Data() << endl;
10  return NULL;
11  }
12  return histogram;
13  }
14 
15  TH2I * Get2DHistogram(const char * plugin, const char * directoryName, const char * name){
16  TH2I * histogram;
17  TString fullName = TString(plugin) + "/" + TString(directoryName) + "/" + TString(name);
18  thisFile->GetObject(fullName, histogram);
19  if (histogram == 0){
20  cout << "Unable to find histogram " << fullName.Data() << endl;
21  return NULL;
22  }
23  return histogram;
24  }
25 };
26 
27 void GetCCDBConstants(TString path, Int_t run, TString variation, vector<double>& container, Int_t column = 1){
28  char command[256];
29 
30  sprintf(command, "ccdb dump %s:%i:%s", path.Data(), run, variation.Data());
31  FILE* inputPipe = gSystem->OpenPipe(command, "r");
32  if(inputPipe == NULL)
33  return;
34  //get the first (comment) line
35  char buff[1024];
36  if(fgets(buff, sizeof(buff), inputPipe) == NULL)
37  return;
38  //get the remaining lines
39  double entry;
40  int counter = 0;
41  while(fgets(buff, sizeof(buff), inputPipe) != NULL){
42  istringstream locConstantsStream(buff);
43  while (locConstantsStream >> entry){
44  counter++;
45  if (counter % column == 0) container.push_back(entry);
46  }
47  }
48  //Close the pipe
49  gSystem->ClosePipe(inputPipe);
50 }
51 
52 //Overload this function to handle the base time offsets
53 void GetCCDBConstants1(TString path, Int_t run, TString variation, double& constant1){
54  char command[256];
55  sprintf(command, "ccdb dump %s:%i:%s", path.Data(), run, variation.Data());
56  FILE* inputPipe = gSystem->OpenPipe(command, "r");
57  if(inputPipe == NULL)
58  return;
59  //get the first (comment) line
60  char buff[1024];
61  if(fgets(buff, sizeof(buff), inputPipe) == NULL)
62  return;
63  //get the line containing the values
64  while(fgets(buff, sizeof(buff), inputPipe) != NULL){
65  istringstream locConstantsStream(buff);
66  locConstantsStream >> constant1;
67  }
68  //Close the pipe
69  gSystem->ClosePipe(inputPipe);
70 }
71 
72 void GetCCDBConstants2(TString path, Int_t run, TString variation, double& constant1, double& constant2){
73  char command[256];
74  sprintf(command, "ccdb dump %s:%i:%s", path.Data(), run, variation.Data());
75  FILE* inputPipe = gSystem->OpenPipe(command, "r");
76  if(inputPipe == NULL)
77  return;
78  //get the first (comment) line
79  char buff[1024];
80  if(fgets(buff, sizeof(buff), inputPipe) == NULL)
81  return;
82  //get the line containing the values
83  while(fgets(buff, sizeof(buff), inputPipe) != NULL){
84  istringstream locConstantsStream(buff);
85  locConstantsStream >> constant1 >> constant2;
86  }
87  //Close the pipe
88  gSystem->ClosePipe(inputPipe);
89 }
90 
91 int GetCCDBIndexTAGM(unsigned int column, unsigned int row){
92  int CCDBIndex = column + row;
93  if (column > 9) CCDBIndex += 5;
94  if (column > 27) CCDBIndex += 5;
95  if (column > 81) CCDBIndex += 5;
96  if (column > 99) CCDBIndex += 5;
97 
98  return CCDBIndex;
99 }
100 
101 int GetF1TDCslotTAGH(int id) {
102  double N = 32.0; // channels per slot
103  if (id >= 132 && id <= 172) throw("TAGH: unknown id in [132,172]");
104  int HVid = (id <= 131) ? id : (id - 274 + 233);
105  return int((HVid-1)/N) + 1;
106 }
107 
108 void ExtractTrackBasedTiming(TString fileName = "hd_root.root", int runNumber = 10390, TString variation = "default", bool verbose = false,TString prefix = ""){
109 
110  // set "prefix" in case you want to ship the txt files elsewhere...
111  cout << "Performing Track Matched timing fits for File: " << fileName.Data() << " Run: " << runNumber << " Variation: " << variation.Data() << endl;
112 
113  ExtractTrackBasedTimingNS::thisFile = TFile::Open( fileName , "UPDATE");
115  cout << "Unable to open file " << fileName.Data() << "...Exiting" << endl;
116  return;
117  }
118 
119  //We need the existing constants, The best we can do here is just read them from the file.
120  vector<double> sc_tdc_time_offsets;
121  vector<double> sc_fadc_time_offsets;
122  vector<double> tof_tdc_time_offsets;
123  vector<double> tof_fadc_time_offsets;
124  vector<double> tagm_tdc_time_offsets;
125  vector<double> tagm_fadc_time_offsets;
126  vector<double> tagh_tdc_time_offsets;
127  vector<double> tagh_fadc_time_offsets;
128  vector<double> tagh_counter_quality;
129 
130  double sc_t_base_fadc, sc_t_base_tdc;
131  double tof_t_base_fadc, tof_t_base_tdc;
132  double bcal_t_base_fadc, bcal_t_base_tdc;
133  double tagm_t_base_fadc, tagm_t_base_tdc;
134  double tagh_t_base_fadc, tagh_t_base_tdc;
135  double fdc_t_base_fadc, fdc_t_base_tdc;
136  double fcal_t_base;
137  double cdc_t_base;
138  double RF_Period;
139 
140  cout << "Grabbing CCDB constants..." << endl;
141  // Base times
142  GetCCDBConstants1("/CDC/base_time_offset" ,runNumber, variation, cdc_t_base);
143  GetCCDBConstants1("/FCAL/base_time_offset",runNumber, variation, fcal_t_base);
144  GetCCDBConstants1("/PHOTON_BEAM/RF/beam_period",runNumber, variation, RF_Period);
145  GetCCDBConstants2("/FDC/base_time_offset" ,runNumber, variation, fdc_t_base_fadc, fdc_t_base_tdc);
146  GetCCDBConstants2("/BCAL/base_time_offset" ,runNumber, variation, bcal_t_base_fadc, bcal_t_base_tdc);
147  GetCCDBConstants2("/PHOTON_BEAM/microscope/base_time_offset" ,runNumber, variation, tagm_t_base_fadc, tagm_t_base_tdc);
148  GetCCDBConstants2("/PHOTON_BEAM/hodoscope/base_time_offset" ,runNumber, variation, tagh_t_base_fadc, tagh_t_base_tdc);
149  GetCCDBConstants2("/START_COUNTER/base_time_offset" ,runNumber, variation, sc_t_base_fadc, sc_t_base_tdc);
150  GetCCDBConstants2("/TOF/base_time_offset" ,runNumber, variation, tof_t_base_fadc, tof_t_base_tdc);
151  // Per channel
152  //GetCCDBConstants("/BCAL/TDC_offsets" ,runNumber, variation, bcal_tdc_offsets);
153  //GetCCDBConstants("/FCAL/timing_offsets" ,runNumber, variation, fcal_adc_offsets);
154  GetCCDBConstants("/START_COUNTER/adc_timing_offsets" ,runNumber, variation, sc_fadc_time_offsets);
155  GetCCDBConstants("/START_COUNTER/tdc_timing_offsets" ,runNumber, variation, sc_tdc_time_offsets);
156  GetCCDBConstants("/PHOTON_BEAM/microscope/fadc_time_offsets" ,runNumber, variation, tagm_fadc_time_offsets,3);// Interested in 3rd column
157  GetCCDBConstants("/PHOTON_BEAM/microscope/tdc_time_offsets" ,runNumber, variation, tagm_tdc_time_offsets,3);
158  GetCCDBConstants("/PHOTON_BEAM/hodoscope/fadc_time_offsets" ,runNumber, variation, tagh_fadc_time_offsets,2);// Interested in 2nd column
159  GetCCDBConstants("/PHOTON_BEAM/hodoscope/tdc_time_offsets" ,runNumber, variation, tagh_tdc_time_offsets,2);
160  GetCCDBConstants("/PHOTON_BEAM/hodoscope/counter_quality" ,runNumber, variation, tagh_counter_quality,2);
161  GetCCDBConstants("/TOF/adc_timing_offsets",runNumber, variation, tof_fadc_time_offsets);
162  GetCCDBConstants("/TOF/timing_offsets",runNumber, variation, tof_tdc_time_offsets);
163 
164  cout << "CDC base times = " << cdc_t_base << endl;
165  cout << "FCAL base times = " << fcal_t_base << endl;
166  cout << "FDC base times = " << fdc_t_base_fadc << ", " << fdc_t_base_tdc << endl;
167  cout << "BCAL base times = " << bcal_t_base_fadc << ", " << bcal_t_base_tdc << endl;
168  cout << "SC base times = " << sc_t_base_fadc << ", " << sc_t_base_tdc << endl;
169  cout << "TOF base times = " << tof_t_base_fadc << ", " << tof_t_base_tdc << endl;
170  cout << "TAGH base times = " << tagh_t_base_fadc << ", " << tagh_t_base_tdc << endl;
171  cout << "TAGM base times = " << tagm_t_base_fadc << ", " << tagm_t_base_tdc << endl;
172 
173  cout << endl;
174  cout << "RF_Period = " << RF_Period << endl;
175  cout << endl;
176 
177  cout << "Done grabbing CCDB constants...Entering fits..." << endl;
178 
179  // Do our final step in the timing alignment with tracking
180 
181  //When the RF is present we can try to simply pick out the correct beam bucket for each of the runs
182  //First just a simple check to see if we have the appropriate data
183  bool useRF = false;
184  TH1I *testHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "TAGH_TDC_RF_Compare","Counter ID 001");
185  if (testHist != NULL){ // Not great since we rely on channel 1 working, but can be craftier later.
186  cout << "Using RF Times for Calibration" << endl;
187  useRF = true;
188  }
189  ofstream outFile;
190  TH2I *thisHist;
191  thisHist = ExtractTrackBasedTimingNS::Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - SC Target Time");
192  if (useRF) thisHist = ExtractTrackBasedTimingNS::Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGM - RFBunch Time");
193  if (thisHist != NULL){
194  //Statistics on these histograms are really quite low we will have to rebin and do some interpolation
195  outFile.open(prefix + "tagm_tdc_timing_offsets.txt", ios::out | ios::trunc);
196  outFile.close(); // clear file
197  outFile.open(prefix + "tagm_adc_timing_offsets.txt", ios::out | ios::trunc);
198  outFile.close(); // clear file
199  int nBinsX = thisHist->GetNbinsX();
200  int nBinsY = thisHist->GetNbinsY();
201  TH1D * selectedTAGMOffset = new TH1D("selectedTAGMOffset", "Selected TAGM Offset; Column; Offset [ns]", nBinsX, 0.5, nBinsX + 0.5);
202  TH1I * TAGMOffsetDistribution = new TH1I("TAGMOffsetDistribution", "TAGM Offset; TAGM Offset [ns]; Entries", 500, -250, 250);
203  for (int i = 1 ; i <= nBinsX; i++){
204  TH1D *projY = thisHist->ProjectionY("temp", i, i);
205  // Scan over the histogram
206  //chose the correct number of bins based on the histogram
207  float nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
208  float timeWindow = 3; //ns (Full Width)
209  int binWindow = int(timeWindow / nsPerBin);
210  double maxEntries = 0;
211  double maxMean = 0;
212  for (int j = 1 ; j <= projY->GetNbinsX();j++){
213  int minBin = j;
214  int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
215  double sum = 0, nEntries = 0;
216  for (int bin = minBin; bin <= maxBin; bin++){
217  sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
218  nEntries += projY->GetBinContent(bin);
219  if (bin == maxBin){
220  if (nEntries > maxEntries) {
221  maxMean = sum / nEntries;
222  maxEntries = nEntries;
223  }
224  }
225  }
226  }
227  //In the case there is RF, our job is to pick just the number of the correct beam bunch, so that's really all we need.
228  if(useRF) {
229  int beamBucket = int((maxMean / RF_Period) + 0.5); // +0.5 to handle rounding correctly
230  selectedTAGMOffset->SetBinContent(i, beamBucket);
231  TAGMOffsetDistribution->Fill(beamBucket);
232  }
233  else{
234  selectedTAGMOffset->SetBinContent(i, maxMean);
235  TAGMOffsetDistribution->Fill(maxMean);
236  }
237  }
238  double meanOffset = TAGMOffsetDistribution->GetMean();
239  // This might be in units of beam bunches, so we need to convert
240  if (useRF) meanOffset *= RF_Period;
241  if (verbose) {
242  cout << "Dumping TAGM results...\n=======================================" << endl;
243  cout << "TAGM mean Offset = " << meanOffset << endl;
244  cout << "fADC Offsets" << endl;
245  }
246 
247  outFile.open(prefix + "tagm_adc_timing_offsets.txt", ios::out);
248  //for (int i = 1 ; i <= nBinsX; i++){
249  // Loop over rows
250  if (verbose) cout << "Column\tRow\tvalueToUse\toldValue\tmeanOffset\tTotal" << endl;
251  for (unsigned int column = 1; column <= 102; column++){
252  int index = GetCCDBIndexTAGM(column, 0);
253  double valueToUse = selectedTAGMOffset->GetBinContent(index);
254  if (useRF) valueToUse *= RF_Period;
255 
256  //if (valueToUse == 0) valueToUse = meanOffset;
257  outFile << "0 " << column << " " << valueToUse + tagm_fadc_time_offsets[index-1] - meanOffset<< endl;
258  if (verbose) printf("0\t%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", column, valueToUse, tagm_fadc_time_offsets[index-1], meanOffset,
259  valueToUse + tagm_fadc_time_offsets[index-1] - meanOffset);
260  if (column == 9 || column == 27 || column == 81 || column == 99){
261  for (unsigned int row = 1; row <= 5; row++){
262  index = GetCCDBIndexTAGM(column, row);
263  valueToUse = selectedTAGMOffset->GetBinContent(index);
264  if (useRF) valueToUse *= RF_Period;
265  //if (valueToUse == 0) valueToUse = meanOffset;
266  outFile << row << " " << column << " " << valueToUse + tagm_fadc_time_offsets[index-1] - meanOffset<< endl;
267  if (verbose) printf("%i\t%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", row, column, valueToUse, tagm_fadc_time_offsets[index-1], meanOffset,
268  valueToUse + tagm_fadc_time_offsets[index-1] - meanOffset);
269  }
270  }
271  }
272  outFile.close();
273 
274  if (verbose) {
275  cout << "TDC Offsets" << endl;
276  cout << "Column\tRow\tvalueToUse\toldValue\tmeanOffset\tTotal" << endl;
277  }
278  outFile.open(prefix + "tagm_tdc_timing_offsets.txt", ios::out);
279  //for (int i = 1 ; i <= nBinsX; i++){
280  // Loop over rows
281  for (unsigned int column = 1; column <= 102; column++){
282  int index = GetCCDBIndexTAGM(column, 0);
283  double valueToUse = selectedTAGMOffset->GetBinContent(index);
284  if (useRF) valueToUse *= RF_Period;
285  //if (valueToUse == 0) valueToUse = meanOffset;
286  outFile << "0 " << column << " " << valueToUse + tagm_tdc_time_offsets[index-1] - meanOffset << endl;
287  if (verbose) printf("0\t%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", column, valueToUse, tagm_tdc_time_offsets[index-1], meanOffset,
288  valueToUse + tagm_tdc_time_offsets[index-1] - meanOffset);
289  if (column == 9 || column == 27 || column == 81 || column == 99){
290  for (unsigned int row = 1; row <= 5; row++){
291  index = GetCCDBIndexTAGM(column, row);
292  valueToUse = selectedTAGMOffset->GetBinContent(index);
293  if (useRF) valueToUse *= RF_Period;
294  //if (valueToUse == 0) valueToUse = meanOffset;
295  outFile << row << " " << column << " " << valueToUse + tagm_tdc_time_offsets[index-1] - meanOffset << endl;
296  if (verbose) printf("%i\t%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", row, column, valueToUse, tagm_tdc_time_offsets[index-1], meanOffset,
297  valueToUse + tagm_tdc_time_offsets[index-1] - meanOffset);
298  }
299  }
300  }
301  outFile.close();
302  outFile.open(prefix + "tagm_base_time.txt", ios::out);
303  if (verbose) {
304  printf("TAGM ADC Base = %f - (%f) = %f\n", tagm_t_base_fadc, meanOffset, tagm_t_base_fadc - meanOffset);
305  printf("TAGM TDC Base = %f - (%f) = %f\n", tagm_t_base_tdc, meanOffset, tagm_t_base_tdc - meanOffset);
306  }
307  outFile << tagm_t_base_fadc - meanOffset << " " << tagm_t_base_tdc - meanOffset << endl;
308  outFile.close();
309 
310  }
311 
312  thisHist = ExtractTrackBasedTimingNS::Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGH - SC Target Time");
313  if (useRF) thisHist = ExtractTrackBasedTimingNS::Get2DHistogram("HLDetectorTiming", "TRACKING", "TAGH - RFBunch Time");
314  if (thisHist != NULL) {
315  outFile.open(prefix + "tagh_tdc_timing_offsets.txt", ios::out | ios::trunc);
316  outFile.close(); // clear file
317  outFile.open(prefix + "tagh_adc_timing_offsets.txt", ios::out | ios::trunc);
318  outFile.close(); // clear file
319 
320  // Setup histogram for determining the most probable change in offset for each F1TDC slot
321  // This is needed to account for the occasional uniform shift in offsets of the 32 counters in a slot
322  const int NtdcSlots = 8;
323  TH1I * tdcDist[NtdcSlots];
324  for (int i = 1; i <= NtdcSlots; i++) {
325  stringstream ss; ss << i;
326  TString s = ss.str();
327  double range = 500.0; double width = 0.1;
328  int Nbins = range/width;
329  double low = -0.5*range - 0.5*width;
330  double high = 0.5*range - 0.5*width;
331  tdcDist[i-1] = new TH1I("TAGHOffsetDistribution_"+s, "TAGH Offset (slot "+s+"); TAGH Offset [ns]; Entries", Nbins, low, high);
332  }
333 
334  int nBinsX = thisHist->GetNbinsX();
335  TH1D * selectedTAGHOffset = new TH1D("selectedTAGHOffset", "Selected TAGH Offset; ID; Offset [ns]", nBinsX, 0.5, nBinsX + 0.5);
336  for (int i = 1 ; i <= nBinsX; i++) {
337  TH1D *projY = thisHist->ProjectionY("temp", i, i);
338  // Scan over histogram to find mean offset in timeWindow with largest integral
339  // Choose the correct number of bins based on the histogram
340  double nsPerBin = (projY->GetBinCenter(projY->GetNbinsX()) - projY->GetBinCenter(1)) / projY->GetNbinsX();
341  double timeWindow = 2.0; // ns (Full Width)
342  int binWindow = int(timeWindow / nsPerBin);
343 
344  double maxEntries = 0;
345  double maxMean = 0;
346  for (int j = 1; j <= projY->GetNbinsX(); j++) {
347  int minBin = j;
348  int maxBin = (j + binWindow) <= projY->GetNbinsX() ? (j + binWindow) : projY->GetNbinsX();
349  double sum = 0;
350  double nEntries = 0;
351  for (int bin = minBin; bin <= maxBin; bin++) {
352  sum += projY->GetBinContent(bin) * projY->GetBinCenter(bin);
353  nEntries += projY->GetBinContent(bin);
354  if (bin == maxBin) {
355  if (nEntries > maxEntries) {
356  maxMean = sum / nEntries;
357  maxEntries = nEntries;
358  }
359  }
360  }
361  }
362 
363  if (tagh_counter_quality[i-1] == 0.0) {
364  selectedTAGHOffset->SetBinContent(i, 0);
365  continue;
366  }
367  int tdc_slot = GetF1TDCslotTAGH(i);
368  if (useRF) {
369  int beamBucket;
370  if (maxMean >= 0) beamBucket = int((maxMean / RF_Period) + 0.5); // +0.5 to handle rounding correctly
371  else beamBucket = int((maxMean / RF_Period) - 0.5);
372  selectedTAGHOffset->SetBinContent(i, beamBucket);
373  if (maxEntries != 0.0) tdcDist[tdc_slot - 1]->Fill(beamBucket);
374  } else {
375  selectedTAGHOffset->SetBinContent(i, maxMean);
376  if (maxEntries != 0.0) tdcDist[tdc_slot - 1]->Fill(maxMean);
377  }
378  }
379  // Most probable change in offset or beam bucket per F1TDC slot
380  double mpDelta[NtdcSlots];
381  for (int i = 1; i <= NtdcSlots; i++) {
382  int mpBin = tdcDist[i-1]->GetMaximumBin();
383  mpDelta[i-1] = (mpBin > 0) ? tdcDist[i-1]->GetBinCenter(mpBin) : 0.0;
384  if (useRF) mpDelta[i-1] *= RF_Period;
385  if (verbose) {
386  cout << "TAGH most probable Offset = " << i << ", " << mpDelta[i-1] << endl;
387  }
388  }
389 
390  if (verbose) {
391  cout << "Dumping TAGH results...\n=======================================" << endl;
392  cout << "Type\tChannel\tvalueToUse\toldValue\tmpDelta\tTotal" << endl;
393  }
394 
395  double limit = 2.5; // ns
396  double ccdb_sum = 0.0;
397  for (int i = 1; i <= nBinsX; i++) ccdb_sum += tagh_tdc_time_offsets[i-1];
398  double c1_tdcOffset = 0.0;
399  outFile.open(prefix + "tagh_tdc_timing_offsets.txt");
400  for (int i = 1; i <= nBinsX; i++) {
401  if (tagh_counter_quality[i-1] == 0.0) {
402  outFile << i << " " << 0 << endl;
403  continue;
404  }
405  int tdc_slot = GetF1TDCslotTAGH(i);
406  double delta = selectedTAGHOffset->GetBinContent(i);
407  if (useRF) delta *= RF_Period;
408  if (ccdb_sum > 0.0 && fabs(delta - mpDelta[tdc_slot-1]) > limit) {
409  delta = mpDelta[tdc_slot-1];
410  }
411  double ccdb = tagh_tdc_time_offsets[i-1];
412  double offset = ccdb + delta;
413  if (i == 1) c1_tdcOffset = offset;
414  offset -= c1_tdcOffset;
415  outFile << i << " " << offset << endl;
416  if (verbose) printf("TDC\t%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", i, delta, ccdb, mpDelta[tdc_slot-1], offset);
417  }
418  outFile.close();
419 
420  ccdb_sum = 0.0;
421  for (int i = 1; i <= nBinsX; i++) ccdb_sum += tagh_fadc_time_offsets[i-1];
422  double c1_adcOffset = 0.0;
423  outFile.open(prefix + "tagh_adc_timing_offsets.txt");
424  for (int i = 1; i <= nBinsX; i++) {
425  if (tagh_counter_quality[i-1] == 0.0) {
426  outFile << i << " " << 0 << endl;
427  continue;
428  }
429  int tdc_slot = GetF1TDCslotTAGH(i);
430  double delta = selectedTAGHOffset->GetBinContent(i);
431  if (useRF) delta *= RF_Period;
432  if (ccdb_sum > 0.0 && fabs(delta - mpDelta[tdc_slot-1]) > limit) {
433  delta = mpDelta[tdc_slot-1];
434  }
435  double ccdb = tagh_fadc_time_offsets[i-1];
436  double offset = ccdb + delta;
437  if (i == 1) c1_adcOffset = offset;
438  offset -= c1_adcOffset;
439  outFile << i << " " << offset << endl;
440  if (verbose) printf("ADC\t%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n", i, delta, ccdb, mpDelta[tdc_slot-1], offset);
441  }
442  outFile.close();
443 
444  outFile.open(prefix + "tagh_base_time.txt");
445  outFile << tagh_t_base_fadc - c1_adcOffset << " " << tagh_t_base_tdc - c1_tdcOffset << endl;
446  if (verbose) {
447  printf("TAGH ADC Base = %f - (%f) = %f\n", tagh_t_base_fadc, c1_adcOffset, tagh_t_base_fadc - c1_adcOffset);
448  printf("TAGH TDC Base = %f - (%f) = %f\n", tagh_t_base_tdc, c1_tdcOffset, tagh_t_base_tdc - c1_tdcOffset);
449  }
450  outFile.close();
451  }
452 
453  // We can use the RF time to calibrate the SC time (Experimental for now)
454  double meanSCOffset = 0.0; // In case we change the time of the SC, we need this in this scope
455  if(useRF){
456  TH1F * selectedSCSectorOffset = new TH1F("selectedSCSectorOffset", "Selected TDC-RF offset;Sector; Time", 30, 0.5, 30.5);
457  TH1F * selectedSCSectorOffsetDistribution = new TH1F("selectedSCSectorOffsetDistribution", "Selected TDC-RF offset;Time;Entries", 100, -3.0, 3.0);
458  TF1* f = new TF1("f","pol0(0)+gaus(1)", -3.0, 3.0);
459  for (int sector = 1; sector <= 30; sector++){
460  TH1I *scRFHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "SC_Target_RF_Compare", Form("Sector %.2i", sector));
461  if (scRFHist == NULL) continue;
462  //Do the fit
463  TFitResultPtr fr = scRFHist->Fit("pol0", "SQ", "", -2, 2);
464  double p0 = fr->Parameter(0);
465 
466  f->FixParameter(0,p0);
467  f->SetParLimits(2, -2, 2);
468  f->SetParLimits(3, 0, 2);
469  f->SetParameter(1, 10);
470  f->SetParameter(2, scRFHist->GetBinCenter(scRFHist->GetMaximumBin()));
471  f->SetParameter(3, 0);
472 
473  fr = scRFHist->Fit(f, "SQ", "", -2, 2);
474  double SCOffset = fr->Parameter(2);
475  selectedSCSectorOffset->SetBinContent(sector, SCOffset);
476  selectedSCSectorOffsetDistribution->Fill(SCOffset);
477  }
478  // Now write out the offsets
479  meanSCOffset = selectedSCSectorOffsetDistribution->GetMean();
480  if (verbose){
481  cout << "Dumping SC results...\n=======================================" << endl;
482  cout << "SC mean Offset = " << meanSCOffset << endl;
483  cout << "TDC Offsets" << endl;
484  cout << "Sector\toldValue\tValueToUse\tmeanOffset\tTotal" << endl;
485  }
486  outFile.open(prefix + "sc_tdc_timing_offsets.txt");
487  for (int sector = 1; sector <= 30; sector++){
488  outFile << sc_tdc_time_offsets[sector-1] + selectedSCSectorOffset->GetBinContent(sector) - meanSCOffset << endl;
489  if (verbose) printf("%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n",sector, sc_tdc_time_offsets[sector-1], selectedSCSectorOffset->GetBinContent(sector), meanSCOffset,
490  sc_tdc_time_offsets[sector-1] + selectedSCSectorOffset->GetBinContent(sector) - meanSCOffset);
491  }
492  outFile.close();
493  if (verbose){
494  cout << "ADC Offsets" << endl;
495  cout << "Sector\tvalueToUse\toldValue\tmeanOffset\tTotal" << endl;
496  }
497  outFile.open(prefix + "sc_adc_timing_offsets.txt");
498  for (int sector = 1; sector <= 30; sector++){
499  outFile << sc_fadc_time_offsets[sector-1] + selectedSCSectorOffset->GetBinContent(sector) - meanSCOffset << endl;
500  if (verbose) printf("%i\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\n",sector,sc_fadc_time_offsets[sector-1], selectedSCSectorOffset->GetBinContent(sector), meanSCOffset,
501  sc_fadc_time_offsets[sector-1] + selectedSCSectorOffset->GetBinContent(sector) - meanSCOffset);
502  }
503  outFile.close();
504 
505  outFile.open(prefix + "sc_base_time.txt");
506  outFile << sc_t_base_fadc - meanSCOffset << " " << sc_t_base_tdc - meanSCOffset << endl;
507  if (verbose) {
508  printf("SC ADC Base = %f - (%f) = %f\n", sc_t_base_fadc, meanSCOffset, sc_t_base_fadc - meanSCOffset);
509  printf("SC TDC Base = %f - (%f) = %f\n", sc_t_base_tdc, meanSCOffset, sc_t_base_tdc - meanSCOffset);
510  }
511  outFile.close();
512  }
513 
514  TH1I *this1DHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "TRACKING", "TOF - RF Time");
515  if(this1DHist != NULL){
516  //Gaussian
517  Double_t maximum = this1DHist->GetBinCenter(this1DHist->GetMaximumBin());
518  TFitResultPtr fr = this1DHist->Fit("gaus", "S", "", maximum - 1.5, maximum + 1.5);
519  float mean = fr->Parameter(1);
520  outFile.open(prefix + "tof_base_time.txt");
521  if (verbose) {
522  printf("TOF ADC Base = %f - (%f) - (%f) = %f\n", tof_t_base_fadc, mean, meanSCOffset, tof_t_base_fadc - mean - meanSCOffset);
523  printf("TOF TDC Base = %f - (%f) - (%f) = %f\n", tof_t_base_tdc, mean, meanSCOffset, tof_t_base_tdc - mean - meanSCOffset);
524  }
525  outFile << tof_t_base_fadc - mean - meanSCOffset<< " " << tof_t_base_tdc - mean - meanSCOffset<< endl;
526  outFile.close();
527  }
528 
529  this1DHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "TRACKING", "BCAL - RF Time");
530  if(this1DHist != NULL){
531  //Gaussian
532  Double_t maximum = this1DHist->GetBinCenter(this1DHist->GetMaximumBin());
533  TFitResultPtr fr = this1DHist->Fit("gaus", "S", "", maximum - 5, maximum + 5);
534  float mean = fr->Parameter(1);
535  outFile.open(prefix + "bcal_base_time.txt");
536  if (verbose) {
537  printf("BCAL ADC Base = %f - (%f) - (%f) = %f\n", bcal_t_base_fadc, mean, meanSCOffset, bcal_t_base_fadc - mean - meanSCOffset);
538  printf("BCAL TDC Base = %f - (%f) - (%f) = %f\n", bcal_t_base_tdc, mean, meanSCOffset, bcal_t_base_tdc - mean - meanSCOffset);
539  }
540  outFile << bcal_t_base_fadc - mean - meanSCOffset << " " << bcal_t_base_tdc - mean - meanSCOffset << endl; // TDC info not used
541  outFile.close();
542  }
543 
544  this1DHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "TRACKING", "FCAL - RF Time");
545  if(this1DHist != NULL){
546  //Gaussian
547  Double_t maximum = this1DHist->GetBinCenter(this1DHist->GetMaximumBin());
548  TFitResultPtr fr = this1DHist->Fit("gaus", "S", "", maximum - 5, maximum + 5);
549  float mean = fr->Parameter(1);
550  outFile.open(prefix + "fcal_base_time.txt");
551  if (verbose) {
552  printf("FCAL ADC Base = %f - (%f) - (%f) = %f\n",fcal_t_base, mean, meanSCOffset, fcal_t_base - mean - meanSCOffset);
553  }
554  outFile << fcal_t_base - mean - meanSCOffset<< endl;
555  outFile.close();
556  }
557 
558  this1DHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "TRACKING", "Earliest CDC Time Minus Matched SC Time");
559  if(this1DHist != NULL){
560  //Gaussian
561  Double_t maximum = this1DHist->GetBinCenter(this1DHist->GetMaximumBin());
562  TFitResultPtr fr = this1DHist->Fit("gaus", "S", "", maximum - 15, maximum + 10);
563  float mean = fr->Parameter(1);
564  outFile.open(prefix + "cdc_base_time.txt");
565  if (verbose) {
566  printf("CDC ADC Base = %f - (%f) - (%f) = %f\n",cdc_t_base, mean, meanSCOffset, cdc_t_base - mean - meanSCOffset);
567  }
568  outFile << cdc_t_base - mean - meanSCOffset << endl;
569  outFile.close();
570  }
571 
572  // We want to account for any residual difference between the cathode and anode times.
573  double FDC_ADC_Offset = 0.0, FDC_TDC_Offset = 0.0;
574  this1DHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "FDC", "FDCHit Cathode time;1");
575  if(this1DHist != NULL){
576  Int_t firstBin = this1DHist->FindFirstBinAbove( 1 , 1); // Find first bin with content above 1 in the histogram
577  for (int i = 0; i <= 16; i++){
578  if ((firstBin + i) > 0) this1DHist->SetBinContent((firstBin + i), 0);
579  }
580  //Fit a gaussian to the left of the main peak
581  Double_t maximum = this1DHist->GetBinCenter(this1DHist->GetMaximumBin());
582  TF1 *f = new TF1("f", "gaus");
583  f->SetParameters(100, maximum, 20);
584  //this1DHist->Rebin(2);
585  TFitResultPtr fr = this1DHist->Fit(f, "S", "", maximum - 10, maximum + 7); // Cant fix value at end of range
586  double mean = fr->Parameter(1);
587  float sigma = fr->Parameter(2);
588  FDC_ADC_Offset = mean;
589  delete f;
590  }
591 
592  this1DHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "FDC", "FDCHit Wire time;1");
593  if(this1DHist != NULL){
594  Int_t firstBin = this1DHist->FindLastBinAbove( 1 , 1); // Find first bin with content above 1 in the histogram
595  for (int i = 0; i <= 25; i++){
596  if ((firstBin + i) > 0) this1DHist->SetBinContent((firstBin + i), 0);
597  }
598  //Fit a gaussian to the left of the main peak
599  Double_t maximum = this1DHist->GetBinCenter(this1DHist->GetMaximumBin());
600  TF1 *f = new TF1("f", "gaus");
601  f->SetParameters(100, maximum, 20);
602  TFitResultPtr fr = this1DHist->Fit(f, "S", "", maximum - 10, maximum + 5); // Cant fix value at end of range
603  double mean = fr->Parameter(1);
604  float sigma = fr->Parameter(2);
605  FDC_TDC_Offset = mean;
606  delete f;
607  }
608  double FDC_ADC_TDC_Offset = FDC_ADC_Offset - FDC_TDC_Offset;
609 
610  this1DHist = ExtractTrackBasedTimingNS::Get1DHistogram("HLDetectorTiming", "TRACKING", "Earliest Flight-time Corrected FDC Time");
611  if(this1DHist != NULL){
612  //Landau
613  Double_t maximum = this1DHist->GetBinCenter(this1DHist->GetMaximumBin());
614  TFitResultPtr fr = this1DHist->Fit("landau", "S", "", maximum - 2.5, maximum + 4);
615  float MPV = fr->Parameter(1);
616  outFile.open(prefix + "fdc_base_time.txt");
617  if (verbose) {
618  printf("FDC ADC Base = %f - (%f) - (%f) - (%f) = %f\n",fdc_t_base_fadc, MPV, meanSCOffset, FDC_ADC_TDC_Offset, fdc_t_base_fadc - MPV - meanSCOffset - FDC_ADC_TDC_Offset);
619  printf("FDC TDC Base = %f - (%f) - (%f) = %f\n",fdc_t_base_tdc, MPV, meanSCOffset, fdc_t_base_tdc - MPV - meanSCOffset);
620  }
621  outFile << fdc_t_base_fadc - MPV - meanSCOffset - FDC_ADC_TDC_Offset << " " << fdc_t_base_tdc - MPV - meanSCOffset << endl;
622  outFile.close();
623  }
624 
626  return;
627 }
void ExtractTrackBasedTiming(TString fileName="hd_root.root", int runNumber=10390, TString variation="default", bool verbose=false, TString prefix="")
TFile * outFile
Definition: FitGains.C:15
sprintf(text,"Post KinFit Cut")
TH1I * Get1DHistogram(const char *plugin, const char *directoryName, const char *name)
static char index(char c)
Definition: base64.cpp:115
TF1 * f
Definition: FitGains.C:21
double counter
Definition: FitGains.C:151
TH2I * Get2DHistogram(const char *plugin, const char *directoryName, const char *name)
void GetCCDBConstants(TString path, Int_t run, TString variation, vector< double > &container, Int_t column=1)
Double_t sigma[NCHANNELS]
Definition: st_tw_resols.C:37
void GetCCDBConstants2(TString path, Int_t run, TString variation, double &constant1, double &constant2)
void GetCCDBConstants1(TString path, Int_t run, TString variation, double &constant1)
int GetCCDBIndexTAGM(unsigned int column, unsigned int row)
printf("string=%s", string)