Fix: cutEnable, tof histogram draw, fill event hist only if there good tracks; Feat: draw option

This commit is contained in:
liuyihui 2024-11-25 13:50:15 +08:00
parent fc40af7463
commit 8fae28b96a
3 changed files with 39 additions and 30 deletions

View File

@ -143,7 +143,7 @@ public:
void readTreeData(vector<string> files); void readTreeData(vector<string> files);
void FillHistograms(); void FillHistograms();
void postOperations(); void postOperations();
void saveHistograms(TString); void saveHistograms(TString filename);
void drawHistograms(); void drawHistograms();
private: private:

24
main.cc
View File

@ -13,22 +13,25 @@ int main(int argc, char **argv) {
string cfgfile; string cfgfile;
string datafile; string datafile;
string outputfile = "bluetAnalyzer.root"; string outputfile = "bluetAnalyzer.root";
bool draw = false;
vector<string> ranges; vector<string> ranges;
vector<string> runfiles; vector<string> runfiles;
Long_t rl, rr; Long_t rl, rr;
auto cfg = value("config file", cfgfile); auto cfg = value("config file", cfgfile);
auto data = (((required("-f", "--file").doc("data file path") & auto data = ((required("-f", "--file").doc("data file path") &
value("data file", datafile)) | value("data file", datafile)) |
(required("-r", "--range").doc("data file number range") & (required("-r", "--range").doc("data file number range") &
values("number1 number2", ranges))), values("number1 number2", ranges)));
option("-o", "--output") auto output = (option("-o", "--output")
.doc("output file path, default: bluetAnalyzer.root") & .doc("output file path, default: bluetAnalyzer.root") &
value("output file", outputfile)); value("output file", outputfile),
option("-d", "--draw").doc("draw histograms").set(draw, true));
auto cmd = (cfg & data & output);
auto fmt = doc_formatting{}.first_column(4); auto fmt = doc_formatting{}.first_column(4);
if (!parse(argc, const_cast<char **>(argv), cfg & data)) { if (!parse(argc, const_cast<char **>(argv), cmd)) {
cerr << make_man_page(cfg & data, "BluetAnalyzer", fmt) cerr << make_man_page(cmd, "BluetAnalyzer", fmt)
.prepend_section("DESCRIPTION", .prepend_section("DESCRIPTION",
" Program for MTPC Data Analysis.") " Program for MTPC Data Analysis.")
<< endl; << endl;
@ -60,7 +63,8 @@ int main(int argc, char **argv) {
bluet->FillHistograms(); bluet->FillHistograms();
bluet->postOperations(); bluet->postOperations();
bluet->saveHistograms(outputfile); bluet->saveHistograms(outputfile);
bluet->drawHistograms(); if (draw)
bluet->drawHistograms();
return 0; return 0;
} }

View File

@ -144,7 +144,7 @@ void BluetAnalyzer::readCutParameters(string cfgfile) {
dv0 = bluet::XMLToQuantity(gas.child("DriftVelocity")) * 1e4; dv0 = bluet::XMLToQuantity(gas.child("DriftVelocity")) * 1e4;
pugi::xml_node cut = root.child("CutParameters"); pugi::xml_node cut = root.child("CutParameters");
cutEnable = bluet::XMLToQuantity(cut.child("CutEnable")); cutEnable = bluet::XMLToQuantity(cut.child("cutEnable"));
EcutMin = bluet::XMLToQuantity(cut.child("ECutMin")); EcutMin = bluet::XMLToQuantity(cut.child("ECutMin"));
EcutMax = bluet::XMLToQuantity(cut.child("ECutMax")); EcutMax = bluet::XMLToQuantity(cut.child("ECutMax"));
tofCutMin = bluet::XMLToQuantity(cut.child("TOFCutMin")); tofCutMin = bluet::XMLToQuantity(cut.child("TOFCutMin"));
@ -165,7 +165,7 @@ void BluetAnalyzer::readCutParameters(string cfgfile) {
lengthCutMax = bluet::XMLToQuantity(cut.child("LengthCutMax")); lengthCutMax = bluet::XMLToQuantity(cut.child("LengthCutMax"));
cosCutMin = bluet::XMLToQuantity(cut.child("cosCutMin")); cosCutMin = bluet::XMLToQuantity(cut.child("cosCutMin"));
CenterPos[0] = -3.03631; CenterPos[0] = 0;
CenterPos[1] = 0.634384; CenterPos[1] = 0.634384;
sliceW = 2.5; // [mm] sliceW = 2.5; // [mm]
rCut = bluet::XMLToQuantity(cut.child("radiusCutEnable")); rCut = bluet::XMLToQuantity(cut.child("radiusCutEnable"));
@ -474,26 +474,13 @@ void BluetAnalyzer::FillHistograms() {
double ts = fEvent->TimeStamp; double ts = fEvent->TimeStamp;
double tof = fEvent->TOF0; double tof = fEvent->TOF0;
double tcath = fEvent->CathodeTime; double tcath = fEvent->CathodeTime;
hts->Fill(ts, tof);
// 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);
htt->Fill(hit->HitT * sampleT / 1e3, hit->HitTau);
}
bool goodTrack = false;
for (size_t itrk = 0; itrk < fEvent->Tracks.size(); ++itrk) { for (size_t itrk = 0; itrk < fEvent->Tracks.size(); ++itrk) {
if (!isGoodTrack(&(fEvent->Tracks[itrk]))) if (!isGoodTrack(&(fEvent->Tracks[itrk])))
continue; continue;
goodTrack = true;
double sumE = fEvent->Tracks[itrk].TrSumEdep; double sumE = fEvent->Tracks[itrk].TrSumEdep;
double length = fEvent->Tracks[itrk].TrRawLength; double length = fEvent->Tracks[itrk].TrRawLength;
double traTmax = fEvent->Tracks[itrk].Tmax; double traTmax = fEvent->Tracks[itrk].Tmax;
@ -540,6 +527,24 @@ void BluetAnalyzer::FillHistograms() {
htofE->Fill(tof, sumE); htofE->Fill(tof, sumE);
} }
} }
if (goodTrack) {
hts->Fill(ts, tof);
// 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);
htt->Fill(hit->HitT * sampleT / 1e3, hit->HitTau);
}
}
} }
printf("\n===> Finish reading data ...\n"); printf("\n===> Finish reading data ...\n");
} }
@ -701,9 +706,9 @@ void BluetAnalyzer::drawTOF() {
pad->cd(); pad->cd();
gPad->SetLogx(); gPad->SetLogx();
gPad->SetLogy(); gPad->SetLogy();
htofr0->Draw("HIST"); htof->Draw("HIST");
htofr0->Draw("same HIST");
htofr1->Draw("same HIST"); htofr1->Draw("same HIST");
htof->Draw("same HIST");
TLegend *leg = new TLegend(0.45, 0.75, 0.8, 0.9); TLegend *leg = new TLegend(0.45, 0.75, 0.8, 0.9);
leg->AddEntry(htofr0, Form("r < %.1f mm", rCutMin)); leg->AddEntry(htofr0, Form("r < %.1f mm", rCutMin));
leg->AddEntry(htofr1, Form("%.1f < r < %.1f mm", rCutMin, rCutMax)); leg->AddEntry(htofr1, Form("%.1f < r < %.1f mm", rCutMin, rCutMax));