DataReader.cpp 14.6 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
#include "Visualizer.h"
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
69

/** @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
70
71
72
  : m_nCh( nCh ), m_nSamp( nSamp ), m_fNameIn( fNameIn ),
    m_runNumber( runNum ), m_readListOfFiles( false ), m_fIn( NULL ){

73
74
75
76
77
78
}


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

  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 );
93
94
}

95
96
97
98
99
100
101

/** @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
102
void DataReader::ReadListOfFiles( std::string listname ){
103
104
105
106
107

    m_readListOfFiles = true;
    m_fListOfFiles = listname;
}

Riccardo Longo's avatar
Riccardo Longo committed
108
109
110
111
112

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

    m_XMLparser = new XMLSettingsReader();

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

122
123
    m_alignment = new Alignment();

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

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

    m_XMLparser = new XMLSettingsReader();

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

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

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

    for (unsigned int i = 0; i < m_XMLparser->getBaseNodeCount("channel"); i++) {
167
        Channel *buffer_ch = new Channel();
168
169
170
171
        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
172
        if(m_runNumber < first_run || m_runNumber > last_run) continue;
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
189
190
191
192
        bool isNew(true);
        for( int k = 0; k < channelEntries.size(); k++){
            if(buffer_ch->name == channelEntries.at(k)->name){
                std::cout << "WARNING!!! Redundancy in your settings file for " << buffer_ch->name << ". Check it carefully. The second entry found will be skipped..." << std::endl;
                isNew = false;
            }
        }
        if(isNew) channelEntries.push_back(buffer_ch);
193
194
195
    }

    std::cout << "Loaded " << channelEntries.size() << " configuration entries " << std::endl;
196
    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;
197
198
199
200
    ZDC* zdc1 = new ZDC(channelEntries,1);
    ZDC* zdc2 = new ZDC(channelEntries,2);
    RPD* rpd = new RPD(channelEntries);

201
202
203
204
    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

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

207
    return;
208
209
}

210
211
212
213
214
215
216
217
218
219
220
221
/**
 * @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
222
223
224
    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);
225

Riccardo Longo's avatar
Riccardo Longo committed
226
227
    std::cout << "WARNING: detector recognition glitch. NULL pointer being returned..." << std::endl;
    return NULL;
228
229
}

230

231
232
233
234
235
236
237
238
239
240
241
242
243
244
/** @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();
}

245
246
247
248
249
250
251
252
253
254
255
256
/** @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 ){
257
258
259
260
261
262
263
264
265
266
267
268
269

    // 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
     */
270
  } else {
Yakov Kulinich's avatar
Yakov Kulinich committed
271
    m_fIn = TFile::Open( m_fNameIn.c_str() );
272
  }
Yakov Kulinich's avatar
Yakov Kulinich committed
273
274
275
276
277
278
279
280

  // 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" );
281
  
Yakov Kulinich's avatar
Yakov Kulinich committed
282
283
284
285
  for( auto& ana : m_ana ){
    ana->Initialize();
    ana->SetupHistograms();
  }
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
}



/** @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
304
305
  
  //Added for debug/demo purposes. Remove after implementation of Visualizer
306
  TCanvas *canvas = new TCanvas( "Diff Demo", "Diff Demo", 200, 10, 1000, 600);
307
308
  TPad *pad = new TPad("pad", "pad",0.15,0.11,0.85,0.79);
  canvas->Divide(4,2);
clantz's avatar
clantz committed
309
  
310
311
312
313
314
315
  std::string style = "ATLAS";
  //Visualizer vis;
  std::vector< TH1* > raw;
  std::vector< TH1* > diff;
  
  
316
  // Processed Raw data to read in as vector of vectors size NxM
Yakov Kulinich's avatar
Yakov Kulinich committed
317
  // Where N = nCh and M = nSamples per channel.
318
319
  std::vector< std::vector< float >  >  vWF;
  std::vector< std::vector< float >* > pvWF;
Yakov Kulinich's avatar
Yakov Kulinich committed
320
321
322
323
324
325

  // 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.
326
327
328
  vWF .resize( m_nCh );
  pvWF.resize( m_nCh );
  vWFH.resize( m_nCh );
Yakov Kulinich's avatar
Yakov Kulinich committed
329

330
331
332
  /** 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
333
  TTree* tree = static_cast< TTree* >( m_fIn->Get( "tree" ) );
334

335
336
337
338
339
340
341
342
343
344
345
346
347
  //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.
348
  // Other items should be implemented in the same way if needed.
349
350
  // Also create the histograms to fill
  for( uint ch = 0; ch < m_nCh; ch++ ){
351
    //Example - here we retrieve the already processed waveform (M.Phipps approach during the data production)
352
    pvWF[ ch ] = &vWF[ ch ];
353
    tree->SetBranchAddress( Form( "C%d", ch ), &pvWF[ ch ] );
Yakov Kulinich's avatar
Yakov Kulinich committed
354
    vWFH[ ch ] = new TH1D( Form( "hWF%d", ch ), Form( "hWF%d;samp;amp", ch ), m_nSamp, 0, m_nSamp );
355
  }
356
357

  std::cout << "File: " << m_fIn->GetName() << " has " << tree->GetEntries() << " events." << std::endl;
358
359
  
  // !! EVENT LOOP
360
361
  for( int ev = 0; ev < tree->GetEntries(); ev++ ){
    
362
    tree->GetEntry( ev );
363
    
364
365
366
    // Fill the raw waveforms
    for( uint detID = 0; detID < (int) m_detectors.size(); detID++ )
        m_detectors.at(detID)->FillHistograms();
367

368
   //Here if you're interested in already processed waveform
369
   for( uint ch = 0; ch < m_nCh; ch++ ) {
370
   // Loop over samples in each channel
371
372
373
374
375
    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
376
    // Now call all analysis and run their AnalyzeEvent.
377
378
    // Can either send a vector of histos, or a 2D vector
    // Of all the data, depending on what you want to do.
379
    // Note that at the moment none of this methods is doing anything
Yakov Kulinich's avatar
Yakov Kulinich committed
380
    for( auto& ana : m_ana ){
381
      //raw data analysis
clantz's avatar
clantz committed
382
383
384
385
386
387
388
389
390
      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() );
391
392
393
394
395
        for(int test = 0; test<16; test ++){
          raw.push_back( rpd->GetChannelsVector().at(test)->WF_histo);
          diff.push_back( rpd->GetChannelsVector().at(test)->FirstDerivative);
        }
        //vis.ManyPadsPlot(raw,diff,4,4,Form("Event %d" ,ev),"Overlay");
clantz's avatar
clantz committed
396
397
398
399
      }else{
      ana->AnalyzeEvent( zdc1->GetChannelsVector() );
      ana->AnalyzeEvent( zdc2->GetChannelsVector() );
      ana->AnalyzeEvent( rpd->GetChannelsVector()  );
400
      //already processed wf analysis
401
402
      //ana->AnalyzeEvent( vWFH );
      //ana->AnalyzeEvent( vWF  );
clantz's avatar
clantz committed
403
      }
Yakov Kulinich's avatar
Yakov Kulinich committed
404
    }
405
  } // End event loop
406
  
clantz's avatar
clantz committed
407
408
409
410
411
412
413
  
  //Added for debug/demo purposes. Remove after implementation of Visualizer
  if(m_debug){
    pad->Update();
    canvas->Draw();
    canvas->Print( "Output.pdf" );
  }
414
415
416
417
418
419
420
421
422
423
424
425

  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
426
  
427
428
429
430
  if( m_fIn ){
    m_fIn->Close();
  }

Yakov Kulinich's avatar
Yakov Kulinich committed
431
432
433
434
435
436
437
438
439
  // 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();
440
441
}