-
Yftach Moyal authored
Added functionality of parsing LHCf driven files and analyzing them.\n Added a check on is_on status for Detector::SetBranches and Detector::FillHistogram.\n Added Dummy channels for UEM in H4 LHCf section.
Yftach Moyal authoredAdded functionality of parsing LHCf driven files and analyzing them.\n Added a check on is_on status for Detector::SetBranches and Detector::FillHistogram.\n Added Dummy channels for UEM in H4 LHCf section.
LHCf_reader.cpp 15.85 KiB
/** @defgroup ana Analysis
* @ingroup ana
* @file LHCf_reader.cpp
* @brief Implementation of LHCf_reader.
*
* Function definitions for LHCf_reader are provided.
* This class reads a rootfile with raw waveforms stored in LHCf format
* that are processed by LHCf DAQ.
* The first scope is just to provide a 1-1 conversion to our pdrf format
*
* @author Riccardo Longo, Yftach Moyal
* @bug No known bugs.
*/
#include "LHCf_reader.h"
#include <string>
#include <stdio.h>
/** @brief Constructor for LHCf_reader.
*/
LHCf_reader::LHCf_reader( const std::string& _fileName,
const unsigned int _nCh,
const unsigned int _nSamples )
: DataReader2021( _nCh, _nSamples, _fileName, 0 )
{
//m_fileName = _fileName;
//m_nCh = 32;
}
/** @brief Destructor for LHCf_reader.
*/
LHCf_reader::~LHCf_reader( ){
}
void LHCf_reader::InitializeHistograms(){
std::string chName;
for( int i = 0; i < m_nCh; i++){
chName = Form( "Channel%d", i );
m_v_histoChannels.push_back( new TH1D(chName.c_str(), chName.c_str(), m_nSamp, 0, m_nSamp));
chName = Form( "ChannelCorr%d", i );
m_v_histoChannelsCorr.push_back( new TH1D(chName.c_str(), chName.c_str(), m_nSamp, 0, m_nSamp));
}
}
void LHCf_reader::OpenFile( void ){
std::cout << "Opening file " << m_fNameIn << std::endl;
if ( m_fNameIn == "ND"){
std::cerr << "File name not defined. Going to crash!" << std::endl;
return;
}
m_filePointer = TFile::Open( m_fNameIn.c_str() );
m_readTrigger1 = new LHCf_Trigger1((TTree*)m_filePointer->Get( "Trigger1" ));
InitializeHistograms();
std::cout << "File " << m_fNameIn << " opened successfully! " << std::endl;
}
void LHCf_reader::OpenFileName( std::string _fileName ){
std::cout << "Opening file " << m_fNameIn << std::endl;
m_fNameIn = _fileName;
m_filePointer = TFile::Open( m_fNameIn.c_str() );
//m_trigger1 =
m_readTrigger1 = new LHCf_Trigger1((TTree*)m_filePointer->Get( "Trigger1" ));
InitializeHistograms();
std::cout << "File " << m_fNameIn << " opened successfully! " << std::endl;
}
// Get 3 32 bit dwords and return 8 Ch sample values of the group Ch 0,1,2…7
std::vector<float> LHCf_reader::GetGroupSample(unsigned int dwrd_0, unsigned int dwrd_1, unsigned int dwrd_2){
/*
Get 3 32 bit dwords and return 8 Ch sample values
of the group Ch 0,1,2…7.
Details:
To shift x by n bits to the left (make it bigger):
x*2^n
===============================================
To read the m bits n positions from the beginning of the word x
x/(2^n)%2^m
===============================================
*/
std::vector<float> sample;
sample.push_back(dwrd_0%4096); // CH0
sample.push_back((dwrd_0/4096)%4096); // CH1
sample.push_back((dwrd_0/16777216)%4096 ^ (dwrd_1%16)*256); // CH2
sample.push_back((dwrd_1/16)%4096); // CH3
sample.push_back((dwrd_1/65536)%4096); // CH4
sample.push_back((dwrd_1/268435456)%4096 ^ (dwrd_2%256)*16); // CH5
sample.push_back((dwrd_2/256)%4096); // CH6
sample.push_back((dwrd_2/1048576)%4096); // CH7
return sample;
}
void LHCf_reader::ShowEvent( int iEvent ){
m_readTrigger1->Show( iEvent );
}
void LHCf_reader::ConvertFile(){
std::string outputDir = Form("%s/ProcessedLHCf",std::getenv("JZCaPA"));
gSystem->Exec( Form("mkdir -p %s", outputDir.c_str() ) );
std::string base_filename = m_fNameIn.substr(m_fNameIn.find_last_of("/\\") + 1);
std::string of = Form("%s/%s",outputDir.c_str() ,base_filename.c_str() );
if (m_debug) printf(" output file:%s\n",of.c_str());
TFile* out_file = TFile::Open(of.c_str(),"RECREATE");
TObject *obj;
TKey *key;
TTree* tree[5];
TIter next( m_filePointer->GetListOfKeys());
int id = 0;
bool readCurrent = 0;
int RunNumber;
int evt;
// Get runNumber
std::string fname_str(base_filename.c_str());
std::string result = fname_str.substr(0, fname_str.find_last_of("."));
RunNumber = atoi(result.substr(result.find_last_not_of("0123456789") + 1).c_str());
for ( int iCh = 0; iCh < m_nCh; iCh++){
std::vector < float > v_buffer, v_buffer2;
v_buffer.resize( m_nSamp );
v_buffer2.resize( m_nSamp );
m_v_Channels.push_back( v_buffer );
m_v_Channels_Corrected.push_back( v_buffer2 );
}
while ((key = (TKey *) next())) {
obj = m_filePointer->Get(key->GetName()); // copy object to memory
if (strcmp(key->GetClassName(),"TTree")) continue;
if (m_debug) printf(" found object:%s\n",key->GetName()); if (!strcmp(((TTree*)obj)->GetName(),"Trigger1")){
if (readCurrent) continue;
readCurrent = 1;
// Do Things with Trigger1
TTree* tree_to_read = (TTree*)obj;
if (m_debug) printf(" Nubmer of entries:%lld\n",tree_to_read->GetEntries());
// Deactivate the branch which needs to be processed
tree_to_read->SetBranchStatus("ZDCH",0);
// Clone skimmed
if (m_debug) printf(" About to clone LHCf branches\n");
tree[id] = tree_to_read->CloneTree(0);
tree[id]->CopyEntries(tree_to_read);
tree[id]->SetBranchStatus("*",0);
if (m_debug) printf(" Successfully cloned LHCf branches\n");
// Reactivate ZDCH
tree_to_read->SetBranchStatus("*",0);
tree_to_read->SetBranchStatus("ZDCH",1);
//Create run number and event number branches
tree[id]->Branch("RunNumber", &RunNumber, "RunNumber/I" );
tree[id]->Branch("EventNumber", &evt, "EventNumber/I" );
// Generate new Branches for unpacking
if (m_debug) printf(" About to generate ZDCH WF branches\n");
for (int ch=0;ch<m_nChPerGroup*m_nGroups;ch++){
tree[id]->Branch( Form( "RawSignal%d", ch), "std::vector<float>", &m_v_Channels_Corrected.at(ch) );
}
if (m_debug) printf(" Finished generating ZDCH WF branches\n");
// Run over all events and unpack, as per the 2 event shift
int MaxEvents = tree_to_read->GetEntries();
for(evt=0;evt<MaxEvents;evt++){
// Get Samples from evt+2
if ((MaxEvents-evt) < 3) continue;// Fill Empty
GetEventWaveforms(evt+2);
tree[id]->Fill();
}
tree[id]->SetEntries(MaxEvents);
tree[id]->SetName("tree"); // Temporary to see if solves JZCaPA issue
}else{
if (m_debug) printf(" Writing to TTree with ID: %d\n",id);
tree[id] = ((TTree*)obj)->CloneTree();
}
id++;
}
if (m_filePointer) m_filePointer->Close();
if (m_debug) printf(" Moving into new file.\n");
out_file->cd();
if (m_debug) printf(" Writing into new file.\n");
out_file->Write();
if (m_debug) printf(" Finished writing.\n");
out_file->Close();
}
void LHCf_reader::GetEventWaveforms(int iEvent){
m_v_Channels.clear();
m_v_Channels_Corrected.clear();
m_readTrigger1->GetEntry( iEvent );
bool trg_beam = false;
bool trg_pede = false;
trg_beam = (m_readTrigger1->GPI0_GPI0[0] >> 8) & 0x1;
trg_pede = (m_readTrigger1->GPI0_GPI0[0] >> 5) & 0x1;
//Checking trigger ADCs
double adc[2];
adc[0] = m_readTrigger1->ADC2_ADC2[0] - m_readTrigger1->ADC2_ADC2[32];// trg.scin 1 adc[1] = m_readTrigger1->ADC2_ADC2[2] - m_readTrigger1->ADC2_ADC2[34];// trg.scin 2
if(trg_pede) {
std::cout << "PEDESTAL TRIGGER - SKIP " << std::endl;
// return;
}
if(trg_beam){
std::cout << "BEAM TRIGGER " << std::endl;
}
for ( int iCh = 0; iCh < m_nCh; iCh++){
std::vector < float > v_buffer, v_buffer2;
v_buffer.resize( m_nSamp );
v_buffer2.resize( m_nSamp );
m_v_Channels.push_back( v_buffer );
m_v_Channels_Corrected.push_back( v_buffer2 );
}
//Retrieving calibration file from Util folder
std::string calibFile = std::getenv("JZCaPA") + std::string("/Utils/calib_0234_2G_0.dat");
caen_correction *cc;
cc = new caen_correction ( calibFile.c_str() );
unsigned int curser = 4;
unsigned int n_words_for_group_data_extraction;
unsigned int c_stop_flag;
//NSI - Noise subtraction implementation
unsigned int index_cell[m_nGroups];
// For group 1
for (int gid = 0; gid < m_nGroups; gid++){
n_words_for_group_data_extraction = m_readTrigger1->ZDCH_ZDCH[curser]%4096;
index_cell[gid] = (m_readTrigger1->ZDCH_ZDCH[curser]/1048576)%4096;
if(m_debug)
std::cout << "GID: " << gid << " n_words_for_group_data_extraction: " << n_words_for_group_data_extraction << std::endl;
//std::cout << "TR: " << (m_readTrigger1->ZDCH_ZDCH[curser]/4096)%2 << std::endl;
curser += 1;
c_stop_flag = curser + n_words_for_group_data_extraction;
// run over the dwords and extract each sample
int iBin = 0;
if(m_debug)
std::cout << "PRE --> CURSER: " << curser << " - CSTOPFLAG: " << c_stop_flag << std::endl;
for (int dummy = 0; curser < c_stop_flag; curser+=3){
std::vector < float > sample_vec = GetGroupSample(m_readTrigger1->ZDCH_ZDCH[curser],m_readTrigger1->ZDCH_ZDCH[curser+1],m_readTrigger1->ZDCH_ZDCH[curser+2]);
// store samples
for ( int iSamp = 0; iSamp < m_nChPerGroup; iSamp++){
m_v_Channels.at( (gid * m_nChPerGroup) + iSamp ).at(iBin) = sample_vec.at(iSamp);
}
iBin++;
}
if(m_debug)
std::cout << "POST --> CURSER: " << curser << " - CSTOPFLAG: " << c_stop_flag << std::endl;
// Move to next group data words
curser += 1;
}
// Perform signal correction using firmware values
cc->init(index_cell, m_nGroups, m_v_Channels);
for(int iCh = 0; iCh < m_nChPerGroup*m_nGroups; iCh++){
for(int iSamp = 0; iSamp < m_nSamp; iSamp++){
m_v_Channels_Corrected.at( iCh ).at( iSamp ) = cc->caen_corrected( iSamp, iCh );
}
}
}
void LHCf_reader::ReadEvents( int iStartEvent, int iEndEvent ){
if( iStartEvent == -1 ) iStartEvent = 0;
if( iEndEvent == -1 ) iEndEvent = m_readTrigger1->GetEntries();
if( iEndEvent-iStartEvent > GetMaxEvents() && iEndEvent != -1
&& GetMaxEvents() != -1 && iStartEvent != -1)
iEndEvent = iStartEvent+GetMaxEvents();
TCanvas* ct = new TCanvas("ct", "ct", 3200, 4800);
ct->Divide(4,m_nGroups*2);
for( int iEvent = iStartEvent ; iEvent < iEndEvent; iEvent++){
m_readTrigger1->GetEntry( iEvent );
bool trg_beam = false;
bool trg_pede = false;
trg_beam = (m_readTrigger1->GPI0_GPI0[0] >> 8) & 0x1;
trg_pede = (m_readTrigger1->GPI0_GPI0[0] >> 5) & 0x1;
//Checking trigger ADCs
double adc[2];
adc[0] = m_readTrigger1->ADC2_ADC2[0] - m_readTrigger1->ADC2_ADC2[32];// trg.scin 1
adc[1] = m_readTrigger1->ADC2_ADC2[2] - m_readTrigger1->ADC2_ADC2[34];// trg.scin 2
if(trg_pede) {
std::cout << "PEDESTAL TRIGGER - SKIP " << std::endl;
continue;
}
if(trg_beam){
std::cout << "BEAM TRIGGER " << std::endl;
}
for ( int iCh = 0; iCh < m_nCh; iCh++){
std::vector < float > v_buffer, v_buffer2;
v_buffer.resize( m_nSamp );
v_buffer2.resize( m_nSamp );
m_v_Channels.push_back( v_buffer );
m_v_Channels_Corrected.push_back( v_buffer2 );
}
//Retrieving calibration file from Util folder
std::string calibFile = std::getenv("JZCaPA") + std::string("/Utils/calib_0234_2G_0.dat");
caen_correction *cc;
cc = new caen_correction ( calibFile.c_str() );
int nPerGroup = 8;
unsigned int curser = 4;
unsigned int n_words_for_group_data_extraction;
unsigned int c_stop_flag;
//NSI - Noise subtraction implementation
unsigned int index_cell[m_nGroups];
// For group 1
for (int gid = 0; gid < m_nGroups; gid++){
n_words_for_group_data_extraction = m_readTrigger1->ZDCH_ZDCH[curser]%4096;
index_cell[gid] = (m_readTrigger1->ZDCH_ZDCH[curser]/1048576)%4096;
if(m_debug)
std::cout << "GID: " << gid << " n_words_for_group_data_extraction: " << n_words_for_group_data_extraction << std::endl;
//std::cout << "TR: " << (m_readTrigger1->ZDCH_ZDCH[curser]/4096)%2 << std::endl;
curser += 1;
c_stop_flag = curser + n_words_for_group_data_extraction;
// run over the dwords and extract each sample
int iBin = 0;
if(m_debug)
std::cout << "PRE --> CURSER: " << curser << " - CSTOPFLAG: " << c_stop_flag << std::endl;
for (int dummy = 0; curser < c_stop_flag; curser+=3){
std::vector < float > sample_vec = GetGroupSample(m_readTrigger1->ZDCH_ZDCH[curser],m_readTrigger1->ZDCH_ZDCH[curser+1],m_readTrigger1->ZDCH_ZDCH[curser+2]); // store samples
for ( int iSamp = 0; iSamp < nPerGroup; iSamp++){
m_v_Channels.at( (gid * nPerGroup) + iSamp ).at(iBin) = sample_vec.at(iSamp);
m_v_histoChannels.at( (gid * nPerGroup) + iSamp )->Fill( iBin, sample_vec.at(iSamp));
}
iBin++;
}
if(m_debug)
std::cout << "POST --> CURSER: " << curser << " - CSTOPFLAG: " << c_stop_flag << std::endl;
// Move to next group data words
curser += 1;
}
cc->init(index_cell, m_nGroups, m_v_Channels);
for(int iCh = 0; iCh < nPerGroup*m_nGroups; iCh++){
for(int iSamp = 0; iSamp < m_nSamp; iSamp++){
m_v_Channels_Corrected.at( iCh ).at( iSamp ) = cc->caen_corrected( iSamp, iCh );
m_v_histoChannelsCorr.at( iCh )->Fill( iSamp, m_v_Channels_Corrected.at( iCh ).at( iSamp )+400);
}
}
for( int iCh = 0; iCh < nPerGroup*m_nGroups; iCh++){
ct->cd( iCh+1 );
m_v_histoChannels.at(iCh)->Draw( "HISTO" );
m_v_histoChannelsCorr.at(iCh)->SetLineColor(kRed);
m_v_histoChannelsCorr.at(iCh)->Draw( "HISTOSAME" );
if(iCh < 8) m_v_histoChannels.at(iCh)->GetYaxis()->SetRangeUser( 3100, 4100 );
if(iCh >= 8) m_v_histoChannels.at(iCh)->GetYaxis()->SetRangeUser( 1800, 2800 );
}
/* The following block works only if you don't set range User above.
// Also - it's so basic to be otrageous
sci1Corr->Fill( adc[0], m_v_histoChannels.at( 0 )->GetMinimum( ) );
sci2Corr->Fill( adc[1], m_v_histoChannels.at( 1 )->GetMinimum( ) );
if(m_debug) std::cout << "In Time WF: " << m_v_histoChannels.at( 0 )->GetMinimum( ) << " - " << m_v_histoChannels.at( 1 )->GetMinimum( ) << std::endl;
if(m_debug) std::cout << "In Time ADC: " << adc[0] << " - " << adc[1] << std::endl;
if(iEvent > 0){
m_readTrigger1->GetEntry( iEvent-1 );
adc[0] = m_readTrigger1->ADC2_ADC2[0] - m_readTrigger1->ADC2_ADC2[32];// trg.scin 1
adc[1] = m_readTrigger1->ADC2_ADC2[2] - m_readTrigger1->ADC2_ADC2[34];// trg.scin 2
p_sci1Corr->Fill( adc[0], m_v_histoChannels.at( 0 )->GetMinimum( ) );
p_sci2Corr->Fill( adc[1], m_v_histoChannels.at( 1 )->GetMinimum( ) );
if(m_debug)
std::cout << "Previous Trigger ADC: " << adc[0] << " - " << adc[1] << std::endl;
}
if(iEvent > 1){
m_readTrigger1->GetEntry( iEvent-2 );
adc[0] = m_readTrigger1->ADC2_ADC2[0] - m_readTrigger1->ADC2_ADC2[32];// trg.scin 1
adc[1] = m_readTrigger1->ADC2_ADC2[2] - m_readTrigger1->ADC2_ADC2[34];// trg.scin 2
pp_sci1Corr->Fill( adc[0], m_v_histoChannels.at( 0 )->GetMinimum( ) );
pp_sci2Corr->Fill( adc[1], m_v_histoChannels.at( 1 )->GetMinimum( ) );
if(m_debug)
std::cout << "Previous x2 Trigger ADC: " << adc[0] << " - " << adc[1] << std::endl;
}
*/
std::string svname;
svname = Form( "Event%d.png", iEvent );
ct->Print(svname.c_str());
for( int iCh = 0; iCh < m_v_histoChannels.size(); iCh++){
m_v_histoChannels.at(iCh)->Reset();
m_v_histoChannelsCorr.at(iCh)->Reset();
}
}
}
void LHCf_reader::VisualizeChannels( ){
TCanvas* cc = new TCanvas("cc", "cc", 3200, 800);
cc->Divide(8,m_nGroups); for( int iCh = 0; iCh < m_v_histoChannels.size(); iCh++){
cc->cd( iCh+1 );
m_v_histoChannels.at(iCh)->Draw( "HISTO" );
m_v_histoChannels.at(iCh)->GetYaxis()->SetRangeUser( 0, 4096 );
}
}