DataReader.cpp 13 KB
Newer Older
1
 /** @file DataReader.cxxs
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 *  @brief Implementation of DataReader.
 *
 *  Function definitions for DataReader are provided. 
 *  This class reads a rootfile with raw waveforms 
 *  that are processed by rcdaqAnalysis running on prdf files.
 *  Then, in the event loop, analysis classes can be called.
 *
 *  @author Yakov Kulinich, Riccardo Longo
 *  @bug No known bugs.
 */

#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TCanvas.h>
17
#include <TChain.h>
18
19
20
21

#include <iostream>

#include "DataReader.h"
Yakov Kulinich's avatar
Yakov Kulinich committed
22
#include "Analysis.h"
23
#include "Containers.h"
24
25
#include "RPD.h"
#include "ZDC.h"
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

/** @brief Default Constructor for DataReader.
 */
DataReader::DataReader() : DataReader( 0, 0, "", 0 ){

}

/** @brief Constructor for DataReader.
 *
 *  @param1 Number of channels being read
 *  @param2 Number of samples per channel
 */
DataReader::DataReader( const unsigned int nCh, const unsigned int nSamp )
  : DataReader( nCh, nSamp, "", 0 ){

  // here say we will read in list of files, maybe
  // there is better way to do it, just an idea for now.
  m_readListOfFiles = true;
}

/** @brief Constructor for DataReader.
 *
 *  @param1 Number of channels being read
 *  @param2 Number of samples per channel
 *  @param3 Input filename.
 */
DataReader::DataReader( const uint nCh, const uint nSamp,
			const std::string& fNameIn )
  : DataReader( nCh, nSamp, fNameIn, 0 ){
  
}


/** @brief Constructor for DataReader.
 *
 *  @param1 Number of channels being read
 *  @param2 Number of samples per channel
 *  @param4 Output file name (custom)
 *  @param3 Run number being used.
 
 */
DataReader::DataReader( const uint nCh, const uint nSamp,
			const std::string& fNameIn, const uint runNum )
Yakov Kulinich's avatar
Yakov Kulinich committed
69
70
71
  : m_nCh( nCh ), m_nSamp( nSamp ), m_fNameIn( fNameIn ),
    m_runNumber( runNum ), m_readListOfFiles( false ), m_fIn( NULL ){

72
73
74
75
76
77
}


/** @brief Destructor for DataReader.
 */
DataReader::~DataReader(){
Yakov Kulinich's avatar
Yakov Kulinich committed
78
79
80
81
82
83
84
85
86
87
88
89
90
91

  for( auto& ana : m_ana ){
    delete ana; ana = NULL;
  }
}

/** @brief Adds an analysis to vector of analysis
 *
 *  @param1 Pointer to an Analysis.
 *
 *  @return none
 */
void DataReader::AddAnalysis( Analysis* ana ){
  m_ana.push_back( ana );
92
93
}

94
95
96
97
98
99
100

/** @brief Enables the read from list of files option for DataReader
 *
 *  @param1 name of the list of files to be processed (with the full path if it's not in the execution folder)
 *
 *  @return none
 */
Yakov Kulinich's avatar
Yakov Kulinich committed
101
void DataReader::ReadListOfFiles( std::string listname ){
102
103
104
105
106

    m_readListOfFiles = true;
    m_fListOfFiles = listname;
}

Riccardo Longo's avatar
Riccardo Longo committed
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143

/**
 * @brief Reads the .xml configuration file and load characteristics for all the channels, immediately sorted into detectors objects
 * @param _inFile
 */
void DataReader::LoadAlignmentFile(std::string _inFile = "$JCaPA/Utils/Alignment_2018.xml"){

    m_XMLparser = new XMLSettingsReader();

    if (!m_XMLparser->parseFile(_inFile)) {
            std::cerr << " Data Reader could not parse file : " << _inFile << std::endl;
            return;
    }

    std::cout << "Loading .xml Alignment File..." << std::endl;
    std::cout << "Found " << m_XMLparser->getBaseNodeCount("Alignment") << " alignment entries " << std::endl;
    std::cout << "Retrieving the information for run " << m_runNumber << std::endl;

    int run;
    for (unsigned int i = 0; i < m_XMLparser->getBaseNodeCount("Alignment"); i++) {
        m_XMLparser->getChildValue("Alignment",i,"run",run);
        if(run != m_runNumber) continue;
        std::cout << "Found Run Entry in Alignment file for run " << m_runNumber << std::endl;
        m_XMLparser->getChildValue("Alignment",i,"x_table",m_alignment->x_table);
        m_XMLparser->getChildValue("Alignment",i,"y_table",m_alignment->y_table);
        m_XMLparser->getChildValue("Alignment",i,"upstream_Det",m_alignment->upstream_Det);
        m_XMLparser->getChildValue("Alignment",i,"mid_Det",m_alignment->mid_Det);
        m_XMLparser->getChildValue("Alignment",i,"downstream_Det",m_alignment->downstream_Det);
        m_XMLparser->getChildValue("Alignment",i,"target_In",m_alignment->target_In);
        m_XMLparser->getChildValue("Alignment",i,"lead_In",m_alignment->lead_In);
        m_XMLparser->getChildValue("Alignment",i,"magnet_On",m_alignment->magnet_On);
    }

    if(m_alignment == NULL) std::cout << "WARNING: ALIGNMENT NOT FOUND!!!" << std::endl;
    return;
}

144
145
146
147
/**
 * @brief Reads the .xml configuration file and load characteristics for all the channels, immediately sorted into detectors objects
 * @param _inFile
 */
148
void DataReader::LoadConfigurationFile(std::string _inFile = "$JCaPA/Utils/ConfigFile2018.xml"){
149
150
151
152
153

    m_XMLparser = new XMLSettingsReader();

    if (!m_XMLparser->parseFile(_inFile)) {
            std::cerr << " Data Reader could not parse file : " << _inFile << std::endl;
154
            return;
155
156
157
158
159
    }

    std::cout << "Loading .xml Configuration File..." << std::endl;
    std::cout << "Found " << m_XMLparser->getBaseNodeCount("channel") << " channel entries " << std::endl;

160
    std::vector < Channel* > channelEntries;
161
162
163
    int first_run, last_run;

    for (unsigned int i = 0; i < m_XMLparser->getBaseNodeCount("channel"); i++) {
164
        Channel *buffer_ch;
165
166
167
168
        m_XMLparser->getChildValue("channel",i,"start_run",first_run);
        m_XMLparser->getChildValue("channel",i,"end_run",last_run);

        //Discard entries for any channel that does not apply to our run
169
        if(m_runNumber < first_run || m_runNumber > last_run) continue;
170
171

        //If the entry applies, we store it in the vector
172
173
174
175
176
177
178
179
180
        m_XMLparser->getChildValue("channel",i,"detector",buffer_ch->detector);
        m_XMLparser->getChildValue("channel",i,"name",buffer_ch->name);
        m_XMLparser->getChildValue("channel",i,"mapping_row",buffer_ch->mapping_row);
        m_XMLparser->getChildValue("channel",i,"mapping_column",buffer_ch->mapping_column);
        m_XMLparser->getChildValue("channel",i,"delay",buffer_ch->delay);
        m_XMLparser->getChildValue("channel",i,"offset",buffer_ch->offset);
        m_XMLparser->getChildValue("channel",i,"HV",buffer_ch->HV);
        m_XMLparser->getChildValue("channel",i,"is_on",buffer_ch->is_on);
        m_XMLparser->getChildValue("channel",i,"Vop",buffer_ch->Vop);
181
182
183
184
185

        channelEntries.push_back(buffer_ch);
    }

    std::cout << "Loaded " << channelEntries.size() << " configuration entries " << std::endl;
186
    if( channelEntries.size() < 18 ) std::cout << "WARNING!!!! Number of Channels < 18. Seems that some entry is missed for this run in the config.xml. BE CAREFUL!" << std::endl;
187
188
189
190
    ZDC* zdc1 = new ZDC(channelEntries,1);
    ZDC* zdc2 = new ZDC(channelEntries,2);
    RPD* rpd = new RPD(channelEntries);

191
192
193
194
    m_detectors.push_back(zdc1); //Position 0 goes for ZDC1
    m_detectors.push_back(zdc2); //Position 1 goes for ZDC2
    m_detectors.push_back(rpd); //Position 2 goes for the RPD

195
    std::cout << "Detector configuration: loading complete! " << std::endl;
Riccardo Longo's avatar
Riccardo Longo committed
196

197
    return;
198
199
}

200
201
202
203
204
205
206
207
208
209
210
211
/**
 * @brief DataReader::GetDetector allows the user to access the detectors objects after loading them at the beginning of the execution
 * @param _detName can be ZDC1 (upstream module), ZDC2 (downstream module) or RPD.
 * @return
 */
Detector* DataReader::GetDetector( std::string _detName ){

    if(_detName != "ZDC1" && _detName != "ZDC2" && _detName != "zdc1" && _detName != "zdc2" && _detName != "RPD" && _detName != "rpd")
    {
        std::cout << "The detector you're looking for is not ZDC1, ZDC2 or RPD. Please check and correct your request " << std::endl;
        return NULL;
    }
Riccardo Longo's avatar
Riccardo Longo committed
212
213
214
    if( _detName == "ZDC1" || _detName == "zdc1" ) return m_detectors.at(0);
    if( _detName == "ZDC2" || _detName == "zdc2" ) return m_detectors.at(1);
    if( _detName == "RPD" || _detName == "rpd" ) return m_detectors.at(2);
215

Riccardo Longo's avatar
Riccardo Longo committed
216
217
    std::cout << "WARNING: detector recognition glitch. NULL pointer being returned..." << std::endl;
    return NULL;
218
219
}

220

221
222
223
224
225
226
227
228
229
230
231
232
233
234
/** @brief Run method for DataReader
 *
 *  Call Initialize, ProcessEvents, and Finalize
 *  Made so the driver class only has to call one method.
 *  
 *  @return none
 */
void DataReader::Run(){

  Initialize();
  ProcessEvents();
  Finalize();
}

235
236
237
238
239
240
241
242
243
244
245
246
/** @brief Initialization method for DataReader
 *
 *  Select which file(s) to read. For now just a single
 *  file, later this can be extended to read many and make
 *  chain of files for example. Also create and initialize
 *  the analysis that will be running.
 *  
 *  @return none
 */
void DataReader::Initialize(){

  if( m_readListOfFiles ){
247
248
249
250
251
252
253
254
255
256
257
258
259

    // Riccardo - 21/01/2019 - TChain implementation
    // The file list name is supposed to be already set
    m_fileChain = new TChain("tree");
    std::ifstream inFile;
    inFile.open(m_fListOfFiles.c_str());
    std::string reader_buff;
    while(inFile >> reader_buff){
        //Let's push all the files in the list into the TChain
        m_fileChain->Add(reader_buff.c_str());
    }
    /** TODO - Add fileChain reading below
     */
260
  } else {
Yakov Kulinich's avatar
Yakov Kulinich committed
261
    m_fIn = TFile::Open( m_fNameIn.c_str() );
262
  }
Yakov Kulinich's avatar
Yakov Kulinich committed
263
264
265
266
267
268
269
270

  // If we are reading a list of files, or have no run number
  // make default name output.root, otherwise make it
  // outputN.root, where N is a run number of a file.
  std::string fNameOut = m_readListOfFiles ?
    "output.root" : Form( "output%d.root", m_runNumber );

  m_fOut = new TFile( fNameOut.c_str(), "RECREATE" );
271
  
Yakov Kulinich's avatar
Yakov Kulinich committed
272
273
274
275
  for( auto& ana : m_ana ){
    ana->Initialize();
    ana->SetupHistograms();
  }
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
}



/** @brief Process Events method for DataReader
 *
 *  Connect to files / chain of files.
 *  Read tree event by event, and fill basic waveforms.
 *  So far (12.12.18) For each event - Read Raw data for all channels,
 *  put it into 2D vector size NxM where N = nCh, M = nSamp
 *  Also, For each event - Fill N 1D histos that have size M.
 *  Then one can send these to any Analysis function they make.
 *  Here for example (12.12.18), we send to WFAnalysis::AnalyzeEvent( .. )
 *  See below - Have fun!
 *  
 *  @return none
 */
void DataReader::ProcessEvents(){

295
  // Processed Raw data to read in as vector of vectors size NxM
Yakov Kulinich's avatar
Yakov Kulinich committed
296
  // Where N = nCh and M = nSamples per channel.
297
298
  std::vector< std::vector< float >  >  vWF;
  std::vector< std::vector< float >* > pvWF;
Yakov Kulinich's avatar
Yakov Kulinich committed
299
300
301
302
303
304

  // Histograms (N of them) for the raw waveforms from each event.
  // They will go to AnalyzeEvent for processing
  std::vector< TH1* > vWFH;

  // Resize these to be of size nCh.
305
306
307
  vWF .resize( m_nCh );
  pvWF.resize( m_nCh );
  vWFH.resize( m_nCh );
Yakov Kulinich's avatar
Yakov Kulinich committed
308

309
310
311
  /** TODO : add reading for list of files
    * Please note that many of the implementations are now for a single-file treatment
    */
Yakov Kulinich's avatar
Yakov Kulinich committed
312
  TTree* tree = static_cast< TTree* >( m_fIn->Get( "tree" ) );
313

314
315
316
317
318
319
320
321
322
323
324
325
326
  //Specific pointers to each detector, if needed afterwards
  ZDC* zdc1 = static_cast< ZDC* >( GetDetector("ZDC1") );
  ZDC* zdc2 = static_cast< ZDC* >( GetDetector("ZDC2") );
  RPD* rpd =  static_cast< RPD* >( GetDetector("RPD") );

  //All the raw channels addresses set for read-out
  for( uint detID = 0; detID < (int) m_detectors.size(); detID++ ){
       m_detectors.at(detID)->SetBranches(tree);
       m_detectors.at(detID)->SetNSamples(m_nSamp);
       m_detectors.at(detID)->DeclareHistograms();
  }

  // Connect raw data to tree. For the moment, the only reading implemented is the raw data from each channel.
327
  // Other items should be implemented in the same way if needed.
328
329
  // Also create the histograms to fill
  for( uint ch = 0; ch < m_nCh; ch++ ){
330
    //Example - here we retrieve the already processed waveform (M.Phipps approach during the data production)
331
    pvWF[ ch ] = &vWF[ ch ];
332
    tree->SetBranchAddress( Form( "C%d", ch ), &pvWF[ ch ] );
Yakov Kulinich's avatar
Yakov Kulinich committed
333
    vWFH[ ch ] = new TH1D( Form( "hWF%d", ch ), Form( "hWF%d;samp;amp", ch ), m_nSamp, 0, m_nSamp );
334
  }
335
336

  std::cout << "File: " << m_fIn->GetName() << " has " << tree->GetEntries() << " events." << std::endl;
337
338
339
340
341
  
  // !! EVENT LOOP
  for( int ev = 0; ev < tree->GetEntries(); ev++ ){
    tree->GetEntry( ev );

342
343
344
    // Fill the raw waveforms
    for( uint detID = 0; detID < (int) m_detectors.size(); detID++ )
        m_detectors.at(detID)->FillHistograms();
345

346
   //Here if you're interested in already processed waveform
347
   for( uint ch = 0; ch < m_nCh; ch++ ) {
348
   // Loop over samples in each channel
349
350
351
352
353
    for( uint samp = 0; samp < m_nSamp; samp++ ){
      vWFH[ ch ]->SetBinContent( samp + 1, vWF[ ch ][ samp ] );
      } // End loop over samples in each channel
    } // End loop over channels

Yakov Kulinich's avatar
Yakov Kulinich committed
354
    // Now call all analysis and run their AnalyzeEvent.
355
356
    // Can either send a vector of histos, or a 2D vector
    // Of all the data, depending on what you want to do.
357
    // Note that at the moment none of this methods is doing anything
Yakov Kulinich's avatar
Yakov Kulinich committed
358
    for( auto& ana : m_ana ){
359
360
361
362
363
      //raw data analysis
      ana->AnalyzeEvent( zdc1->GetChannelsVector()  );
      ana->AnalyzeEvent( zdc2->GetChannelsVector() );
      ana->AnalyzeEvent( rpd->GetChannelsVector() );
      //already processed wf analysis
Yakov Kulinich's avatar
Yakov Kulinich committed
364
365
366
      ana->AnalyzeEvent( vWFH );
      ana->AnalyzeEvent( vWF  );
    }
367
368
369
370
371
372
373
374
375
376
377
378
379
  } // End event loop

  for( auto& h : vWFH ){ delete h; }
}

/** @brief Finalize method for DataReader
 *
 *  Close any input files 
 *  Call analysis finalize methods
 *  
 *  @return none
 */
void DataReader::Finalize(){
Yakov Kulinich's avatar
Yakov Kulinich committed
380
  
381
382
383
384
  if( m_fIn ){
    m_fIn->Close();
  }

Yakov Kulinich's avatar
Yakov Kulinich committed
385
386
387
388
389
390
391
392
393
  // enter the output file since
  // we will be writing to it now.
  m_fOut->cd();

  for( auto& ana : m_ana ){
    ana->Finalize();
  }

  m_fOut->Close();
394
395
}