source: trunk/Mars/hawc/synchron.C@ 19953

Last change on this file since 19953 was 19951, checked in by giangdo, 5 years ago
script to synchronise the data from HAWC's Eye and HAWC
  • Property svn:executable set to *
File size: 23.9 KB
Line 
1#include <iostream>
2#include <algorithm>
3#include <vector>
4#include <cstdio>
5#include <iterator>
6
7#include "TTree.h"
8#include "TChain.h"
9#include "TGraphErrors.h"
10#include "TGraph.h"
11#include "TCanvas.h"
12#include "TFile.h"
13#include "TF1.h"
14#include "TH1.h"
15#include "TH2.h"
16#include "TH1F.h"
17#include "TLegend.h"
18#include <TStyle.h>
19
20#include "TROOT.h"
21#include "TBrowser.h"
22
23#include "MTime.h"
24//#include "MSignalCam.h"
25#include "MRawEvtHeader.h"
26//#include "MSoftwareTrigger.h"
27#include "MRawBoardsFACT.h"
28
29using namespace std;
30/*
31 C : a character string terminated by the 0 character
32 B : an 8 bit signed integer (Char_t)
33 b : an 8 bit unsigned integer (UChar_t)
34 S : a 16 bit signed integer (Short_t)
35 s : a 16 bit unsigned integer (UShort_t)
36 I : a 32 bit signed integer (Int_t)
37 i : a 32 bit unsigned integer (UInt_t)
38 F : a 32 bit floating point (Float_t)
39 D : a 64 bit floating point (Double_t)
40 L : a 64 bit signed integer (Long64_t)
41 l : a 64 bit unsigned integer (ULong64_t)
42 O : [the letter o, not a zero] a boolean (Bool_t)
43*/
44
45vector<Double_t> calibration(TString HEye_path)
46{
47 /* INPUT: HEye_path: Path to HAWC's Eye data (calibrated or cleaned)
48 OUTPUT: calibrated time in Unix_s
49 */
50 // Vectors to store the Unix time after calibration
51 vector<Double_t> unix_s = {};
52
53 TFile *tHEye = new TFile(HEye_path);
54
55 cout << "------------------------------------------------------------------------" << endl;
56 cout << "Do time calibration with file: " << HEye_path << endl;
57
58 TTree *tEvents = (TTree*) tHEye->Get("Events");
59
60 MTime *mtime_utc = 0; // PC's timestamp in UNIX
61 MRawBoardsFACT *mtime_board = 0; // board's counter
62
63 tEvents->SetBranchAddress("MTime.", &mtime_utc);
64 tEvents->SetBranchAddress("MRawBoardsFACT.", &mtime_board);
65
66 Long64_t nEvents = tEvents->GetEntries();
67 //cout << "Number of Events: " << nEvents << endl;
68
69 Double_t Mjd; // time in Mjd
70 int Mjd_unix = 40587; // Mjd at 01/01/1970
71 Double_t unix_timestamp;
72
73 // For fitting and plotting
74 Double_t Unix_fit[nEvents], Board_fit[nEvents], UnixErr_fit[nEvents];
75
76 for (Long64_t i = 0; i < nEvents; i++)
77 {
78 tEvents->GetEntry(i);
79
80 //mtime_utc->Print();
81 Mjd = mtime_utc->GetMjd();
82 UInt_t board_timestamp = mtime_board->GetFadTime(0);
83
84 // calculate the PC's unix timestamp based on the PC UTC's time
85 Double_t unix_timestamp = (Mjd - Mjd_unix)*24*60*60;
86
87 Unix_fit[i] = unix_timestamp;
88 Board_fit[i] = board_timestamp;
89 UnixErr_fit[i] = 0.001; // jitter of 1ms
90
91 }
92
93 // Use a linear fit to determine the offset and the conversion factor as the fit's slope
94 auto *board_unix = new TGraphErrors(nEvents, Board_fit, Unix_fit, 0, UnixErr_fit);
95 TF1 *fit_timestamp = new TF1("fit_timestamp", "[0] + [1]*x", Board_fit[0], Board_fit[nEvents-1]);
96 board_unix->Fit("fit_timestamp", "S");
97 Double_t offset = fit_timestamp->GetParameter(0);
98 Double_t conv_timestamp = fit_timestamp->GetParameter(1);
99
100 // Calculate the Unix timestamp based on board's counter and the calibration data
101 vector<Double_t> unix_time_s = {};
102 for (Long64_t i = 0; i < nEvents; i++)
103 {
104 tEvents->GetEntry(i);
105 UInt_t board_timestamp = mtime_board->GetFadTime(0);
106 Double_t unix_time = board_timestamp * conv_timestamp + offset;
107 unix_time_s.push_back(unix_time);
108 }
109
110 // Calculate the mean time difference between calibrated unix_timestamp and PC's timestamp
111 Double_t time_diff = 0.0;
112 for (Long64_t i = 0; i < nEvents; i++)
113 {
114 time_diff += TMath::Abs(Unix_fit[i] - unix_time_s[i]);
115 }
116 Double_t mean_time_diff = time_diff/nEvents;
117 printf("Mean time difference between PC's timestamps and calibrated counter's timestamp: %f s\n", mean_time_diff);
118
119/*
120 TCanvas *c = new TCanvas("time_calibration", "Board's time and PC's UTC time",1200,800);
121 c->SetGridx();
122 c->SetGridy();
123
124 board_unix->SetTitle("FAD's timestamp vs. PC's UTC timestamp");
125 board_unix->GetXaxis()->SetTitle("Board's counter");
126 board_unix->GetYaxis()->SetTitle("Unix timestamp [s]");
127 board_unix->SetMarkerStyle(8);
128 board_unix->SetMarkerColor(1);
129 gStyle->SetOptFit(0111);
130
131 board_unix->Draw("AP");
132 */
133
134 return unix_time_s;
135}
136
137// structure of HAWC reconstructed data
138struct rec_struct
139{ Int_t status ;
140 Int_t version ;
141 Int_t eventID ;
142 Int_t runID ;
143 Int_t timeSliceID ;
144 Int_t trigger_flags;
145 Int_t event_flags ;
146 Int_t gtc_flags ;
147 Int_t gpsSec ;
148 Int_t gpsNanosec ;
149 Int_t nChTot ;
150 Int_t nChAvail ;
151 Int_t nHitTot ;
152 Int_t nHit ;
153 Int_t nHitSP10 ;
154 Int_t nHitSP20 ;
155 Int_t nTankTot ;
156 Int_t nTankAvail ;
157 Int_t nTankHitTot ;
158 Int_t nTankHit ;
159 Int_t windowHits ;
160 Int_t angleFitStatus;
161 Int_t planeNDOF ;
162 Int_t coreFitStatus;
163 Int_t CxPE40PMT ;
164 Int_t CxPE40XnCh ;
165 Double_t mPFnHits;
166 Double_t mPFnPlanes;
167 Double_t mPFp0nAssign;
168 Double_t mPFp1nAssign;
169 Double_t nHitMPF ;
170 Int_t coreFiduScale;
171 Double_t coreFiduScaleOR ;
172 Double_t coreLoc ;
173 Double_t zenithAngle ;
174 Double_t azimuthAngle ;
175 Double_t dec ;
176 Double_t ra ;
177 Double_t planeChi2 ;
178 Double_t coreX ;
179 Double_t coreY ;
180 Double_t logCoreAmplitude ;
181 Double_t coreFitUnc ;
182 Double_t logNNEnergy ;
183 Double_t fAnnulusCharge0;
184 Double_t fAnnulusCharge1;
185 Double_t fAnnulusCharge2;
186 Double_t fAnnulusCharge3;
187 Double_t fAnnulusCharge4;
188 Double_t fAnnulusCharge5;
189 Double_t fAnnulusCharge6;
190 Double_t fAnnulusCharge7;
191 Double_t fAnnulusCharge8;
192 Double_t protonlheEnergy;
193 Double_t protonlheLLH ;
194 Double_t gammalheEnergy ;
195 Double_t gammalheLLH ;
196 Int_t chargeFiduScale50 ;
197 Int_t chargeFiduScale70 ;
198 Int_t chargeFiduScale90 ;
199 Double_t logMaxPE;
200 Double_t logNPE ;
201 Int_t CxPE40 ;
202 Double_t CxPE40SPTime;
203 Double_t LDFAge ;
204 Double_t LDFAmp ;
205 Double_t LDFChi2 ;
206 Double_t logGP ;
207 Double_t mPFp0Weight;
208 Double_t mPFp0toangleFit;
209 Double_t mPFp1Weight ;
210 Double_t mPFp1toangleFit;
211 Double_t PINC ;
212 Double_t disMax ;
213 Double_t TankLHR ;
214 Double_t LHLatDistFitXmax ;
215 Double_t LHLatDistFitEnergy;
216 Double_t LHLatDistFitGoF ;
217 Double_t probaGBT ;
218
219};
220
221struct tel_struct
222{
223 MTime *mTime = 0;
224 MRawEvtHeader *mRawEvtHeader = 0;
225 //MSignalCam *mSignalCam = 0;
226 //MRawEvtHeader *mRawEvtHeader = 0;
227 //MSoftwareTrigger *mSoftwareTrigger = 0;
228 //MRawBoardsFACT *mRawBoardsFACT = 0;
229 Double_t cal_unix; // calibrated unix timestamp based on board's counter
230 UInt_t DAQ_id; // event's DAQ ID (fDAQEvtNumber)
231
232};
233
234struct sync_struct
235{
236 Double_t ts;
237 UInt_t id;
238};
239
240//Earlier timestamps first in vector:
241bool sort_time (sync_struct str_1, sync_struct str_2) { return (str_1.ts<str_2.ts); }
242
243int synchron(TString HAWC_path, TString HEye_path, TString Outpath)
244{
245
246 /* INPUT: - path to reconstructed HAWC's data
247 - path to calibrated HAWC's Eye data
248 - output path
249 OUTPUT:- 1 plot file contains the plot "Time Difference vs. Time Shift"
250 - 1 root file contains the sync data and sync informations
251 */
252
253
254 const Double_t HAWC_toffset = 315964781.0-0.048100; //HAWC Offset to UTC (01.01.1980), - best shift
255 const Double_t TimeDiff_Cut = 0.010; // Maximal time [s] to be accepted as synchonised
256 const Double_t Start_Shift = -2.0; // [s]
257 const Double_t Stop_Shift = 2.0; // [s]
258 const Double_t shift_step = 0.0001; // [s]
259
260 TChain tHawc("XCDF");
261 TChain tHEye("Events");
262
263 tHawc.Add(HAWC_path);
264 tHEye.Add(HEye_path);
265
266
267 //--------------------- HAWC's data ---------------------
268 rec_struct rec;
269
270 tHawc.SetBranchAddress("rec.status" , &rec.status) ;
271 tHawc.SetBranchAddress("rec.version", &rec.version) ;
272 tHawc.SetBranchAddress("rec.eventID", &rec.eventID) ;
273 tHawc.SetBranchAddress("rec.runID", &rec.runID) ;
274 tHawc.SetBranchAddress("rec.timeSliceID", &rec.timeSliceID);
275 tHawc.SetBranchAddress("rec.trigger_flags", &rec.trigger_flags);
276 tHawc.SetBranchAddress("rec.event_flags",&rec.event_flags) ;
277 tHawc.SetBranchAddress("rec.gtc_flags",&rec.gtc_flags) ;
278 tHawc.SetBranchAddress("rec.gpsSec",&rec.gpsSec) ;
279 tHawc.SetBranchAddress("rec.gpsNanosec",&rec.gpsNanosec) ;
280 tHawc.SetBranchAddress("rec.nChTot",&rec.nChTot) ;
281 tHawc.SetBranchAddress("rec.nChAvail",&rec.nChAvail) ;
282 tHawc.SetBranchAddress("rec.nHitTot",&rec.nHitTot) ;
283 tHawc.SetBranchAddress("rec.nHit",&rec.nHit) ;
284 tHawc.SetBranchAddress("rec.nHitSP10",&rec.nHitSP10) ;
285 tHawc.SetBranchAddress("rec.nHitSP20",&rec.nHitSP20) ;
286 tHawc.SetBranchAddress("rec.nTankTot",&rec.nTankTot) ;
287 tHawc.SetBranchAddress("rec.nTankAvail",&rec.nTankAvail) ;
288 tHawc.SetBranchAddress("rec.nTankHitTot",&rec.nTankHitTot) ;
289 tHawc.SetBranchAddress("rec.nTankHit",&rec.nTankHit) ;
290 tHawc.SetBranchAddress("rec.windowHits",&rec.windowHits) ;
291 tHawc.SetBranchAddress("rec.angleFitStatus",&rec.angleFitStatus);
292 tHawc.SetBranchAddress("rec.planeNDOF",&rec.planeNDOF) ;
293 tHawc.SetBranchAddress("rec.coreFitStatus",&rec.coreFitStatus);
294 tHawc.SetBranchAddress("rec.CxPE40PMT",&rec.CxPE40PMT) ;
295 tHawc.SetBranchAddress("rec.CxPE40XnCh",&rec.CxPE40XnCh);
296 tHawc.SetBranchAddress("rec.mPFnHits",&rec.mPFnHits);
297 tHawc.SetBranchAddress("rec.mPFnPlanes",&rec.mPFnPlanes);
298 tHawc.SetBranchAddress("rec.mPFp0nAssign",&rec.mPFp0nAssign);
299 tHawc.SetBranchAddress("rec.mPFp1nAssign",&rec.mPFp1nAssign);
300 tHawc.SetBranchAddress("rec.nHitMPF",&rec.nHitMPF) ;
301 tHawc.SetBranchAddress("rec.coreFiduScale",&rec.coreFiduScale);
302 tHawc.SetBranchAddress("rec.coreFiduScaleOR",&rec.coreFiduScaleOR) ;
303 tHawc.SetBranchAddress("rec.coreLoc",&rec.coreLoc) ;
304 tHawc.SetBranchAddress("rec.zenithAngle",&rec.zenithAngle) ;
305 tHawc.SetBranchAddress("rec.azimuthAngle",&rec.azimuthAngle) ;
306 tHawc.SetBranchAddress("rec.dec",&rec.dec) ;
307 tHawc.SetBranchAddress("rec.ra",&rec.ra) ;
308 tHawc.SetBranchAddress("rec.planeChi2",&rec.planeChi2) ;
309 tHawc.SetBranchAddress("rec.coreX",&rec.coreX) ;
310 tHawc.SetBranchAddress("rec.coreY",&rec.coreY) ;
311 tHawc.SetBranchAddress("rec.logCoreAmplitude",&rec.logCoreAmplitude) ;
312 tHawc.SetBranchAddress("rec.coreFitUnc",&rec.coreFitUnc) ;
313 tHawc.SetBranchAddress("rec.logNNEnergy",&rec.logNNEnergy) ;
314 tHawc.SetBranchAddress("rec.fAnnulusCharge0",&rec.fAnnulusCharge0);
315 tHawc.SetBranchAddress("rec.fAnnulusCharge1",&rec.fAnnulusCharge1);
316 tHawc.SetBranchAddress("rec.fAnnulusCharge2",&rec.fAnnulusCharge2);
317 tHawc.SetBranchAddress("rec.fAnnulusCharge3",&rec.fAnnulusCharge3);
318 tHawc.SetBranchAddress("rec.fAnnulusCharge4",&rec.fAnnulusCharge4);
319 tHawc.SetBranchAddress("rec.fAnnulusCharge5",&rec.fAnnulusCharge5);
320 tHawc.SetBranchAddress("rec.fAnnulusCharge6",&rec.fAnnulusCharge6);
321 tHawc.SetBranchAddress("rec.fAnnulusCharge7",&rec.fAnnulusCharge7);
322 tHawc.SetBranchAddress("rec.fAnnulusCharge8",&rec.fAnnulusCharge8);
323 tHawc.SetBranchAddress("rec.protonlheEnergy",&rec.protonlheEnergy);
324 tHawc.SetBranchAddress("rec.protonlheLLH",&rec.protonlheLLH) ;
325 tHawc.SetBranchAddress("rec.gammalheEnergy",&rec.gammalheEnergy) ;
326 tHawc.SetBranchAddress("rec.gammalheLLH",&rec.gammalheLLH) ;
327 tHawc.SetBranchAddress("rec.chargeFiduScale50",&rec.chargeFiduScale50) ;
328 tHawc.SetBranchAddress("rec.chargeFiduScale70",&rec.chargeFiduScale70) ;
329 tHawc.SetBranchAddress("rec.chargeFiduScale90",&rec.chargeFiduScale90) ;
330 tHawc.SetBranchAddress("rec.logMaxPE",&rec.logMaxPE);
331 tHawc.SetBranchAddress("rec.logNPE",&rec.logNPE) ;
332 tHawc.SetBranchAddress("rec.CxPE40",&rec.CxPE40) ;
333 tHawc.SetBranchAddress("rec.CxPE40SPTime",&rec.CxPE40SPTime);
334 tHawc.SetBranchAddress("rec.LDFAge",&rec.LDFAge) ;
335 tHawc.SetBranchAddress("rec.LDFAmp",&rec.LDFAmp) ;
336 tHawc.SetBranchAddress("rec.LDFChi2",&rec.LDFChi2) ;
337 tHawc.SetBranchAddress("rec.logGP",&rec.logGP) ;
338 tHawc.SetBranchAddress("rec.mPFp0Weight",&rec.mPFp0Weight);
339 tHawc.SetBranchAddress("rec.mPFp0toangleFit",&rec.mPFp0toangleFit);
340 tHawc.SetBranchAddress("rec.mPFp1Weight",&rec.mPFp1Weight) ;
341 tHawc.SetBranchAddress("rec.mPFp1toangleFit",&rec.mPFp1toangleFit);
342 tHawc.SetBranchAddress("rec.PINC",&rec.PINC) ;
343 tHawc.SetBranchAddress("rec.disMax",&rec.disMax) ;
344 tHawc.SetBranchAddress("rec.TankLHR",&rec.TankLHR) ;
345 tHawc.SetBranchAddress("rec.LHLatDistFitXmax",&rec.LHLatDistFitXmax) ;
346 tHawc.SetBranchAddress("rec.LHLatDistFitEnergy",&rec.LHLatDistFitEnergy);
347 tHawc.SetBranchAddress("rec.LHLatDistFitGoF",&rec.LHLatDistFitGoF) ;
348 tHawc.SetBranchAddress("rec.probaGBT",&rec.probaGBT) ;
349
350 //--------------------- HawcEye's data ---------------------
351 tel_struct tel;
352
353 tHEye.SetBranchAddress("MTime.", &(tel.mTime));
354 //tHEye.SetBranchAddress("MSignalCam.", &(tel.mSignalCam));
355 tHEye.SetBranchAddress("MRawEvtHeader.", &(tel.mRawEvtHeader));
356 //tHEye.SetBranchAddress("MSoftwareTrigger.", &(tel.mSoftwareTrigger));
357 //tHEye.SetBranchAddress("MRawBoardsFACT.", &(tel.mRawBoardsFACT));
358
359 // Calibrated Unix timestamp based on FAD's counter
360 vector<Double_t> unix_s = calibration(HEye_path);
361
362
363 //--------------------------------------------------------------------------------------------
364 // Read in the ID list and timestamps and sort them
365 vector<sync_struct> vts_HEye;
366 vector<sync_struct> vts_Hawc;
367 sync_struct tsync;
368
369 printf("Reading telescope data....\n");
370
371 Long64_t nHE_events = tHEye.GetEntries();
372 cout << "Number of HawcEye's Events: " << nHE_events << endl;
373
374 for(Int_t i = 0; i < nHE_events; i++)
375 {
376 tHEye.GetEntry(i);
377 tsync.id = i;
378 tsync.ts = unix_s[i];
379 vts_HEye.push_back(tsync);
380 }
381
382 printf("Sorting telescope data....\n");
383 sort(vts_HEye.begin(), vts_HEye.end(),sort_time);
384
385
386 printf("Reading HAWC data....\n");
387
388 Long64_t nHawc_events = tHawc.GetEntries();
389 cout << "Number of Hawc's Events: " << nHawc_events << endl;
390 for(Int_t i=0;i < nHawc_events; i++)
391 {
392 tHawc.GetEntry(i);
393 tsync.id = i;
394 tsync.ts = HAWC_toffset + rec.gpsSec + ((Double_t)rec.gpsNanosec)/1000000000.0;
395 vts_Hawc.push_back(tsync);
396 }
397 printf("Sorting HAWC data....\n");
398 sort(vts_Hawc.begin(), vts_Hawc.end(),sort_time);
399
400 // Find the best sync offset
401 printf("HAWC Data:\nFirst event: %f s Last event: %f s (Trigger: %ld, Rate: %f)\n\n",
402 vts_Hawc[0].ts, vts_Hawc[vts_Hawc.size()-1].ts, vts_Hawc.size(),
403 vts_Hawc.size() / (vts_Hawc[vts_Hawc.size()-1].ts - vts_Hawc[0].ts));
404
405 printf("HEye Data:\nFirst event: %f s Last event: %f s (Trigger: %ld, Rate: %f)\n\n",
406 vts_HEye[0].ts, vts_HEye[vts_HEye.size()-1].ts, vts_HEye.size(),
407 vts_HEye.size()/(vts_HEye[vts_HEye.size()-1].ts-vts_HEye[0].ts));
408
409
410 if (vts_Hawc[0].ts>vts_HEye[nHE_events-1].ts)
411 {
412 printf("NO OVERLAP: HAWC data later, choose later HEye file!\n");
413 return 1;
414 }
415 if (vts_Hawc[nHawc_events-1].ts<vts_HEye[0].ts)
416 {
417 printf("NO OVERLAP: HAWC data earlier, choose earlier HEye file!\n");
418 return 2;
419 }
420
421 printf("Cutting HAWC data to fit shift....\n");
422
423 Int_t first_shift = 0; // index of first Hawc's event with time > first HE's time - Start_Shift
424 Int_t last_shift = 0; // index of last Hawc's event with time < last HE's time - Stop_Shift
425 Int_t first = 0; // index of first Hawc's event with trigger time > first HE's time
426 Int_t last = 0 ; // index of last Hawc's event with trigger time < last HE's time
427
428 for (Long64_t i = 0; i < nHawc_events; i++)
429 {
430 if (vts_Hawc[i].ts > vts_HEye[0].ts)
431 {
432 first = i;
433 }
434 if (vts_Hawc[i].ts > vts_HEye[0].ts - Start_Shift)
435 {
436 first_shift = i;
437 break;
438 }
439 }
440
441 for(Long64_t i = nHawc_events-1; i > 0; i--)
442 {
443 if (vts_Hawc[i].ts < vts_HEye[nHE_events-1].ts)
444 {
445 last = i;
446 }
447 if (vts_Hawc[i].ts < vts_HEye[nHE_events-1].ts - Stop_Shift)
448 {
449 last_shift = i;
450 break;
451 }
452 }
453
454 printf("HAWC Data (reduced):\nFirst event: %f s Last event: %f s\n\n",
455 vts_Hawc[first_shift].ts, vts_Hawc[last_shift].ts);
456
457 printf("HEye Data:\nFirst event: %f s Last event: %f s\n\n",
458 vts_HEye[0].ts, vts_HEye[vts_HEye.size()-1].ts);
459
460 //---------------------- Define the Output TTree ----------------------
461 TFile fOut(Outpath, "RECREATE");
462
463 TTree tOut1("Hawc_Data", "Hawc_Data"); // Hawc's data
464 TTree tOut2("HE_Data", "HE_Data"); // Telescope's data
465 TTree tOut3("timediff_tot", "timediff_tot"); // Synchronisation's info
466 TTree tOut4("timediff_pass", "timediff_pass");
467 TTree tOut5("SyncInfo", "SyncInfo");
468
469 // Hawc Data
470 tOut1.Branch("rec.status" , &rec.status) ;
471 tOut1.Branch("rec.version", &rec.version) ;
472 tOut1.Branch("rec.eventID", &rec.eventID) ;
473 tOut1.Branch("rec.runID", &rec.runID) ;
474 tOut1.Branch("rec.timeSliceID", &rec.timeSliceID);
475 tOut1.Branch("rec.trigger_flags", &rec.trigger_flags);
476 tOut1.Branch("rec.event_flags",&rec.event_flags) ;
477 tOut1.Branch("rec.gtc_flags",&rec.gtc_flags) ;
478 tOut1.Branch("rec.gpsSec",&rec.gpsSec) ;
479 tOut1.Branch("rec.gpsNanosec",&rec.gpsNanosec) ;
480 tOut1.Branch("rec.nChTot",&rec.nChTot) ;
481 tOut1.Branch("rec.nChAvail",&rec.nChAvail) ;
482 tOut1.Branch("rec.nHitTot",&rec.nHitTot) ;
483 tOut1.Branch("rec.nHit",&rec.nHit) ;
484 tOut1.Branch("rec.nHitSP10",&rec.nHitSP10) ;
485 tOut1.Branch("rec.nHitSP20",&rec.nHitSP20) ;
486 tOut1.Branch("rec.nTankTot",&rec.nTankTot) ;
487 tOut1.Branch("rec.nTankAvail",&rec.nTankAvail) ;
488 tOut1.Branch("rec.nTankHitTot",&rec.nTankHitTot) ;
489 tOut1.Branch("rec.nTankHit",&rec.nTankHit) ;
490 tOut1.Branch("rec.windowHits",&rec.windowHits) ;
491 tOut1.Branch("rec.angleFitStatus",&rec.angleFitStatus);
492 tOut1.Branch("rec.planeNDOF",&rec.planeNDOF) ;
493 tOut1.Branch("rec.coreFitStatus",&rec.coreFitStatus);
494 tOut1.Branch("rec.CxPE40PMT",&rec.CxPE40PMT) ;
495 tOut1.Branch("rec.CxPE40XnCh",&rec.CxPE40XnCh);
496 tOut1.Branch("rec.mPFnHits",&rec.mPFnHits);
497 tOut1.Branch("rec.mPFnPlanes",&rec.mPFnPlanes);
498 tOut1.Branch("rec.mPFp0nAssign",&rec.mPFp0nAssign);
499 tOut1.Branch("rec.mPFp1nAssign",&rec.mPFp1nAssign);
500 tOut1.Branch("rec.nHitMPF",&rec.nHitMPF) ;
501 tOut1.Branch("rec.coreFiduScale",&rec.coreFiduScale);
502 tOut1.Branch("rec.coreFiduScaleOR",&rec.coreFiduScaleOR) ;
503 tOut1.Branch("rec.coreLoc",&rec.coreLoc) ;
504 tOut1.Branch("rec.zenithAngle",&rec.zenithAngle) ;
505 tOut1.Branch("rec.azimuthAngle",&rec.azimuthAngle) ;
506 tOut1.Branch("rec.dec",&rec.dec) ;
507 tOut1.Branch("rec.ra",&rec.ra) ;
508 tOut1.Branch("rec.planeChi2",&rec.planeChi2) ;
509 tOut1.Branch("rec.coreX",&rec.coreX) ;
510 tOut1.Branch("rec.coreY",&rec.coreY) ;
511 tOut1.Branch("rec.logCoreAmplitude",&rec.logCoreAmplitude) ;
512 tOut1.Branch("rec.coreFitUnc",&rec.coreFitUnc) ;
513 tOut1.Branch("rec.logNNEnergy",&rec.logNNEnergy) ;
514 tOut1.Branch("rec.fAnnulusCharge0",&rec.fAnnulusCharge0);
515 tOut1.Branch("rec.fAnnulusCharge1",&rec.fAnnulusCharge1);
516 tOut1.Branch("rec.fAnnulusCharge2",&rec.fAnnulusCharge2);
517 tOut1.Branch("rec.fAnnulusCharge3",&rec.fAnnulusCharge3);
518 tOut1.Branch("rec.fAnnulusCharge4",&rec.fAnnulusCharge4);
519 tOut1.Branch("rec.fAnnulusCharge5",&rec.fAnnulusCharge5);
520 tOut1.Branch("rec.fAnnulusCharge6",&rec.fAnnulusCharge6);
521 tOut1.Branch("rec.fAnnulusCharge7",&rec.fAnnulusCharge7);
522 tOut1.Branch("rec.fAnnulusCharge8",&rec.fAnnulusCharge8);
523 tOut1.Branch("rec.protonlheEnergy",&rec.protonlheEnergy);
524 tOut1.Branch("rec.protonlheLLH",&rec.protonlheLLH) ;
525 tOut1.Branch("rec.gammalheEnergy",&rec.gammalheEnergy) ;
526 tOut1.Branch("rec.gammalheLLH",&rec.gammalheLLH) ;
527 tOut1.Branch("rec.chargeFiduScale50",&rec.chargeFiduScale50) ;
528 tOut1.Branch("rec.chargeFiduScale70",&rec.chargeFiduScale70) ;
529 tOut1.Branch("rec.chargeFiduScale90",&rec.chargeFiduScale90) ;
530 tOut1.Branch("rec.logMaxPE",&rec.logMaxPE);
531 tOut1.Branch("rec.logNPE",&rec.logNPE) ;
532 tOut1.Branch("rec.CxPE40",&rec.CxPE40) ;
533 tOut1.Branch("rec.CxPE40SPTime",&rec.CxPE40SPTime);
534 tOut1.Branch("rec.LDFAge",&rec.LDFAge) ;
535 tOut1.Branch("rec.LDFAmp",&rec.LDFAmp) ;
536 tOut1.Branch("rec.LDFChi2",&rec.LDFChi2) ;
537 tOut1.Branch("rec.logGP",&rec.logGP) ;
538 tOut1.Branch("rec.mPFp0Weight",&rec.mPFp0Weight);
539 tOut1.Branch("rec.mPFp0toangleFit",&rec.mPFp0toangleFit);
540 tOut1.Branch("rec.mPFp1Weight",&rec.mPFp1Weight) ;
541 tOut1.Branch("rec.mPFp1toangleFit",&rec.mPFp1toangleFit);
542 tOut1.Branch("rec.PINC",&rec.PINC) ;
543 tOut1.Branch("rec.disMax",&rec.disMax) ;
544 tOut1.Branch("rec.TankLHR",&rec.TankLHR) ;
545 tOut1.Branch("rec.LHLatDistFitXmax",&rec.LHLatDistFitXmax) ;
546 tOut1.Branch("rec.LHLatDistFitEnergy",&rec.LHLatDistFitEnergy);
547 tOut1.Branch("rec.LHLatDistFitGoF",&rec.LHLatDistFitGoF) ;
548 tOut1.Branch("rec.probaGBT",&rec.probaGBT) ;
549
550 // Telescope's Data
551 tOut2.Branch("cal_unix", &tel.cal_unix);
552 tOut2.Branch("DAQ_id", &tel.DAQ_id);
553
554 // Synchronisation Info
555 Double_t min_miss = 100; // [s]
556 Double_t best_shift = 0.0; // best time shift
557 Double_t timediff_tot; // time difference after best shift - for all events
558 Double_t timediff_passed; // time difference after best shift - for events within TimeDiff_Cut
559
560 tOut3.Branch("timediff_tot", &timediff_tot); // time diff between neighbor events after shift for all events
561 tOut4.Branch("timediff_passed", &timediff_passed); // time diff between neighbors after shift for sync events
562 tOut5.Branch("min_miss", &min_miss); // mean time diff with best shift
563 tOut5.Branch("best_shift", &best_shift); // best time shift
564
565 //----------------------------------------------------------------------------
566
567 Double_t total_miss = 0.0; // total time difference sum( |t_Hawc(i) - t_HE(i)|) [s]
568 Double_t mean_miss = 0.0; // total_miss/n_Hawc
569
570 Double_t diff_p = 0.0;
571 Double_t diff_n = 0.0;
572 Double_t diff_min = 0.0;
573 Int_t temp_ind = 0;
574
575 Int_t last_index = 0;
576
577
578 // For the plot Mean Time Difference vs. Time Shift
579 vector<Double_t> gr_tShift, gr_tDiff;
580 Int_t gr_n = (Stop_Shift-Start_Shift)/shift_step;
581
582 printf("Doing shift match....\n");
583
584 // Search for best shift with min. mean timem difference
585 for (Double_t shift = Start_Shift; shift < Stop_Shift; shift += shift_step)
586 {
587 for (Int_t i = first_shift; i < last_shift; i++)
588 {
589
590 //1. Find the nearest neighbor in vts_HEye to vts_Hawc[i]+shift:
591 while (1)
592 {
593 if(vts_Hawc[i].ts + shift > vts_HEye[last_index + 1].ts)
594 last_index++;
595 else
596 break;
597 }
598 //2. calculate the difference
599 diff_n = vts_Hawc[i].ts+shift - vts_HEye[last_index].ts;
600 diff_p = vts_HEye[last_index+1].ts - vts_Hawc[i].ts-shift;
601 diff_min = diff_n > diff_p ? diff_p : diff_n;
602
603 //3. Add to total_miss
604 total_miss += diff_min;
605 }
606 mean_miss = total_miss/(last_shift-first_shift);
607
608 gr_tShift.push_back(shift);
609 gr_tDiff.push_back(mean_miss);
610
611 if (min_miss > mean_miss)
612 {
613 min_miss = mean_miss;
614 best_shift = shift;
615 }
616 total_miss = 0.0;
617 last_index = 0;
618
619 }
620 tOut5.Fill();
621
622 printf("->Minimum sync miss: %f @ a shift of %fs\n", min_miss, best_shift);
623
624
625 TH1F timediff_neighbor("timediff_neighbor","timediff_neighbor",1000,-0.01,0.01);
626
627 // Shift all events by the found best shift
628 for(Int_t i = first; i < last; i++)
629 {
630 //1. Find the nearest neighbor in vts_HEye to vts_Hawc[i]+best_shift:
631 while(1)
632 {
633 if(vts_Hawc[i].ts + best_shift > vts_HEye[last_index+1].ts)
634 last_index++;
635 else
636 break;
637 }
638
639 //2. calculate the difference
640 diff_n = vts_Hawc[i].ts + best_shift - vts_HEye[last_index].ts;
641 diff_p = vts_HEye[last_index+1].ts - vts_Hawc[i].ts-best_shift;
642
643 if (diff_p > diff_n)
644 {
645 diff_min = diff_n;
646 temp_ind = last_index;
647 }
648 else
649 {
650 diff_min = diff_p;
651 temp_ind = last_index+1;
652 }
653
654 timediff_tot = diff_min;
655
656 if (diff_min < TimeDiff_Cut)
657 {
658 tHEye.GetEntry(vts_HEye[temp_ind].id);
659 tHawc.GetEntry(vts_Hawc[i].id);
660
661 tel.cal_unix = vts_HEye[temp_ind].ts;
662 tel.DAQ_id = tel.mRawEvtHeader->GetDAQEvtNumber();
663
664 timediff_passed = diff_min;
665
666 tOut1.Fill();
667 tOut2.Fill();
668 tOut4.Fill();
669
670 }
671 timediff_neighbor.Fill(diff_min); //Fill historgramm with diff
672 tOut3.Fill();
673 }
674
675 //---------------------------Plots and Output files---------------------------
676
677 tOut1.Write();
678 tOut2.Write();
679 tOut3.Write();
680 tOut4.Write();
681 tOut5.Write();
682 fOut.Close();
683
684 TCanvas *c1 = new TCanvas("timediff_shift", "Synchronisation", 1200,800);
685 c1->SetGridx();
686 c1->SetGridy();
687
688 TGraph *gr_shift = new TGraph(gr_n, &gr_tShift[0], &gr_tDiff[0]);
689 gr_shift->SetTitle("Synchronisation");
690 gr_shift->GetXaxis()->SetTitle("Time Shift [s]");
691 gr_shift->GetYaxis()->SetTitle("Mean time difference [s]");
692 gr_shift->Draw();
693
694 TCanvas *c2 = new TCanvas("timediff_neighbor", "timediff_neighbor", 1200, 800);
695 c2->SetGridx();
696 c2->SetGridy();
697 timediff_neighbor.GetXaxis()->SetTitle("Time Difference [s]");
698 timediff_neighbor.GetYaxis()->SetTitle("Counts");
699 timediff_neighbor.DrawClone();
700
701 // Write out files
702
703 TFile fOutPlots(Outpath+"_syncplots.root","RECREATE");
704 c1->Write();
705 c2->Write();
706 fOutPlots.Close();
707
708 return 0;
709
710}
Note: See TracBrowser for help on using the repository browser.