Index: trunk/Mars/hawc/callisto.C =================================================================== --- trunk/Mars/hawc/callisto.C (revision 19950) +++ trunk/Mars/hawc/callisto.C (revision 19951) @@ -27,5 +27,5 @@ To run the macro from the command line (assuming you are in a directory - Mars/build where you have built your Mars environment) ou can do + Mars/build where you have built your Mars environment) you can do root ../hawc/callisto.C\(\"20191002_018.fits.fz\",\"20191002_013.drs.fits\"\) @@ -54,5 +54,5 @@ // mapping file (found in Mars/hawc) - const char *mmap = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL; + const char *mmap = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL; //const char *pulse_template = "template-pulse.root"; @@ -76,5 +76,6 @@ // Calibration constant (for the moment a single constant to // convert the extracted charge to photo-electrons) - double scale = 1./22.553;//0.2; + //double scale = 1./22.553;//0.2; + double scale = 1./17.557; // ------------------------------------------------------ @@ -317,4 +318,5 @@ write5.AddContainer("MSimSourcePos", "Events", kFALSE); write5.AddContainer("MTime", "Events", kFALSE); + write5.AddContainer("MRawBoardsFACT", "Events", kFALSE); write5.AddContainer("MSrcPosCam", "Events", kFALSE); Index: trunk/Mars/hawc/extract_pulse.C =================================================================== --- trunk/Mars/hawc/extract_pulse.C (revision 19950) +++ trunk/Mars/hawc/extract_pulse.C (revision 19951) @@ -60,5 +60,5 @@ }; -//Storage class to kKeep a list of the extracted single +//Storage class to keep a list of the extracted single class MSingles : public MParContainer, public MCamEvent { @@ -68,8 +68,8 @@ public: - MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30) + MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(40) { fName = name ? name : "MSingles"; - fName = title ? title : "Storeage for vector of singles"; + fName = title ? title : "Storage for vector of singles"; } @@ -109,5 +109,5 @@ public: - MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10) + MHBaseline() : fPedestal(0), fSkipStart(5), fSkipEnd(10) { fName = "MHBaseline"; @@ -185,4 +185,12 @@ for (int x=0; x<64; x++) { + //exclude the pixels without cones + if (x == 61) + continue; + if (x == 62) + continue; + if (x == 63) + continue; + // Get the corresponding slice from the histogram TH1D *hist = fBaseline.ProjectionY("proj", x+1, x+1); @@ -286,8 +294,11 @@ class MHSingles : public MH { - TH2F fSignalAll; + TH2F fSignals; // all signals, without normalization + TH2F fOsicllation; + + TH2F fSignalAll; // Histogram of all signals TH2F fSignalOut; - TH2F fSignal1; - TH2F fSignal2; + TH2F fSignal1; // Histogram of signals with 1pe + TH2F fSignal2; // Histogram of signals with 2pe TH2F fSignal3; TH2F fSignal4; @@ -320,4 +331,25 @@ // Setup histograms + + fSignals.SetName("Signals"); + fSignals.SetTitle(""); + fSignals.SetXTitle("Time [ns]"); + fSignals.SetYTitle("Amplitude [mV]"); + fSignals.GetXaxis()->CenterTitle(); + fSignals.GetXaxis()->SetTitleSize(0.045); + fSignals.GetYaxis()->SetTitleSize(0.045); + fSignals.SetDirectory(NULL); + fSignals.SetStats(kFALSE); + + fOsicllation.SetName("Oscillations"); + fOsicllation.SetTitle(""); + fOsicllation.SetXTitle("Oscillation's Amplitude [mV]"); + fOsicllation.SetYTitle("Pulse's Amplitude [mV]"); + fOsicllation.GetXaxis()->CenterTitle(); + fOsicllation.GetXaxis()->SetTitleSize(0.045); + fOsicllation.GetYaxis()->SetTitleSize(0.045); + fOsicllation.SetDirectory(NULL); + fOsicllation.SetStats(kFALSE); + fSignalAll.SetName("SignalAll"); fSignalAll.SetTitle(""); @@ -498,7 +530,16 @@ //const Int_t w = fSingles->GetIntegrationWindow(); - MBinning binsx, binsy; + MBinning binsx, binsy, bin_time, bin_pulse, binosc, binpulse; binsx.SetEdges(1024, -512, 512); binsy.SetEdges( 800, -20, 20); + + bin_time.SetEdges(1024, -512, 512); + bin_pulse.SetEdges(240, -60, 60); //bin width is about 1 dac count = 0.5mV + + binosc.SetEdges(140, -65, 5); + binpulse.SetEdges(140, -5, 65); + + MH::SetBinning(fSignals, bin_time, bin_pulse); + MH::SetBinning(fOsicllation, binosc, binpulse); MH::SetBinning(fSignalAll, binsx, binsy); @@ -579,8 +620,8 @@ continue; - if (idx!=0) + if (idx!=43) continue; - const double gain = hgain->GetBinContent(idx+1)/30; + const double gain = hgain->GetBinContent(idx+1)/40; const double base = hbase->GetBinContent(idx+1); @@ -594,19 +635,58 @@ // Fill pulse shape const double tm = result[0].fTime; - const double sig = result[0].fSignal/30 - base; + const double sig = result[0].fSignal/40 - base; bool ok = true;//tm>3 && (ptr[int(tm)-3]-ped-base)/gain > -5; - double mean = 0; - - for (int s=20; s max_pulse) + max_pulse = ptr[s]-ped-base; + max_pulse_pos = s; + } + //cout << "________Check if the filled histograms are correct________" << endl; + //cout << "Arrival time: " << tm << ",Amplitude: " << ptr[t_arrival]-ped-base << endl; + //cout << "Pulse's position: "<< max_pulse_pos << ",Max. Pulse's Amplitude: " << ptr[max_pulse_pos]-ped-base << endl; + + // Find maximal oscillation's amplitude + double max_osc = 0; + for (int osc = max_pulse_pos; osc < max_pulse_pos+40; osc++) + { + if (ptr[osc]-ped-base < max_osc) + { + max_osc = ptr[osc]-ped-base; + } + else + continue; + } + fOsicllation.Fill(max_osc, max_pulse); + + // Set conditions for 1pe signals + for (int s=5; s-100 && x<-10) - mean += y/90; - - if (y<-3) + //if (x>-100 && x<-10) + // mean += y/90; + + if (x>5 && y>2.5+x*-0.01) + ok = false; + if (y<-1.5) ok = false; if (x>5 && x<40 && y<=x/40*-0.5) @@ -614,10 +694,13 @@ if (x>-150 && x<-5 && y>1.5) ok = false; - if (y>2.5) + if (x<0 && y > 1.2) ok = false; + } - // Nothing at or bleow zero from 5 to 40machen - + // Fill nPE signals into histograms + + // Set conditions for signals with n PE + // Nothing at or below zero from 5 to 40machen bool ok1 = sig>gain*0.5 && siggain*1.5 && sigSetFrameBorderMode(0); fSignal8.Draw("colz"); - //fProf8.Draw("same"); + fProf8.Draw("same"); } @@ -836,6 +922,6 @@ MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0), fExtractionRange(64), - fNumAverage(10), fStartSkip(20), fEndSkip(10), - fIntegrationSize(3*10), fMaxSearchWindow(30) + fNumAverage(10), fStartSkip(10), fEndSkip(10), + fIntegrationSize(4*10), fMaxSearchWindow(20) { } @@ -877,5 +963,5 @@ const size_t navg = fNumAverage; - const float threshold = fThreshold; + const float threshold = fThreshold; const UInt_t start_skip = fStartSkip; const UInt_t end_skip = fEndSkip; @@ -934,4 +1020,13 @@ if (k_max == i+5 || k_max == i + max_search_window-1) continue; + + //Make sure the pulse is isolated (area after the maximum is empty) + UInt_t k_falling = k_max+200; + for (UInt_t k_after = k_falling; k_after < k_max +250; k_after++) + { + if (val[k_after] > val[k_falling] && + val[k_after + fNumAverage/2] > val[k_falling]) + continue; + } //search for half maximum before maximum @@ -946,4 +1041,5 @@ } } + // Check if the finding makes sense if (k_half_max+integration_size > roi-navg-end_skip) @@ -954,6 +1050,8 @@ continue; - // FIXME: Evaluate arrival time more precisely!!! - // FIXME: Get a better integration window + //cout << "________Check if the finding makes sense________" << endl; + //cout << "Pulse's position: " << k_max << ", Max. Pulse's Amplitude: " << val[k_max] << endl; + //cout << "Half maximum: " << k_half_max << ", Pulse at half maximum: " << val[k_half_max] << endl; + // Compile "single" @@ -970,4 +1068,5 @@ // Add single to list result.push_back(single); + break; //accept only one pulse per event // skip falling edge @@ -985,15 +1084,15 @@ Int_t max_dist = 14, Double_t threshold = 5, - const char *runfile = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_", + const char *runfile = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_", int firstRunID = 6, int lastRunID = 6, - const char *drsfile = "/media/tbretz/d84f26d3-fc9b-4b30-8cfd-85246846dd1a/hawcseye01/raw/2019/10/27/20191027_003.drs.fits", + const char *drsfile = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/20191027_004.drs.fits", const char *outpath = "." ) { // ====================================================== - if (!pmap.Read("../hawc/FAMOUSmap171218.txt")) - { - gLog << err << "../hawc/FAMOUSmap171218.txt not found." << endl; + if (!pmap.Read("/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt")) + { + gLog << err << "HAWCsEyemap181214.txt not found." << endl; return 1; } @@ -1005,6 +1104,6 @@ // built output filename and path - TString outfile = "20191027_006_006-pulse.root"; - TString infile = "20190227_006_006-fit.root"; + TString outfile = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_pulse.root"; + TString infile = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_fit.root"; //if (!gSystem->AccessPathName(outfile)) @@ -1048,10 +1147,10 @@ // map file to use (get that from La Palma!) - const char *map = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL; + const char *map = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL; // ------------------------------------------------------ - Long_t max0 = 0; - Long_t max1 = 0; + Long_t max0 = 10000; + Long_t max1 = 10000; // ====================================================== @@ -1097,6 +1196,5 @@ MBadPixelsCam badpixels; badpixels.InitSize(64); - badpixels[7].SetUnsuitable(MBadPixelsPix::kUnsuitableRun); - badpixels[17].SetUnsuitable(MBadPixelsPix::kUnsuitableRun); + badpixels[43].SetUnsuitable(MBadPixelsPix::kUnsuitableRun); MH::SetPalette(); @@ -1219,7 +1317,14 @@ gStyle->SetOptStat(10); + TH2 *hsignals = h->GetSignals(); + TH2 *hoscillation = h->GetOscillation(); + + TH2 *hall = h->GetSignalAll(); + TH1 *h1 = h->GetSignal1(); TH1 *p1 = h->GetProf1(); + cout << "Entries:\t" << (p1->GetEntries()) << endl; + TH1 *h2 = h->GetSignal2(); TH1 *p2 = h->GetProf2(); @@ -1243,7 +1348,14 @@ TH1 *p8 = h->GetProf8(); - TF1 pulse("pulse", "[0]*(1-exp(-x*x/[1]))*exp(-x/[2])", 0, 300); - pulse.SetLineColor(kBlack); - pulse.SetParameters(0, 1./0.094, 1./0.05); + // This is the old fit + //TF1 pulse("pulse", "[0]*(1-exp(-x*x/[1]))*exp(-x/[2])", 0, 300); + //pulse.SetLineColor(kBlack); + //pulse.SetParameters(0, 1./0.094, 1./0.05); + + TF1 pulse("pulse", "[0]*(1-1/(1+exp((x-[1])/[2])))*exp(-(x-[1])/[3])", 0, 250); + pulse.SetLineColor(kRed); + pulse.SetParLimits(1, 1.0, 3); + pulse.SetParLimits(2, 0.5, 2.0); + pulse.SetParLimits(3, 50, 110); // ===== @@ -1259,16 +1371,39 @@ // ===== - + d->AddTab("Extracted Signals"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); + gPad->SetRightMargin(0.15); + hsignals->DrawCopy("colz"); + + d->AddTab("OscillationAmpl"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); + gPad->SetRightMargin(0.15); + hoscillation->DrawCopy("colz"); + + + d->AddTab("All Signals"); + gPad->SetGrid(); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); + gPad->SetRightMargin(0.15); + hall->DrawCopy("colz"); + + + // ===== d->AddTab("1"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); h1->DrawCopy("colz"); p1->DrawCopy("same"); - pulse.SetParameter(0, 1*1.3); - p1->Fit(&pulse, "N0", "", 1.5, 25); - TF1 *f1 = pulse.DrawCopy("same"); + //pulse.FixParameter(0, 1*1.5); + //p1->Fit(&pulse, "N0", "", 1, 70); + //TF1 *f1 = pulse.DrawCopy("same"); // ===== @@ -1276,6 +1411,6 @@ d->AddTab("2"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); h2->DrawCopy("colz"); @@ -1290,6 +1425,6 @@ d->AddTab("3"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); h3->DrawCopy("colz"); @@ -1304,6 +1439,6 @@ d->AddTab("4"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); h4->DrawCopy("colz"); @@ -1318,6 +1453,6 @@ d->AddTab("5"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); h5->DrawCopy("colz"); @@ -1332,6 +1467,6 @@ d->AddTab("6"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); h6->DrawCopy("colz"); @@ -1346,6 +1481,6 @@ d->AddTab("7"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); h7->DrawCopy("colz"); @@ -1360,6 +1495,6 @@ d->AddTab("8"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); h8->DrawCopy("colz"); @@ -1374,6 +1509,6 @@ d->AddTab("Pulse"); gPad->SetGrid(); - gPad->SetLeftMargin(0.07); - gPad->SetTopMargin(0.01); + gPad->SetLeftMargin(0.1); + gPad->SetTopMargin(0.1); gPad->SetRightMargin(0.15); p8->SetStats(kFALSE); @@ -1394,5 +1529,5 @@ p1->DrawCopy("same"); - f1->Draw("same"); + //f1->Draw("same"); f2->Draw("same"); f3->Draw("same"); @@ -1406,9 +1541,9 @@ + + d->SaveAs(outfile); + + TFile f(outfile, "UPDATE"); /* - d->SaveAs(outfile); - - TFile f(outfile, "UPDATE"); - MParameterI par("NumEvents"); par.SetVal(fill.GetNumExecutions()); Index: trunk/Mars/hawc/extract_singles.C =================================================================== --- trunk/Mars/hawc/extract_singles.C (revision 19950) +++ trunk/Mars/hawc/extract_singles.C (revision 19951) @@ -57,9 +57,8 @@ { Int_t fIntegrationWindow; - vector> fData; public: - MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(30) + MSingles(const char *name=NULL, const char *title=NULL) : fIntegrationWindow(40) { fName = name ? name : "MSingles"; @@ -88,5 +87,4 @@ ClassDef(MSingles, 1) }; - // Histogram class to extract the baseline class MHBaseline : public MH @@ -102,5 +100,5 @@ public: - MHBaseline() : fPedestal(0), fSkipStart(20), fSkipEnd(10) + MHBaseline() : fPedestal(0), fSkipStart(10), fSkipEnd(10) { fName = "MHBaseline"; @@ -108,7 +106,7 @@ // Setup the histogram fBaseline.SetName("Baseline"); - fBaseline.SetTitle("Median spectrum"); + fBaseline.SetTitle("Baseline Spectrum"); fBaseline.SetXTitle("Pixel [idx]"); - fBaseline.SetYTitle("Median baseline [mV]"); + fBaseline.SetYTitle("Amplitude [mV]"); fBaseline.SetDirectory(NULL); } @@ -282,4 +280,5 @@ TH2F fTime; TProfile2D fPulse; + TH2F fPulse2D; //plot the extracted singles in a 2D plot vs samples(time) UInt_t fNumEntries; @@ -296,5 +295,5 @@ // Setup histograms fSignal.SetName("Signal"); - fSignal.SetTitle("pulse integral spectrum"); + fSignal.SetTitle("Pulse Integral Spectrum"); fSignal.SetXTitle("pixel [SoftID]"); fSignal.SetYTitle("time [a.u.]"); @@ -302,5 +301,5 @@ fTime.SetName("Time"); - fTime.SetTitle("arival time spectrum"); + fTime.SetTitle("Arrival Time Spectrum"); fTime.SetXTitle("pixel [SoftID]"); fTime.SetYTitle("time [a.u.]"); @@ -308,8 +307,14 @@ fPulse.SetName("Pulse"); - fPulse.SetTitle("average pulse shape spectrum"); + fPulse.SetTitle("Average Pulse Shape Spectrum"); fPulse.SetXTitle("pixel [SoftID]"); fPulse.SetYTitle("time [a.u.]"); fPulse.SetDirectory(NULL); + + fPulse2D.SetName("Pulse2D"); + fPulse2D.SetTitle("Extracted Pulses"); + fPulse2D.SetXTitle("time [a.u.]"); + fPulse2D.SetYTitle("amplitude [a.u.]"); + fPulse2D.SetDirectory(NULL); } @@ -363,9 +368,19 @@ // Correct for DRS timing!! + + //Setup binning for the average time spectrum MBinning binst(roi, -0.5, roi-0.5); MH::SetBinning(fTime, binsx, binst); + //Setup binning for the average pulse spectrum MBinning binspy(2*roi, -roi-0.5, roi-0.5); + //MBinning binspy(1024, -522-0.5, 522-0.5); MH::SetBinning(fPulse, binsx, binspy); + + //Setup binning for the 2D plot amplitude vs time + MBinning bins_time(2*roi, -roi-0.5, roi-0.5); + //MBinning bins_time(1024, -522-0.5, 522-0.5); + MBinning bins_ampl(160, -30, 50); //N_bin is chosen such that 1 DAC = 0.5mV + MH::SetBinning(fPulse2D, bins_time, bins_ampl); return kTRUE; @@ -402,4 +417,7 @@ for (int s=0; sSetFrameBorderMode(0); fPulse.Draw("colz"); + + pad->cd(4); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fPulse2D.Draw("colz"); } @@ -465,4 +489,9 @@ gPad->SetFrameBorderMode(0); fPulse.DrawCopy("colz"); + + pad->cd(4); + gPad->SetBorderMode(0); + gPad->SetFrameBorderMode(0); + fPulse2D.Draw("colz"); } @@ -497,6 +526,6 @@ MExtractSingles() : fSingles(0), fPedestal(0), fEvent(0), fExtractionRange(64), - fNumAverage(10), fStartSkip(20), fEndSkip(10), - fIntegrationSize(3*10), fMaxSearchWindow(30) + fNumAverage(10), fStartSkip(10), fEndSkip(10), + fIntegrationSize(4*10), fMaxSearchWindow(20) { } @@ -538,5 +567,5 @@ const size_t navg = fNumAverage; - const float threshold = fThreshold; + const float threshold = fThreshold; const UInt_t start_skip = fStartSkip; const UInt_t end_skip = fEndSkip; @@ -595,4 +624,13 @@ if (k_max == i+5 || k_max == i + max_search_window-1) continue; + + //Make sure the pulse is isolated (area after the maximum is empty) + UInt_t k_falling = k_max+200; + for (UInt_t k_after = k_falling; k_after < k_max +250; k_after++) + { + if (val[k_after] > val[k_falling] && + val[k_after + fNumAverage/2] > val[k_falling]) + continue; + } //search for half maximum before maximum @@ -623,5 +661,5 @@ single.fTime = k_half_max; - // Crossing of the threshold is at 5 + // Crossing of the threshold is at 7 for (UInt_t j=k_half_max; jGetTime(); const TH2 *hpulse = h->GetPulse(); + const TH2F *ppulse2d = h->GetPulse2D(); d->AddTab("Time"); @@ -869,4 +909,7 @@ hp->Draw(); + d->AddTab("2DPulse"); + ppulse2d->DrawCopy("colz"); + d->SaveAs(outfile); Index: trunk/Mars/hawc/fit_singles.C =================================================================== --- trunk/Mars/hawc/fit_singles.C (revision 19950) +++ trunk/Mars/hawc/fit_singles.C (revision 19951) @@ -93,5 +93,5 @@ // Print function parameters -void PrintFunc(TF1 &f, float integration_window=30) +void PrintFunc(TF1 &f, float integration_window=40) { //cout << "--------------------" << endl; @@ -113,25 +113,23 @@ // The parameters for the function are the filenam from the gain extraction // and the output filename -int fit_singles(const char *filename = "20191027_006_006.root", - const char *outfile = "20190227_006_006-fit.root", +int fit_singles(const char *filename = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006.root", + const char *outfile = "/home/giangdo/Documents/Master/MA/pulse/data_single_pe/extracted/20191027_006_006_fit.root", bool fixednoise=false) { // Read the mapping between bias channels and hardware channels PixelMap pmap; - if (!pmap.Read("../hawc/FAMOUSmap171218.txt")) - { - cout << "FAMOUSmap171218.txt not found." << endl; + if (!pmap.Read("/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt")) + { + cout << "HAWCsEyemap181214.txt not found." << endl; return 1; } // It is more correct to fit the integral, but this is super - // slow, especially for 1440 pixel. To get that, one would have + // slow. To get that, one would have // to analytically integrate and fit this function. - // (Currently the integral is switched off for the 1440 individual - // spectra in all cases) bool fast = false; // Switch off using integral // Values which should be read from the file but not available atm - Int_t integration_window = 30; + Int_t integration_window = 40; // Map for which pixel shall be plotted and which not @@ -266,6 +264,6 @@ TH1F hCoeffR2 ("CoeffR2", "Coefficient R", 90, -1, 2); - // Histigram for the sum of all spectrums - TH1F hSum("Sum1", "Signal spectrum of all pixels", + // Histogram for the sum of all spectrums + TH1F hSum("Sum1", "Signal spectrum of all pixels (excl. blind pixels)", hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax()); hSum.SetXTitle("pulse integral [mV sample]"); @@ -275,5 +273,5 @@ // Histogram for the sum of all pixels excluding the ones with faulty fits - TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (incl TM)", + TH1F hSumClear1("SumC1", "Signal spectrum of all pixels (excl. faulty and blind pixels)", hsignal->GetNbinsY(), hsignal->GetYaxis()->GetXmin(), hsignal->GetYaxis()->GetXmax()); hSumClear1.SetXTitle("pulse integral [mV sample]"); @@ -303,6 +301,6 @@ hPulse.SetStats(false); - // Spektrum for the normalized individual spectra - TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (incl TM)", 120, 0.05, 12.05); + // Spectrum for the normalized individual spectra + TH1F hSumScale1("SumScale1", "Signal spectrum of all pixels (excl. blind pixels)", 120, 0.05, 12.05); hSumScale1.SetXTitle("pulse integral [pe]"); hSumScale1.SetYTitle("Rate"); @@ -332,4 +330,15 @@ if (usePixel[pixel]==0) continue; + + //exclude the pixels without cones + if (pixel == 61) + continue; + if (pixel == 62) + continue; + if (pixel == 63) + continue; + + + TH1D *hist = hsignal->ProjectionY("proj", pixel+1, pixel+1); @@ -434,4 +443,12 @@ if (usePixel[pixel]==0) continue; + + // Exclude pixels without cones + if (pixel == 61) + continue; + if (pixel == 62) + continue; + if (pixel == 63) + continue; //Projectipon of a certain Pixel out of Ampl.Spectrum Index: trunk/Mars/hawc/plot_trace.C =================================================================== --- trunk/Mars/hawc/plot_trace.C (revision 19950) +++ trunk/Mars/hawc/plot_trace.C (revision 19951) @@ -1,2 +1,44 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include "Getline.h" +#include "MH.h" +#include "MArrayI.h" +#include "MDirIter.h" +#include "MFillH.h" +#include "MEvtLoop.h" +#include "MCamEvent.h" +#include "MGeomApply.h" +#include "MTaskList.h" +#include "MParList.h" +#include "MContinue.h" +#include "MBinning.h" +#include "MDrsCalibApply.h" +#include "MDrsCalibration.h" +#include "MRawFitsRead.h" +#include "MStatusDisplay.h" +#include "MTaskInteractive.h" +#include "MPedestalSubtractedEvt.h" +#include "MGeomCamFAMOUS.h" +#include "MRawRunHeader.h" +#include "MPedestalCam.h" +#include "MPedestalPix.h" +#include "MParameters.h" + // ========================================================================== // ============ see plot_callisto function at the end of the file =========== @@ -21,5 +63,6 @@ // Create a histogram for the trace -TH1D hist("Trace", "Waveform", 1024, -0.25, 511.75); +//TH1D hist("Trace", "Waveform", 1024, -0.25, 511.75); //range = [-0.25,511.75]ns +TH1D hist("Trace", "Waveform (Pixel 5, Event 1)", 1024, -0.5, 1024-0.5); // All 'data members' that are required globally @@ -34,16 +77,20 @@ { fEvt = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt"); - if (!fEvt) - { - gfLog << err << "MPedestalSubtractedEvt not found... aborting." << endl; - return kALSE; - } + + //if (!fEvt) + //{ + // gfLog << err << "MPedestalSubtractedEvt not found... aborting." << endl; + // return kFALSE; + //} // Create a canvas and plot the histogram into it - c = new TCanvas; + c = new TCanvas(); hist.SetStats(kFALSE); hist.Draw(); - - TF1 constant("zero", "[0]", -0.25, 511.75); + hist.SetXTitle("Samples"); + hist.SetYTitle("Amplitude [a.u.]"); + + //TF1 constant("zero", "[0]", -0.25, 511.75); + TF1 constant("zero", "[0]", -0.5, 1024); constant.SetParameter(0, 0); constant.SetLineColor(kBlack); @@ -59,4 +106,5 @@ constant.DrawCopy("same"); + return kTRUE; } @@ -64,5 +112,7 @@ Int_t Process() { - int pixel_index = 0; + int pixel_index = 6; //plot trace only of this pixel + int skipstart = 5; + int skipend = 10; // Number of samples per pixel (usually 1024) @@ -73,19 +123,33 @@ hist.Reset(); for (int i=0; i-50) - hist.SetMinimum(-50); - if (hist.GetMaximum()<250) - hist.SetMaximum(250); + hist.SetMinimum(-50); + hist.SetMaximum(50); + + //// Set a reasonable range + //if (hist.GetMinimum()>-30) + // hist.SetMinimum(-20); + //if (hist.GetMaximum()<60) + // hist.SetMaximum(50); + + + + TLine *line = new TLine(skipstart, 50, skipstart, -50); + line->SetLineStyle(kSolid); + line->SetLineColor(kRed); + line->Draw(); + + TLine *line2 = new TLine(1024-skipend, 50, 1024-skipend, -50); + line2->SetLineStyle(kSolid); + line2->SetLineColor(kRed); + line2->Draw(); // Signal root to update the canvas/pad c->Modified(); c->Update(); + //c->SaveAs("/home/giangdo/Documents/Master/MA/pulse/data_single_pe/traces/pixel3/event_0.png"); + bool rc = HandleInput(); @@ -131,5 +195,5 @@ // mapping file (found in Mars/hawc) - const char *mmap = usemap ? "../hawc/FAMOUSmap171218.txt" : NULL; + const char *mmap = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL; // ====================================================== @@ -193,6 +257,6 @@ MEvtLoop loop(gSystem->BaseName(datafile)); loop.SetParList(&plist); - - if (!loop.Eventloop()) + + if (!loop.Eventloop(100)) return 4; Index: trunk/Mars/hawc/star.C =================================================================== --- trunk/Mars/hawc/star.C (revision 19950) +++ trunk/Mars/hawc/star.C (revision 19951) @@ -18,5 +18,5 @@ To run the macro from the command line (assuming you are in a directory - Mars/build where you have built your Mars environment) ou can do + Mars/build where you have built your Mars environment) you can do root ../hawc/star.C\(\"00000015.003_Y_MonteCarlo003_Events.root\"\) @@ -206,4 +206,5 @@ write.AddContainer("MPointingPos", "Events", kFALSE); write.AddContainer("MTime", "Events", kFALSE); + write.AddContainer("MRawBoardsFACT", "Events", kFALSE); write.AddContainer("MSrcPosCam", "Events", kFALSE); Index: trunk/Mars/hawc/synchron.C =================================================================== --- trunk/Mars/hawc/synchron.C (revision 19951) +++ trunk/Mars/hawc/synchron.C (revision 19951) @@ -0,0 +1,710 @@ +#include +#include +#include +#include +#include + +#include "TTree.h" +#include "TChain.h" +#include "TGraphErrors.h" +#include "TGraph.h" +#include "TCanvas.h" +#include "TFile.h" +#include "TF1.h" +#include "TH1.h" +#include "TH2.h" +#include "TH1F.h" +#include "TLegend.h" +#include + +#include "TROOT.h" +#include "TBrowser.h" + +#include "MTime.h" +//#include "MSignalCam.h" +#include "MRawEvtHeader.h" +//#include "MSoftwareTrigger.h" +#include "MRawBoardsFACT.h" + +using namespace std; +/* + C : a character string terminated by the 0 character + B : an 8 bit signed integer (Char_t) + b : an 8 bit unsigned integer (UChar_t) + S : a 16 bit signed integer (Short_t) + s : a 16 bit unsigned integer (UShort_t) + I : a 32 bit signed integer (Int_t) + i : a 32 bit unsigned integer (UInt_t) + F : a 32 bit floating point (Float_t) + D : a 64 bit floating point (Double_t) + L : a 64 bit signed integer (Long64_t) + l : a 64 bit unsigned integer (ULong64_t) + O : [the letter o, not a zero] a boolean (Bool_t) +*/ + +vector calibration(TString HEye_path) +{ + /* INPUT: HEye_path: Path to HAWC's Eye data (calibrated or cleaned) + OUTPUT: calibrated time in Unix_s + */ + // Vectors to store the Unix time after calibration + vector unix_s = {}; + + TFile *tHEye = new TFile(HEye_path); + + cout << "------------------------------------------------------------------------" << endl; + cout << "Do time calibration with file: " << HEye_path << endl; + + TTree *tEvents = (TTree*) tHEye->Get("Events"); + + MTime *mtime_utc = 0; // PC's timestamp in UNIX + MRawBoardsFACT *mtime_board = 0; // board's counter + + tEvents->SetBranchAddress("MTime.", &mtime_utc); + tEvents->SetBranchAddress("MRawBoardsFACT.", &mtime_board); + + Long64_t nEvents = tEvents->GetEntries(); + //cout << "Number of Events: " << nEvents << endl; + + Double_t Mjd; // time in Mjd + int Mjd_unix = 40587; // Mjd at 01/01/1970 + Double_t unix_timestamp; + + // For fitting and plotting + Double_t Unix_fit[nEvents], Board_fit[nEvents], UnixErr_fit[nEvents]; + + for (Long64_t i = 0; i < nEvents; i++) + { + tEvents->GetEntry(i); + + //mtime_utc->Print(); + Mjd = mtime_utc->GetMjd(); + UInt_t board_timestamp = mtime_board->GetFadTime(0); + + // calculate the PC's unix timestamp based on the PC UTC's time + Double_t unix_timestamp = (Mjd - Mjd_unix)*24*60*60; + + Unix_fit[i] = unix_timestamp; + Board_fit[i] = board_timestamp; + UnixErr_fit[i] = 0.001; // jitter of 1ms + + } + + // Use a linear fit to determine the offset and the conversion factor as the fit's slope + auto *board_unix = new TGraphErrors(nEvents, Board_fit, Unix_fit, 0, UnixErr_fit); + TF1 *fit_timestamp = new TF1("fit_timestamp", "[0] + [1]*x", Board_fit[0], Board_fit[nEvents-1]); + board_unix->Fit("fit_timestamp", "S"); + Double_t offset = fit_timestamp->GetParameter(0); + Double_t conv_timestamp = fit_timestamp->GetParameter(1); + + // Calculate the Unix timestamp based on board's counter and the calibration data + vector unix_time_s = {}; + for (Long64_t i = 0; i < nEvents; i++) + { + tEvents->GetEntry(i); + UInt_t board_timestamp = mtime_board->GetFadTime(0); + Double_t unix_time = board_timestamp * conv_timestamp + offset; + unix_time_s.push_back(unix_time); + } + + // Calculate the mean time difference between calibrated unix_timestamp and PC's timestamp + Double_t time_diff = 0.0; + for (Long64_t i = 0; i < nEvents; i++) + { + time_diff += TMath::Abs(Unix_fit[i] - unix_time_s[i]); + } + Double_t mean_time_diff = time_diff/nEvents; + printf("Mean time difference between PC's timestamps and calibrated counter's timestamp: %f s\n", mean_time_diff); + +/* + TCanvas *c = new TCanvas("time_calibration", "Board's time and PC's UTC time",1200,800); + c->SetGridx(); + c->SetGridy(); + + board_unix->SetTitle("FAD's timestamp vs. PC's UTC timestamp"); + board_unix->GetXaxis()->SetTitle("Board's counter"); + board_unix->GetYaxis()->SetTitle("Unix timestamp [s]"); + board_unix->SetMarkerStyle(8); + board_unix->SetMarkerColor(1); + gStyle->SetOptFit(0111); + + board_unix->Draw("AP"); + */ + + return unix_time_s; +} + +// structure of HAWC reconstructed data +struct rec_struct +{ Int_t status ; + Int_t version ; + Int_t eventID ; + Int_t runID ; + Int_t timeSliceID ; + Int_t trigger_flags; + Int_t event_flags ; + Int_t gtc_flags ; + Int_t gpsSec ; + Int_t gpsNanosec ; + Int_t nChTot ; + Int_t nChAvail ; + Int_t nHitTot ; + Int_t nHit ; + Int_t nHitSP10 ; + Int_t nHitSP20 ; + Int_t nTankTot ; + Int_t nTankAvail ; + Int_t nTankHitTot ; + Int_t nTankHit ; + Int_t windowHits ; + Int_t angleFitStatus; + Int_t planeNDOF ; + Int_t coreFitStatus; + Int_t CxPE40PMT ; + Int_t CxPE40XnCh ; + Double_t mPFnHits; + Double_t mPFnPlanes; + Double_t mPFp0nAssign; + Double_t mPFp1nAssign; + Double_t nHitMPF ; + Int_t coreFiduScale; + Double_t coreFiduScaleOR ; + Double_t coreLoc ; + Double_t zenithAngle ; + Double_t azimuthAngle ; + Double_t dec ; + Double_t ra ; + Double_t planeChi2 ; + Double_t coreX ; + Double_t coreY ; + Double_t logCoreAmplitude ; + Double_t coreFitUnc ; + Double_t logNNEnergy ; + Double_t fAnnulusCharge0; + Double_t fAnnulusCharge1; + Double_t fAnnulusCharge2; + Double_t fAnnulusCharge3; + Double_t fAnnulusCharge4; + Double_t fAnnulusCharge5; + Double_t fAnnulusCharge6; + Double_t fAnnulusCharge7; + Double_t fAnnulusCharge8; + Double_t protonlheEnergy; + Double_t protonlheLLH ; + Double_t gammalheEnergy ; + Double_t gammalheLLH ; + Int_t chargeFiduScale50 ; + Int_t chargeFiduScale70 ; + Int_t chargeFiduScale90 ; + Double_t logMaxPE; + Double_t logNPE ; + Int_t CxPE40 ; + Double_t CxPE40SPTime; + Double_t LDFAge ; + Double_t LDFAmp ; + Double_t LDFChi2 ; + Double_t logGP ; + Double_t mPFp0Weight; + Double_t mPFp0toangleFit; + Double_t mPFp1Weight ; + Double_t mPFp1toangleFit; + Double_t PINC ; + Double_t disMax ; + Double_t TankLHR ; + Double_t LHLatDistFitXmax ; + Double_t LHLatDistFitEnergy; + Double_t LHLatDistFitGoF ; + Double_t probaGBT ; + +}; + +struct tel_struct +{ + MTime *mTime = 0; + MRawEvtHeader *mRawEvtHeader = 0; + //MSignalCam *mSignalCam = 0; + //MRawEvtHeader *mRawEvtHeader = 0; + //MSoftwareTrigger *mSoftwareTrigger = 0; + //MRawBoardsFACT *mRawBoardsFACT = 0; + Double_t cal_unix; // calibrated unix timestamp based on board's counter + UInt_t DAQ_id; // event's DAQ ID (fDAQEvtNumber) + +}; + +struct sync_struct +{ + Double_t ts; + UInt_t id; +}; + +//Earlier timestamps first in vector: +bool sort_time (sync_struct str_1, sync_struct str_2) { return (str_1.ts unix_s = calibration(HEye_path); + + + //-------------------------------------------------------------------------------------------- + // Read in the ID list and timestamps and sort them + vector vts_HEye; + vector vts_Hawc; + sync_struct tsync; + + printf("Reading telescope data....\n"); + + Long64_t nHE_events = tHEye.GetEntries(); + cout << "Number of HawcEye's Events: " << nHE_events << endl; + + for(Int_t i = 0; i < nHE_events; i++) + { + tHEye.GetEntry(i); + tsync.id = i; + tsync.ts = unix_s[i]; + vts_HEye.push_back(tsync); + } + + printf("Sorting telescope data....\n"); + sort(vts_HEye.begin(), vts_HEye.end(),sort_time); + + + printf("Reading HAWC data....\n"); + + Long64_t nHawc_events = tHawc.GetEntries(); + cout << "Number of Hawc's Events: " << nHawc_events << endl; + for(Int_t i=0;i < nHawc_events; i++) + { + tHawc.GetEntry(i); + tsync.id = i; + tsync.ts = HAWC_toffset + rec.gpsSec + ((Double_t)rec.gpsNanosec)/1000000000.0; + vts_Hawc.push_back(tsync); + } + printf("Sorting HAWC data....\n"); + sort(vts_Hawc.begin(), vts_Hawc.end(),sort_time); + + // Find the best sync offset + printf("HAWC Data:\nFirst event: %f s Last event: %f s (Trigger: %ld, Rate: %f)\n\n", + vts_Hawc[0].ts, vts_Hawc[vts_Hawc.size()-1].ts, vts_Hawc.size(), + vts_Hawc.size() / (vts_Hawc[vts_Hawc.size()-1].ts - vts_Hawc[0].ts)); + + printf("HEye Data:\nFirst event: %f s Last event: %f s (Trigger: %ld, Rate: %f)\n\n", + vts_HEye[0].ts, vts_HEye[vts_HEye.size()-1].ts, vts_HEye.size(), + vts_HEye.size()/(vts_HEye[vts_HEye.size()-1].ts-vts_HEye[0].ts)); + + + if (vts_Hawc[0].ts>vts_HEye[nHE_events-1].ts) + { + printf("NO OVERLAP: HAWC data later, choose later HEye file!\n"); + return 1; + } + if (vts_Hawc[nHawc_events-1].ts first HE's time - Start_Shift + Int_t last_shift = 0; // index of last Hawc's event with time < last HE's time - Stop_Shift + Int_t first = 0; // index of first Hawc's event with trigger time > first HE's time + Int_t last = 0 ; // index of last Hawc's event with trigger time < last HE's time + + for (Long64_t i = 0; i < nHawc_events; i++) + { + if (vts_Hawc[i].ts > vts_HEye[0].ts) + { + first = i; + } + if (vts_Hawc[i].ts > vts_HEye[0].ts - Start_Shift) + { + first_shift = i; + break; + } + } + + for(Long64_t i = nHawc_events-1; i > 0; i--) + { + if (vts_Hawc[i].ts < vts_HEye[nHE_events-1].ts) + { + last = i; + } + if (vts_Hawc[i].ts < vts_HEye[nHE_events-1].ts - Stop_Shift) + { + last_shift = i; + break; + } + } + + printf("HAWC Data (reduced):\nFirst event: %f s Last event: %f s\n\n", + vts_Hawc[first_shift].ts, vts_Hawc[last_shift].ts); + + printf("HEye Data:\nFirst event: %f s Last event: %f s\n\n", + vts_HEye[0].ts, vts_HEye[vts_HEye.size()-1].ts); + + //---------------------- Define the Output TTree ---------------------- + TFile fOut(Outpath, "RECREATE"); + + TTree tOut1("Hawc_Data", "Hawc_Data"); // Hawc's data + TTree tOut2("HE_Data", "HE_Data"); // Telescope's data + TTree tOut3("timediff_tot", "timediff_tot"); // Synchronisation's info + TTree tOut4("timediff_pass", "timediff_pass"); + TTree tOut5("SyncInfo", "SyncInfo"); + + // Hawc Data + tOut1.Branch("rec.status" , &rec.status) ; + tOut1.Branch("rec.version", &rec.version) ; + tOut1.Branch("rec.eventID", &rec.eventID) ; + tOut1.Branch("rec.runID", &rec.runID) ; + tOut1.Branch("rec.timeSliceID", &rec.timeSliceID); + tOut1.Branch("rec.trigger_flags", &rec.trigger_flags); + tOut1.Branch("rec.event_flags",&rec.event_flags) ; + tOut1.Branch("rec.gtc_flags",&rec.gtc_flags) ; + tOut1.Branch("rec.gpsSec",&rec.gpsSec) ; + tOut1.Branch("rec.gpsNanosec",&rec.gpsNanosec) ; + tOut1.Branch("rec.nChTot",&rec.nChTot) ; + tOut1.Branch("rec.nChAvail",&rec.nChAvail) ; + tOut1.Branch("rec.nHitTot",&rec.nHitTot) ; + tOut1.Branch("rec.nHit",&rec.nHit) ; + tOut1.Branch("rec.nHitSP10",&rec.nHitSP10) ; + tOut1.Branch("rec.nHitSP20",&rec.nHitSP20) ; + tOut1.Branch("rec.nTankTot",&rec.nTankTot) ; + tOut1.Branch("rec.nTankAvail",&rec.nTankAvail) ; + tOut1.Branch("rec.nTankHitTot",&rec.nTankHitTot) ; + tOut1.Branch("rec.nTankHit",&rec.nTankHit) ; + tOut1.Branch("rec.windowHits",&rec.windowHits) ; + tOut1.Branch("rec.angleFitStatus",&rec.angleFitStatus); + tOut1.Branch("rec.planeNDOF",&rec.planeNDOF) ; + tOut1.Branch("rec.coreFitStatus",&rec.coreFitStatus); + tOut1.Branch("rec.CxPE40PMT",&rec.CxPE40PMT) ; + tOut1.Branch("rec.CxPE40XnCh",&rec.CxPE40XnCh); + tOut1.Branch("rec.mPFnHits",&rec.mPFnHits); + tOut1.Branch("rec.mPFnPlanes",&rec.mPFnPlanes); + tOut1.Branch("rec.mPFp0nAssign",&rec.mPFp0nAssign); + tOut1.Branch("rec.mPFp1nAssign",&rec.mPFp1nAssign); + tOut1.Branch("rec.nHitMPF",&rec.nHitMPF) ; + tOut1.Branch("rec.coreFiduScale",&rec.coreFiduScale); + tOut1.Branch("rec.coreFiduScaleOR",&rec.coreFiduScaleOR) ; + tOut1.Branch("rec.coreLoc",&rec.coreLoc) ; + tOut1.Branch("rec.zenithAngle",&rec.zenithAngle) ; + tOut1.Branch("rec.azimuthAngle",&rec.azimuthAngle) ; + tOut1.Branch("rec.dec",&rec.dec) ; + tOut1.Branch("rec.ra",&rec.ra) ; + tOut1.Branch("rec.planeChi2",&rec.planeChi2) ; + tOut1.Branch("rec.coreX",&rec.coreX) ; + tOut1.Branch("rec.coreY",&rec.coreY) ; + tOut1.Branch("rec.logCoreAmplitude",&rec.logCoreAmplitude) ; + tOut1.Branch("rec.coreFitUnc",&rec.coreFitUnc) ; + tOut1.Branch("rec.logNNEnergy",&rec.logNNEnergy) ; + tOut1.Branch("rec.fAnnulusCharge0",&rec.fAnnulusCharge0); + tOut1.Branch("rec.fAnnulusCharge1",&rec.fAnnulusCharge1); + tOut1.Branch("rec.fAnnulusCharge2",&rec.fAnnulusCharge2); + tOut1.Branch("rec.fAnnulusCharge3",&rec.fAnnulusCharge3); + tOut1.Branch("rec.fAnnulusCharge4",&rec.fAnnulusCharge4); + tOut1.Branch("rec.fAnnulusCharge5",&rec.fAnnulusCharge5); + tOut1.Branch("rec.fAnnulusCharge6",&rec.fAnnulusCharge6); + tOut1.Branch("rec.fAnnulusCharge7",&rec.fAnnulusCharge7); + tOut1.Branch("rec.fAnnulusCharge8",&rec.fAnnulusCharge8); + tOut1.Branch("rec.protonlheEnergy",&rec.protonlheEnergy); + tOut1.Branch("rec.protonlheLLH",&rec.protonlheLLH) ; + tOut1.Branch("rec.gammalheEnergy",&rec.gammalheEnergy) ; + tOut1.Branch("rec.gammalheLLH",&rec.gammalheLLH) ; + tOut1.Branch("rec.chargeFiduScale50",&rec.chargeFiduScale50) ; + tOut1.Branch("rec.chargeFiduScale70",&rec.chargeFiduScale70) ; + tOut1.Branch("rec.chargeFiduScale90",&rec.chargeFiduScale90) ; + tOut1.Branch("rec.logMaxPE",&rec.logMaxPE); + tOut1.Branch("rec.logNPE",&rec.logNPE) ; + tOut1.Branch("rec.CxPE40",&rec.CxPE40) ; + tOut1.Branch("rec.CxPE40SPTime",&rec.CxPE40SPTime); + tOut1.Branch("rec.LDFAge",&rec.LDFAge) ; + tOut1.Branch("rec.LDFAmp",&rec.LDFAmp) ; + tOut1.Branch("rec.LDFChi2",&rec.LDFChi2) ; + tOut1.Branch("rec.logGP",&rec.logGP) ; + tOut1.Branch("rec.mPFp0Weight",&rec.mPFp0Weight); + tOut1.Branch("rec.mPFp0toangleFit",&rec.mPFp0toangleFit); + tOut1.Branch("rec.mPFp1Weight",&rec.mPFp1Weight) ; + tOut1.Branch("rec.mPFp1toangleFit",&rec.mPFp1toangleFit); + tOut1.Branch("rec.PINC",&rec.PINC) ; + tOut1.Branch("rec.disMax",&rec.disMax) ; + tOut1.Branch("rec.TankLHR",&rec.TankLHR) ; + tOut1.Branch("rec.LHLatDistFitXmax",&rec.LHLatDistFitXmax) ; + tOut1.Branch("rec.LHLatDistFitEnergy",&rec.LHLatDistFitEnergy); + tOut1.Branch("rec.LHLatDistFitGoF",&rec.LHLatDistFitGoF) ; + tOut1.Branch("rec.probaGBT",&rec.probaGBT) ; + + // Telescope's Data + tOut2.Branch("cal_unix", &tel.cal_unix); + tOut2.Branch("DAQ_id", &tel.DAQ_id); + + // Synchronisation Info + Double_t min_miss = 100; // [s] + Double_t best_shift = 0.0; // best time shift + Double_t timediff_tot; // time difference after best shift - for all events + Double_t timediff_passed; // time difference after best shift - for events within TimeDiff_Cut + + tOut3.Branch("timediff_tot", &timediff_tot); // time diff between neighbor events after shift for all events + tOut4.Branch("timediff_passed", &timediff_passed); // time diff between neighbors after shift for sync events + tOut5.Branch("min_miss", &min_miss); // mean time diff with best shift + tOut5.Branch("best_shift", &best_shift); // best time shift + + //---------------------------------------------------------------------------- + + Double_t total_miss = 0.0; // total time difference sum( |t_Hawc(i) - t_HE(i)|) [s] + Double_t mean_miss = 0.0; // total_miss/n_Hawc + + Double_t diff_p = 0.0; + Double_t diff_n = 0.0; + Double_t diff_min = 0.0; + Int_t temp_ind = 0; + + Int_t last_index = 0; + + + // For the plot Mean Time Difference vs. Time Shift + vector gr_tShift, gr_tDiff; + Int_t gr_n = (Stop_Shift-Start_Shift)/shift_step; + + printf("Doing shift match....\n"); + + // Search for best shift with min. mean timem difference + for (Double_t shift = Start_Shift; shift < Stop_Shift; shift += shift_step) + { + for (Int_t i = first_shift; i < last_shift; i++) + { + + //1. Find the nearest neighbor in vts_HEye to vts_Hawc[i]+shift: + while (1) + { + if(vts_Hawc[i].ts + shift > vts_HEye[last_index + 1].ts) + last_index++; + else + break; + } + //2. calculate the difference + diff_n = vts_Hawc[i].ts+shift - vts_HEye[last_index].ts; + diff_p = vts_HEye[last_index+1].ts - vts_Hawc[i].ts-shift; + diff_min = diff_n > diff_p ? diff_p : diff_n; + + //3. Add to total_miss + total_miss += diff_min; + } + mean_miss = total_miss/(last_shift-first_shift); + + gr_tShift.push_back(shift); + gr_tDiff.push_back(mean_miss); + + if (min_miss > mean_miss) + { + min_miss = mean_miss; + best_shift = shift; + } + total_miss = 0.0; + last_index = 0; + + } + tOut5.Fill(); + + printf("->Minimum sync miss: %f @ a shift of %fs\n", min_miss, best_shift); + + + TH1F timediff_neighbor("timediff_neighbor","timediff_neighbor",1000,-0.01,0.01); + + // Shift all events by the found best shift + for(Int_t i = first; i < last; i++) + { + //1. Find the nearest neighbor in vts_HEye to vts_Hawc[i]+best_shift: + while(1) + { + if(vts_Hawc[i].ts + best_shift > vts_HEye[last_index+1].ts) + last_index++; + else + break; + } + + //2. calculate the difference + diff_n = vts_Hawc[i].ts + best_shift - vts_HEye[last_index].ts; + diff_p = vts_HEye[last_index+1].ts - vts_Hawc[i].ts-best_shift; + + if (diff_p > diff_n) + { + diff_min = diff_n; + temp_ind = last_index; + } + else + { + diff_min = diff_p; + temp_ind = last_index+1; + } + + timediff_tot = diff_min; + + if (diff_min < TimeDiff_Cut) + { + tHEye.GetEntry(vts_HEye[temp_ind].id); + tHawc.GetEntry(vts_Hawc[i].id); + + tel.cal_unix = vts_HEye[temp_ind].ts; + tel.DAQ_id = tel.mRawEvtHeader->GetDAQEvtNumber(); + + timediff_passed = diff_min; + + tOut1.Fill(); + tOut2.Fill(); + tOut4.Fill(); + + } + timediff_neighbor.Fill(diff_min); //Fill historgramm with diff + tOut3.Fill(); + } + + //---------------------------Plots and Output files--------------------------- + + tOut1.Write(); + tOut2.Write(); + tOut3.Write(); + tOut4.Write(); + tOut5.Write(); + fOut.Close(); + + TCanvas *c1 = new TCanvas("timediff_shift", "Synchronisation", 1200,800); + c1->SetGridx(); + c1->SetGridy(); + + TGraph *gr_shift = new TGraph(gr_n, &gr_tShift[0], &gr_tDiff[0]); + gr_shift->SetTitle("Synchronisation"); + gr_shift->GetXaxis()->SetTitle("Time Shift [s]"); + gr_shift->GetYaxis()->SetTitle("Mean time difference [s]"); + gr_shift->Draw(); + + TCanvas *c2 = new TCanvas("timediff_neighbor", "timediff_neighbor", 1200, 800); + c2->SetGridx(); + c2->SetGridy(); + timediff_neighbor.GetXaxis()->SetTitle("Time Difference [s]"); + timediff_neighbor.GetYaxis()->SetTitle("Counts"); + timediff_neighbor.DrawClone(); + + // Write out files + + TFile fOutPlots(Outpath+"_syncplots.root","RECREATE"); + c1->Write(); + c2->Write(); + fOutPlots.Close(); + + return 0; + +}