DataReader.cpp 13.9 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

/**
 * @brief Reads the .xml configuration file and load characteristics for all the channels, immediately sorted into detectors objects
 * @param _inFile
 */
112
void DataReader::LoadAlignmentFile(std::string _inFile ){
Riccardo Longo's avatar
Riccardo Longo committed
113
114
115
116
117
118
119
120

    m_XMLparser = new XMLSettingsReader();

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

121
122
    m_alignment = new Alignment();

Riccardo Longo's avatar
Riccardo Longo committed
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
    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;
}

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

    m_XMLparser = new XMLSettingsReader();

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

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

162
    std::vector < Channel* > channelEntries;
163
164
165
    int first_run, last_run;

    for (unsigned int i = 0; i < m_XMLparser->getBaseNodeCount("channel"); i++) {
166
        Channel *buffer_ch = new Channel();
167
168
169
170
        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
171
        if(m_runNumber < first_run || m_runNumber > last_run) continue;
172

173

174
        //If the entry applies, we store it in the vector
175
176
177
178
179
180
181
182
183
        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);
184
185
186
187
188

        channelEntries.push_back(buffer_ch);
    }

    std::cout << "Loaded " << channelEntries.size() << " configuration entries " << std::endl;
189
    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;
190
191
192
193
    ZDC* zdc1 = new ZDC(channelEntries,1);
    ZDC* zdc2 = new ZDC(channelEntries,2);
    RPD* rpd = new RPD(channelEntries);

194
195
196
197
    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

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

200
    return;
201
202
}

203
204
205
206
207
208
209
210
211
212
213
214
/**
 * @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
215
216
217
    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);
218

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

223

224
225
226
227
228
229
230
231
232
233
234
235
236
237
/** @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();
}

238
239
240
241
242
243
244
245
246
247
248
249
/** @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 ){
250
251
252
253
254
255
256
257
258
259
260
261
262

    // 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
     */
263
  } else {
Yakov Kulinich's avatar
Yakov Kulinich committed
264
    m_fIn = TFile::Open( m_fNameIn.c_str() );
265
  }
Yakov Kulinich's avatar
Yakov Kulinich committed
266
267
268
269
270
271
272
273

  // 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" );
274
  
Yakov Kulinich's avatar
Yakov Kulinich committed
275
276
277
278
  for( auto& ana : m_ana ){
    ana->Initialize();
    ana->SetupHistograms();
  }
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
}



/** @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(){
clantz's avatar
clantz committed
297
298
  
  //Added for debug/demo purposes. Remove after implementation of Visualizer
299
  TCanvas *canvas = new TCanvas( "Diff Demo", "Diff Demo", 200, 10, 1000, 600);
300
301
  TPad *pad = new TPad("pad", "pad",0.15,0.11,0.85,0.79);
  canvas->Divide(4,2);
clantz's avatar
clantz committed
302
  
303
  // Processed Raw data to read in as vector of vectors size NxM
Yakov Kulinich's avatar
Yakov Kulinich committed
304
  // Where N = nCh and M = nSamples per channel.
305
306
  std::vector< std::vector< float >  >  vWF;
  std::vector< std::vector< float >* > pvWF;
Yakov Kulinich's avatar
Yakov Kulinich committed
307
308
309
310
311
312

  // 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.
313
314
315
  vWF .resize( m_nCh );
  pvWF.resize( m_nCh );
  vWFH.resize( m_nCh );
Yakov Kulinich's avatar
Yakov Kulinich committed
316

317
318
319
  /** 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
320
  TTree* tree = static_cast< TTree* >( m_fIn->Get( "tree" ) );
321

322
323
324
325
326
327
328
329
330
331
332
333
334
  //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.
335
  // Other items should be implemented in the same way if needed.
336
337
  // Also create the histograms to fill
  for( uint ch = 0; ch < m_nCh; ch++ ){
338
    //Example - here we retrieve the already processed waveform (M.Phipps approach during the data production)
339
    pvWF[ ch ] = &vWF[ ch ];
340
    tree->SetBranchAddress( Form( "C%d", ch ), &pvWF[ ch ] );
Yakov Kulinich's avatar
Yakov Kulinich committed
341
    vWFH[ ch ] = new TH1D( Form( "hWF%d", ch ), Form( "hWF%d;samp;amp", ch ), m_nSamp, 0, m_nSamp );
342
  }
343
344

  std::cout << "File: " << m_fIn->GetName() << " has " << tree->GetEntries() << " events." << std::endl;
345
346
  
  // !! EVENT LOOP
347
348
  for( int ev = 0; ev < tree->GetEntries(); ev++ ){
    
clantz's avatar
clantz committed
349

350
  
351
    tree->GetEntry( ev );
352
    
353

354
355
356
    // Fill the raw waveforms
    for( uint detID = 0; detID < (int) m_detectors.size(); detID++ )
        m_detectors.at(detID)->FillHistograms();
357

358
   //Here if you're interested in already processed waveform
359
   for( uint ch = 0; ch < m_nCh; ch++ ) {
360
   // Loop over samples in each channel
361
362
363
364
365
    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
366
    // Now call all analysis and run their AnalyzeEvent.
367
368
    // Can either send a vector of histos, or a 2D vector
    // Of all the data, depending on what you want to do.
369
    // Note that at the moment none of this methods is doing anything
Yakov Kulinich's avatar
Yakov Kulinich committed
370
    for( auto& ana : m_ana ){
371
      //raw data analysis
clantz's avatar
clantz committed
372
373
374
375
376
377
378
379
380
381
382
383
384
385
      if(m_debug){    
        // Uncomment to run a few events at a time
        //if(ev==8) break;
    
        // Uncomment to run a single event
        if(ev!=8) continue;
        
        //ana->AnalyzeEvent( zdc1->GetChannelsVector(), canvas->cd() );
        //ana->AnalyzeEvent( zdc2->GetChannelsVector(), canvas->cd() );
        ana->AnalyzeEvent( rpd->GetChannelsVector() , canvas->cd() );
      }else{
      ana->AnalyzeEvent( zdc1->GetChannelsVector() );
      ana->AnalyzeEvent( zdc2->GetChannelsVector() );
      ana->AnalyzeEvent( rpd->GetChannelsVector()  );
386
      //already processed wf analysis
387
388
      //ana->AnalyzeEvent( vWFH );
      //ana->AnalyzeEvent( vWF  );
clantz's avatar
clantz committed
389
      }
Yakov Kulinich's avatar
Yakov Kulinich committed
390
    }
391
  } // End event loop
392
  
clantz's avatar
clantz committed
393
394
395
396
397
398
399
  
  //Added for debug/demo purposes. Remove after implementation of Visualizer
  if(m_debug){
    pad->Update();
    canvas->Draw();
    canvas->Print( "Output.pdf" );
  }
400
401
402
403
404
405
406
407
408
409
410
411

  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
412
  
413
414
415
416
  if( m_fIn ){
    m_fIn->Close();
  }

Yakov Kulinich's avatar
Yakov Kulinich committed
417
418
419
420
421
422
423
424
425
  // 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();
426
427
}