ARA ROOT Test BEd Software

analysis/makeL2Files_For_Testing.cxx

//Includes
#include <iostream>

//AraRoot Includes
#include "UsefulAraEvent.h"
#include "RawAraEvent.h"
#include "RawAraHeader.h"
#include "AraHkData.h"
#include "araDefines.h"
#include "FFTtools.h"
#include "RFWindow.h"


//ROOT Includes
#include "TTree.h"
#include "TFile.h"
#include "TGraph.h"

double calcEventRate(double lastTime, double thisTime, double lastEvNum, double thisEvNum, double interval, double lastRate);

double getXMax(TGraph *gr);

int main(int argc, char **argv)
{

  if(argc<3) {
    std::cout << "Usage\n" << argv[0] << " <input file> <output file>/n";
    std::cout << "e.g.\n" << argv[0] << " http://www.hep.ucl.ac.uk/uhen/ara/monitor/root/run1841/event1841.root output.root\n";
    return 0;
  }

  TFile *fpIn = TFile::Open(argv[1]);
  if(!fpIn) {
    std::cerr << "Can't open file\n";
     return -1;
   }
   TTree *eventTree = (TTree*) fpIn->Get("eventTree");
   if(!eventTree) {
     std::cerr << "Can't find eventTree\n";
     return -1;
   }

   RawAraEvent *evPtr=0;
   eventTree->SetBranchAddress("event",&evPtr);


   char outFileName[FILENAME_MAX];
   sprintf(outFileName, "%s", argv[2]);   //FIXME
   TFile *fpOut = new TFile(outFileName, "RECREATE");
   TTree *L1Tree = new TTree("L1Tree", "L1Tree");

   //Create some useful variables
   Double_t rms[RFCHANS_PER_STATION]={0};
   Double_t rmsFft[RFCHANS_PER_STATION]={0};
   Double_t peakFreq[RFCHANS_PER_STATION]={0};
   Double_t power[RFCHANS_PER_STATION]={0};
   Double_t eventRate=0;
   Double_t eventRate60=0;
   Double_t lastEventRate=0;
   Double_t lastEventRate60=0;
   Double_t lastTime=0;
   Double_t firstTime=0;
   Int_t lastEventNum=0;
   Int_t eventNum=0;
   Int_t runNumber=0;//FIXME
   Int_t unixTime=0; //FIXME  could use full unix time in us as well

   // TGraph *rfChan[RFCHANS_PER_STATION]={0};
   // TGraph *rfChanFft[RFCHANS_PER_STATION]={0};
   
   RawAraHeader *head=0;
   AraHkData *hk=0;
   UsefulAraEvent *realEvent=0;
   RFWindow *window=0;

   L1Tree->Branch("rms", rms, "rms[16]/D");//FIXME
   L1Tree->Branch("rmsFft", rmsFft, "rmsFft[16]/D");
   L1Tree->Branch("peakFreq", peakFreq, "peakFreq[16]/D");
   L1Tree->Branch("power", power, "power[16]/D");//FIXME
   L1Tree->Branch("eventRate", &eventRate, "eventRate/D");
   L1Tree->Branch("eventRate60", &eventRate60, "eventRate60/D");
   L1Tree->Branch("lastTime", &lastTime, "lastTime/D");
   L1Tree->Branch("firstTime", &firstTime, "firstTime/D");
   L1Tree->Branch("lastEventNum", &lastEventNum, "lastEventNum/I");

   L1Tree->Branch("eventNum", &eventNum, "eventNum/I");
   L1Tree->Branch("runNumber", &runNumber, "runNumber/I");
   L1Tree->Branch("unixTime", &unixTime, "unixTime/I");
   
   // L1Tree->Branch("rfChan","TGraph[16]",rfChan);
   // L1Tree->Branch("rfChanFft","TGraph[16]",rfChanFft);

   L1Tree->Branch("realEvent", "UsefulAraEvent", &realEvent);
   L1Tree->Branch("head","RawAraHeader",&head);
   L1Tree->Branch("hk","AraHkData",&hk);
   L1Tree->Branch("window","RFWindow",&window);

   //Now loop through the tree

   TGraph *rfChan;
   TGraph *rfChanFft;
   window = new RFWindow();

   Long64_t numEntries=eventTree->GetEntries();
   Long64_t starEvery=numEntries/80;
   if(starEvery==0) starEvery++;


   for(Long64_t event=0;event<numEntries;event++) {
     if(event%starEvery==0) {
       std::cerr << "*";       
     }

     eventTree->GetEntry(event);
     realEvent = new UsefulAraEvent(evPtr,AraCalType::kLatestCalib);

     eventNum = evPtr->head.eventNumber;
     unixTime = evPtr->head.unixTime;
     if(event==0) firstTime=unixTime;
     head=&(evPtr->head);
     hk=&(evPtr->hk);
     //     eventRate=calcEventRate(lastTime, unixTime, lastEventNum, eventNum, 1, lastEventRate);
     eventRate60=calcEventRate(lastTime, unixTime, lastEventNum, eventNum, 60, lastEventRate60);
     window->addGraphToWindow(realEvent->getFFTForRFChan(0));

     for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
       
       // rfChan[antenna] = realEvent->getGraphFromRFChan(antenna);
       // rfChanFft[antenna] = realEvent->getFFTForRFChan(antenna);
       // rms[antenna]=rfChan[antenna]->GetRMS(2);
       // Double_t tempPeak=0;
       // FFTtools::getPeakRmsSqVal(rfChan[antenna], tempPeak, rmsFft[antenna]);
       // peakFreq[antenna]=FFTtools::getPeakVal(rfChanFft[antenna]);
       // power[antenna]=0; //FIXME
       rfChan = realEvent->getGraphFromRFChan(antenna);
       rfChanFft = realEvent->getFFTForRFChan(antenna);
       rms[antenna]=rfChan->GetRMS(2);
       Double_t tempPeak=0;
       FFTtools::getPeakRmsSqVal(rfChan, tempPeak, rmsFft[antenna]);
       peakFreq[antenna]=getXMax(rfChanFft);
       power[antenna]=0; //FIXME

       delete rfChan;
       delete rfChanFft;
       

     }
     L1Tree->Fill();
     if(unixTime>lastTime+60){
       lastEventNum = eventNum;
       lastTime = unixTime;
       lastEventRate = eventRate;
       lastEventRate60 = eventRate60;
     }

     // for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
     //   delete rfChan[antenna];
     //   delete rfChanFft[antenna];
     // }

     delete realEvent;
   
   }
   
   delete window;
   
   fpIn->Close();

   fpOut->Write();
   fpOut->Close();

   std::cerr << "\n";

}

double calcEventRate(double lastTime, double thisTime, double lastEvNum, double thisEvNum, double interval, double lastRate){
     Double_t dRate = 0;
     Double_t eventRate = 0;
     if(thisEvNum > lastEvNum)
       {
         Double_t dEventNum = thisEvNum - lastEvNum;
         if(thisTime > lastTime+interval)
           dRate = (thisTime-lastTime)/dEventNum;
       }
     if(dRate>0)
       eventRate = 1./dRate;
     else 
       eventRate = lastRate;
  
     return eventRate;

}

double getXMax(TGraph *gr){
  int bin=FFTtools::getPeakBin(gr);
  return gr->GetX()[bin];


}


//Things to Do

//Average FFT over an n event window

---

UCL HEP Home   

UHE Neutrino Home   

ARA ROOT Software

Main  

Download  

Branches

3.13  

3.12  

3.11  

3.10  

3.9  

3.8  

3.7  

3.6  

3.5  

3.4  

3.3  

3.2  

TestBed