-
Administrator authored
Apply code formatting to all source/header files and root macros.
Administrator authoredApply code formatting to all source/header files and root macros.
standalone_event_build.C 8.77 KiB
/** @file MCBM PSD Energy profile
** @author Nikolay Karpushkin <karpushkin@inr.ru>
** @date 15.12.2019
** ROOT macro to visualize psd energy in psd digi
*/
void standalone_event_build(TString inFileName) {
TString srcDir = gSystem->Getenv("VMCWORKDIR");
TFile* inFile = TFile::Open(inFileName.Data(), "READONLY");
TTree* inTree = (TTree*) inFile->Get("cbmsim");
if (inTree == nullptr) {
cout << "tree not found" << endl;
return;
}
std::vector<CbmPsdDigi>* fPsdDigis = new std::vector<CbmPsdDigi>;
inTree->SetBranchAddress("PsdDigi", &fPsdDigis);
if (fPsdDigis == nullptr) {
cout << "no psd digis vector in tree" << endl;
return;
}
const Float_t event_time_gate = 200.; //in ns
const Int_t total_channels = 9;
const Int_t total_cells = 9;
Double_t ch_energy[total_channels];
for (UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++)
ch_energy[ch_iter] = 0.;
Double_t ch_energy_overEv[total_channels];
for (UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++)
ch_energy_overEv[ch_iter] = 0.;
Int_t Ev_counter = 0;
Int_t Fired_psd_channels_in_ev = 0;
TH1* th1_hist_ptr = nullptr;
TH2* th2_hist_ptr = nullptr;
th1_hist_ptr =
new TH1F("Energy_profile",
"PSD energy profile per event; cell # []; Energy [adc counts]",
total_channels,
0,
total_channels);
th1_hist_ptr = new TH1F("Energy_distrib",
"PSD energy distribution; Energy over event all "
"channels [adc counts]; Counts []",
1000,
0,
10000);
th1_hist_ptr = new TH1F("Spill_structure",
"PSD time spill structure; Time [s]; Counts []",
600,
20,
40);
th1_hist_ptr = new TH1F("Digi_time_diff",
"PSD digi time difference; Time [ns]; Counts []",
(Int_t) 500,
0,
500);
th1_hist_ptr =
new TH1F("Fired_psd_chs",
"Fired PSD channels in event; fired channels []; Counts []",
total_cells,
0,
total_cells);
th1_hist_ptr = new TH1F("Cosmic_check",
"Half a sum from two hodo channels; Energy []; Counts []",
100,
0,
100);
th2_hist_ptr = new TH2F("PsdVsHodo",
"Energy in Psd vs energy in Hodo; Psd energy [adc "
"counts]; Hodo energy [adc counts]",
1000,
0,
2000,
140,
0,
70);
for (UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++)
th1_hist_ptr = new TH1F(
Form("Energy_distrib_ch%i", ch_iter),
Form("Energy_distrib_ch%i; Energy [adc counts]; Counts []", ch_iter),
100,
0,
2000);
for (UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++)
th1_hist_ptr = new TH1F(
Form("Energy_distrib_good_ev_ch%i", ch_iter),
Form("Energy_distrib_good_ev_ch%i; Energy [adc counts]; Counts []",
ch_iter),
100,
0,
2000);
Double_t event_time = -1.;
Double_t prev_event_time = -1.;
Bool_t IsFirstEvent = true;
for (UInt_t ev = 0; ev < inTree->GetEntries(); ev++) {
inTree->GetEntry(ev);
Int_t nPsdDigis = fPsdDigis->size();
if (ev % 10 == 0)
std::cout << "Entry " << ev << " / " << inTree->GetEntries()
<< " Psd Digis: " << nPsdDigis << std::endl;
if (nPsdDigis == 0) continue;
CbmPsdDigi digi;
Bool_t IsNewEvent = false;
for (Int_t iDigi = 0; iDigi < nPsdDigis; iDigi++) {
digi = (*fPsdDigis)[iDigi];
//if(! (digi.GetSectionID() == 0 ) ) continue;
//if(digi.GetSectionID() > 8) continue;
if (digi.GetTime() - event_time > event_time_gate) {
prev_event_time = event_time;
event_time = digi.GetTime();
IsNewEvent = true;
} else
IsNewEvent = false;
if (IsNewEvent) {
if (IsFirstEvent) {
IsFirstEvent = false;
continue;
} //skip info before first event
Bool_t GoodEvent = true;
if (GoodEvent) {
th1_hist_ptr =
((TH1*) (gDirectory->FindObjectAny("Spill_structure")));
th1_hist_ptr->Fill(event_time / 1e9);
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Digi_time_diff")));
th1_hist_ptr->Fill(event_time - prev_event_time);
Double_t psd_energy = 0.;
Fired_psd_channels_in_ev = 0;
for (UInt_t ch_iter = 0; ch_iter < total_cells; ch_iter++) {
ch_energy_overEv[ch_iter] += ch_energy[ch_iter];
psd_energy += ch_energy[ch_iter];
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny(
Form("Energy_distrib_good_ev_ch%i", ch_iter))));
if (ch_energy[ch_iter] > 1.) th1_hist_ptr->Fill(ch_energy[ch_iter]);
if (ch_energy[ch_iter] > 1.) Fired_psd_channels_in_ev++;
ch_energy[ch_iter] = 0.;
}
//if(Fired_psd_channels_in_ev < 2) continue;
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Energy_distrib")));
th1_hist_ptr->Fill(psd_energy);
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Fired_psd_chs")));
th1_hist_ptr->Fill(Fired_psd_channels_in_ev);
/*
Double_t hodo_energy = 0.5*(ch_energy[17] + ch_energy[19]);
ch_energy_overEv[17] += ch_energy[17];
ch_energy_overEv[19] += ch_energy[19];
th1_hist_ptr = ((TH1*)(gDirectory->FindObjectAny( "Cosmic_check" )));
if(ch_energy[17] > 1. && ch_energy[19] > 1.) th1_hist_ptr->Fill(0.5*(ch_energy[17] + ch_energy[19]));
ch_energy[17] = 0.;
ch_energy[19] = 0.;
th2_hist_ptr = ((TH2*)(gDirectory->FindObjectAny( "PsdVsHodo" )));
th2_hist_ptr->Fill(psd_energy, hodo_energy);
*/
Ev_counter++;
} //is good event
else {
for (UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++)
ch_energy[ch_iter] = 0.;
}
} //is new event
Double_t Edep = digi.GetEdep();
UInt_t SectionID = digi.GetSectionID();
//if(ch_energy[SectionID] > 1.) cout<<"PileUp candidate in channel "<<SectionID<<endl;
ch_energy[SectionID] = Edep;
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny(
Form("Energy_distrib_ch%i", SectionID))));
th1_hist_ptr->Fill(Edep);
} //digi loop
} //event loop
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Energy_profile")));
for (UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++)
th1_hist_ptr->Fill(ch_iter, ch_energy_overEv[ch_iter] / Ev_counter);
TCanvas* canvas_time_spill = new TCanvas("c_time_spill", "c_time_spill");
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Spill_structure"))); th1_hist_ptr->Draw();
TCanvas* canvas_time_diff_ev =
new TCanvas("c_time_diff_ev", "c_time_diff_ev");
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Digi_time_diff")));
th1_hist_ptr->Draw();
TCanvas* canvas_energy_pr = new TCanvas("c_energy_pr", "c_energy_pr");
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Energy_profile")));
th1_hist_ptr->Draw("hist");
TCanvas* canvas_energy_dstr = new TCanvas("c_energy_dstr", "c_energy_dstr");
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Energy_distrib")));
th1_hist_ptr->Draw();
TCanvas* canvas_fired_chs = new TCanvas("c_fired_ch", "c_fired_ch");
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Fired_psd_chs")));
th1_hist_ptr->Draw();
TCanvas* canvas_psd_hodo = new TCanvas("c_psd_hodo", "c_psd_hodo");
th2_hist_ptr = ((TH2*) (gDirectory->FindObjectAny("PsdVsHodo")));
th2_hist_ptr->Draw("colz");
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny("Cosmic_check")));
th1_hist_ptr->Draw();
TCanvas* canvas_energy_dstr_ch =
new TCanvas("c_energy_dstr_ch", "c_energy_dstr_ch");
canvas_energy_dstr_ch->DivideSquare(total_channels);
for (UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) {
canvas_energy_dstr_ch->cd(ch_iter + 1);
th1_hist_ptr = ((
TH1*) (gDirectory->FindObjectAny(Form("Energy_distrib_ch%i", ch_iter))));
th1_hist_ptr->Draw();
}
TCanvas* canvas_energy_dstr_gEv_ch =
new TCanvas("c_energy_dstr_gEv_ch", "c_energy_dstr_gEv_ch");
canvas_energy_dstr_gEv_ch->DivideSquare(total_channels);
for (UInt_t ch_iter = 0; ch_iter < total_channels; ch_iter++) {
canvas_energy_dstr_gEv_ch->cd(ch_iter + 1);
th1_hist_ptr = ((TH1*) (gDirectory->FindObjectAny(
Form("Energy_distrib_good_ev_ch%i", ch_iter))));
th1_hist_ptr->Draw();
}
}