Feat: save histograms; calculate neutron energy

include/BluetAnalyzer.hh

@@ -25,8 +25,8 @@
 
 #include "BluetDataModel.hh"
 
-const int sizeX = 900;
+const int sizeX = 1200;
-const int sizeY = 900;
+const int sizeY = 1200;
 const double rBeam = 35. / 2;
 const double rTarget = 73. / 2;
 const double rSubstrate = 89. / 2;
@@ -59,6 +59,8 @@ private:
   double sidelength;
   int Nlayer;
   double driftL;
+  double kickdT;
+  double flighL;
   double sampleT;
 
   double EcutMin;
@@ -128,6 +130,8 @@ private:
   TH2F *han;
   TH2F *hta;
   TH2F *htt;
+  TH2F *hte;
+  TH1F *hec;
 
 public:
   void readCutParameters(string cfgfile);
@@ -138,6 +142,7 @@ public:
   void readTreeData(vector<string> files);
   void FillHistograms();
   void postOperations();
+  void saveHistograms();
   void drawHistograms();
 
 private:
@@ -148,4 +153,5 @@ private:
   void draw1DFeatures();
   void draw2DFeatures();
   void drawPadFeatures();
+  void drawExcitationCurve();
 };

main.cc

@@ -22,6 +22,7 @@ int main(int argc, char **argv) {
   bluet->readTreeData(runfiles);
   bluet->FillHistograms();
   bluet->postOperations();
+  bluet->saveHistograms();
   bluet->drawHistograms();
 
   return 0;

src/BluetAnalyzer.cc

@@ -43,6 +43,8 @@ double dgaus(double *x, double *par) {
 static double dv0;
 static int nInner = 0;
 static int nRing = 0;
+static const double vc = 3 * 1e8;     // [m/s]
+static const double mn = 939.5654133; // [MeV/c^2]
 
 void fitdgaus(TH1F *hist, TF1 *func, double *&par, double m0) {
   double mean = m0 >= 0 ? m0 : hist->GetMean();
@@ -61,8 +63,8 @@ void fitdgaus(TH1F *hist, TF1 *func, double *&par, double m0) {
   double amp2 = pars[3];
   double sigma2 = pars[4];
   sigma = sqrt(sigma1 * sigma1 + sigma2 * sigma2);
-  xl = mean - 3. * sigma;
+  xl = max(mean - 3. * sigma, m0 * 0.6);
-  xr = mean + 3. * sigma;
+  xr = min(mean + 3. * sigma, m0 * 1.4);
   par = new double[7]{mean, amp1, sigma1, amp2, sigma2, xl, xr};
 }
 
@@ -132,9 +134,11 @@ void BluetAnalyzer::readCutParameters(string cfgfile) {
   sidelength = bluet::XMLToQuantity(detector.child("PadSideLength")); // mm
   Nlayer = bluet::XMLToQuantity(detector.child("PadArrayNlayer"));
   driftL = bluet::XMLToQuantity(detector.child("DriftLength"));
+  flighL = bluet::XMLToQuantity(detector.child("FlightLength"));
 
   pugi::xml_node electronics = root.child("ElectronicsParameters");
   sampleT = bluet::XMLToQuantity(electronics.child("WaveformSamplingPeriod"));
+  kickdT = bluet::XMLToQuantity(electronics.child("KickerDelay"));
 
   pugi::xml_node gas = root.child("GasParameters");
   dv0 = bluet::XMLToQuantity(gas.child("DriftVelocity")) * 1e4;
@@ -316,13 +320,13 @@ void BluetAnalyzer::defineHistograms() {
   htpad->SetMarkerColor(kBlack);
   htpad->SetMarkerStyle(20);
 
-  hv0 = new TH1F("hv0", "hv0;Drift Veloctiy [mm/#mus];Count", 300, 0., 300.);
+  hv0 = new TH1F("hv0", "hv0;Drift Veloctiy [mm/#mus];Count", 150, 0., 150.);
   hv0->GetXaxis()->CenterTitle();
   hv0->GetYaxis()->CenterTitle();
   hv0->SetMarkerColor(kBlack);
   hv0->SetMarkerStyle(20);
 
-  hv1 = new TH1F("hv1", "hv1;Drift Veloctiy [mm/#mus];Count", 300, 0., 300.);
+  hv1 = new TH1F("hv1", "hv1;Drift Veloctiy [mm/#mus];Count", 150, 0., 150.);
   hv1->GetXaxis()->CenterTitle();
   hv1->GetYaxis()->CenterTitle();
   hv1->SetMarkerColor(kBlack);
@@ -380,17 +384,28 @@ void BluetAnalyzer::defineHistograms() {
   han->GetYaxis()->CenterTitle();
   han->GetYaxis()->SetTitleOffset(1.5);
 
-  hta =
+  hta = new TH2F("hta", "hta;#tau [us];Amplitude [ADC]", 100, 0, 5, 450, 0,
-      new TH2F("hta", "hta;#tau [us];Amplitude [ADC]", 100, 0, 5, 300, 0,
+                 4.5e3);
-      3e3);
   hta->GetXaxis()->CenterTitle();
   hta->GetYaxis()->CenterTitle();
   hta->GetYaxis()->SetTitleOffset(1.5);
 
-  htt = new TH2F("htt", "htt;T_0 [us];#tau [us]", 100, 0, 10, 100, 0, 5);
+  htt = new TH2F("htt", "htt;T_{0} [us];#tau [us]", 150, 0, 15, 100, 0, 5);
   htt->GetXaxis()->CenterTitle();
   htt->GetYaxis()->CenterTitle();
   htt->GetYaxis()->SetTitleOffset(1.5);
+
+  const double estep = (Log10(1e6) - Log10(1e-1)) / Nbin;
+  double epos[Nbin + 1];
+  for (int ie = 0; ie < (Nbin + 1); ie++)
+    epos[ie] = 1e-1 * TMath::Power(10., ie * estep);
+  hte = new TH2F("hte", "hte;TOF [ns];Neutron Energy [ev]", Nbin, tpos, Nbin,
+                 epos);
+  hte->GetXaxis()->CenterTitle();
+  hte->GetYaxis()->CenterTitle();
+  hte->GetYaxis()->SetMaxDigits(2);
+
+  hec = new TH1F("hec", "hec;Neutron Energy [ev];Count", Nbin, epos);
 }
 
 bool BluetAnalyzer::isGoodEvent(BluetEvent *event) {
@@ -451,6 +466,14 @@ void BluetAnalyzer::FillHistograms() {
     if (!isGoodEvent(fEvent))
       continue;
 
+    // convert time of flight to neutron energy using relativistic formula
+    double tof1 = fEvent->TOF0 - kickdT;
+    double vn = flighL / tof1 * 1e9;
+    double gamma = 1. / Sqrt(1. - vn * vn / (vc * vc));
+    double en = (gamma - 1) * mn * 1e6; // [MeV]->[eV]
+    hec->Fill(en);
+    hte->Fill(tof1, en);
+
     for (size_t ihit = 0; ihit < fEvent->Hits.size(); ++ihit) {
       BluetHit *hit = &(fEvent->Hits[ihit]);
       hta->Fill(hit->HitTau, hit->HitAmp);
@@ -918,29 +941,91 @@ void BluetAnalyzer::drawPadFeatures() {
   pad->cd();
   adjustAxisRange(hta);
   hta->Draw("colz");
-  c0->cd();
 
+  c0->cd();
   pad = new TPad("tn", "", 0.52, 0.02, 0.98, 0.48);
   pad->SetRightMargin(0.15);
   pad->Draw();
   pad->cd();
   adjustAxisRange(htt);
   htt->Draw("colz");
-  c0->cd();
 
+  c0->cd();
   pad = new TPad("amp", "", 0.02, 0.52, 0.48, 0.98);
   pad->Draw();
   pad->cd();
-  c0->cd();
 
+  c0->cd();
   pad = new TPad("an", "", 0.52, 0.52, 0.98, 0.98);
   pad->Draw();
   pad->cd();
-  c0->cd();
 
   c0->Update();
 }
 
+void BluetAnalyzer::drawExcitationCurve() {
+  TGCompositeFrame *tab = mTab->AddTab("Excitation Curve");
+  TRootEmbeddedCanvas *ec =
+      new TRootEmbeddedCanvas("Canva 2d", tab, sizeX, sizeY);
+  tab->AddFrame(ec, new TGLayoutHints(kLHintsExpandX | kLHintsExpandY));
+  TCanvas *c0 = ec->GetCanvas();
+  TPad *pad;
+
+  c0->cd();
+  pad = new TPad("te", "", 0.02, 0.02, 0.98, 0.48);
+  pad->SetRightMargin(0.15);
+  pad->Draw();
+  pad->SetLogx();
+  pad->SetLogy();
+  pad->cd();
+  hte->Draw("colz");
+
+  c0->cd();
+  pad = new TPad("ec", "", 0.02, 0.52, 0.98, 0.98);
+  pad->SetRightMargin(0.15);
+  pad->Draw();
+  pad->SetLogx();
+  pad->SetLogy();
+  pad->cd();
+  hec->Draw("colz");
+
+  c0->Update();
+}
+
+void BluetAnalyzer::saveHistograms() {
+  TFile *fout = new TFile("bluetAnalyzer.root", "RECREATE");
+  fout->cd();
+  hprofile->Write();
+  hslicex->Write();
+  hslicey->Write();
+  hbeamr->Write();
+  hrinte->Write();
+  htofE->Write();
+  htof->Write();
+  htofr0->Write();
+  htofr1->Write();
+  hpull->Write();
+  hratio->Write();
+  hLengthE->Write();
+  htcath->Write();
+  htpad->Write();
+  hv0->Write();
+  hv1->Write();
+  hcos->Write();
+  hphi->Write();
+  hdeltacos1->Write();
+  hdeltacos2->Write();
+  hts->Write();
+  htn->Write();
+  hamp->Write();
+  han->Write();
+  hta->Write();
+  htt->Write();
+  hte->Write();
+  hec->Write();
+  fout->Close();
+}
+
 void BluetAnalyzer::drawHistograms() {
   gROOT->SetStyle("ATLAS");
   gStyle->SetOptFit(1);
@@ -954,6 +1039,7 @@ void BluetAnalyzer::drawHistograms() {
   draw1DFeatures();
   draw2DFeatures();
   drawPadFeatures();
+  drawExcitationCurve();
 
   mFrame->SetWindowName("Bluet Analyzer");
   mFrame->MapSubwindows();