ARA ROOT Test BEd Software

analysis/makeL2Files.cxx

//Includes
#include <iostream>

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

//ROOT Includes
#include "TChain.h"
#include "TTree.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH1D.h"
#include "TH2.h"
#include "TH2D.h"

//Global input file variables
TChain *feventTreeIn=0;
RawAraEvent *evPtr=0;

//Global variables used for the analysis
TH1D *fftWindow[RFCHANS_PER_STATION];
TH2D *fftAvTime[RFCHANS_PER_STATION];
Int_t numEventsInWindow[RFCHANS_PER_STATION]={0};
Double_t lastEventRate=0;
Int_t eventNumLastHour=0;
Double_t lastUnixTime=0;
Double_t lastUnixTimeHours[RFCHANS_PER_STATION]={0};
Double_t unixTimeFirstEvent=0;


void makeHistos();

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

int openRootFile(char *fileDir, int minRun, int maxRun);

double getXMax(TGraph *gr);

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

  if(argc<5) {
    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 minRun maxRun output.root\n";
    return 0;
  }

  if(openRootFile(argv[1], atoi(argv[2]), atoi(argv[3]))==0){
    std::cout << "Opened input file\n";
  }
  else return -1;

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

  char histoName[260];
  Int_t numBinsF=256;
  Double_t minF=0.0;
  Double_t maxF=1000.0;
  minF = minF - ( (maxF-minF)/numBinsF/2.0 ); // adjust histogram edges so that the bin centers
  maxF = maxF + ( (maxF-minF)/numBinsF/2.0 ); // of the histograms are aligned with graph [add bdf]
  numBinsF++;
  
  Int_t numBinsT=1001;
  Double_t minT=0.0;
  Double_t maxT=1000.0;
  
  
  for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
    sprintf(histoName, "average_fft_antenna_%i", antenna);
    fftWindow[antenna] = new TH1D(histoName, histoName, numBinsF, minF, maxF);
    sprintf(histoName, "average_vTime_fft_antenna_%i", antenna);
    fftAvTime[antenna] = new TH2D(histoName, histoName, numBinsT, minT, maxT, numBinsF, minF, maxF);
    
  }

  TTree *L1Tree = new TTree("L1Tree", "L1Tree");

  // Create variables to store in our output tree
  Double_t rms[RFCHANS_PER_STATION]={0};
  Double_t rmsFft[RFCHANS_PER_STATION]={0};
  Double_t peakFreq[RFCHANS_PER_STATION]={0};
  Double_t eventRate=0;
  Int_t eventNum=0;

  Int_t lastEventNum=0;
  Int_t runNumber=0;//FIXME
  Int_t unixTime=0; //FIXME  could use full unix time in us as well

  UsefulAraEvent *realEvent=0;
  RawAraHeader *head=0;
  TGraph *rfChan;
  TGraph *rfChanFft;
  TH1D* tempFft;


  L1Tree->Branch("rms", rms, "rms[16]/D");//FIXME
  L1Tree->Branch("rmsFft", rmsFft, "rmsFft[16]/D");
  L1Tree->Branch("peakFreq", peakFreq, "peakFreq[16]/D");
  L1Tree->Branch("eventRate", &eventRate, "eventRate/D");
  L1Tree->Branch("eventNum", &eventNum, "eventNum/I");
  L1Tree->Branch("runNumber", &runNumber, "runNumber/I");
  L1Tree->Branch("unixTime", &unixTime, "unixTime/I");
  //  L1Tree->Branch("realEvent", "UsefulAraEvent", &realEvent);
  //  L1Tree->Branch("head", "AraEventHeader", &head);


  //Now loop through the event Tree

   Long64_t numEntries=feventTreeIn->GetEntries();

   Long64_t starEvery=numEntries/80;
   if(starEvery==0) starEvery++;
   
   //   numEntries=10;//FIXME


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

     feventTreeIn->GetEntry(event);
     realEvent = new UsefulAraEvent(evPtr,AraCalType::kLatestCalib);
     eventNum = evPtr->head.eventNumber;
     unixTime = evPtr->head.unixTime;
     head = &(evPtr->head);

     eventRate=calcEventRate(lastUnixTime, unixTime, lastEventNum, eventNum, 60, lastEventRate);

     for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
       
       //std::cerr << "antenna " << antenna << "\n";
       rfChan = realEvent->getGraphFromRFChan(antenna);
       rfChanFft = realEvent->getFFTForRFChan(antenna);
       tempFft = realEvent->getFFTHistForRFChan(antenna);
       rms[antenna]=rfChan->GetRMS(2);
       rmsFft[antenna]=rfChanFft->GetRMS(2);
       Double_t tempPeak=0;
       FFTtools::getPeakRmsSqVal(rfChan, tempPeak, rmsFft[antenna]);
       peakFreq[antenna]=getXMax(rfChanFft);


       //std::cerr << "Populated the boring stuff\n";

       //Fill our histograms for the average FFT

       if(unixTime<lastUnixTimeHours[antenna]+60){
         //This case we must add the data to our TH1D    
         //std::cerr << "Adding to the TH1\n";
         fftWindow[antenna]->Add(tempFft);
         numEventsInWindow[antenna]++;
         //std::cerr << "Added to the TH1\n";
       }

       else{
         //This case we must extract the information and fill the TH2, then clear the TH1

         //std::cerr << "populating the TH2\n";

         //first scale the TH1

         if(numEventsInWindow[antenna]>0) 
           fftWindow[antenna]->Scale(1./numEventsInWindow[antenna]);
         
         //Then fill the TH2
         for(int bin=1; bin<fftWindow[antenna]->GetNbinsX()+1;bin++){
           Double_t binContent = fftWindow[antenna]->GetBinContent(bin);
           Double_t binCenter =  fftWindow[antenna]->GetBinCenter(bin);
           fftAvTime[antenna]->Fill((unixTime-unixTimeFirstEvent)/60, binCenter, binContent);
         }
         //Now clear the TH1
         fftWindow[antenna]->Reset();

         numEventsInWindow[antenna]=0;
         lastUnixTimeHours[antenna]=unixTime;

       }


       //std::cerr << "deleting graphs and TH1s\n";
       delete rfChan;
       delete rfChanFft;
       delete tempFft;

     }

     L1Tree->Fill();

     if(unixTime>lastUnixTime+60){
       lastEventNum = eventNum;
       lastUnixTime = unixTime;
       lastEventRate = eventRate;
     }
     

     delete realEvent;

   }

   //Now make sure that we save everything and close the files

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

   std::cerr << "\n";



}

void makeHistos(){

  char histoName[260];
   Int_t numBinsF=256;
   Double_t minF=0.0;
   Double_t maxF=1000.0;
   minF = minF - ( (maxF-minF)/numBinsF/2.0 ); // adjust histogram edges so that the bin centers
   maxF = maxF + ( (maxF-minF)/numBinsF/2.0 ); // of the histograms are aligned with graph [add bdf]
   numBinsF++;

   Int_t numBinsT=1000;
   Double_t minT=0.0;
   Double_t maxT=10.0;


  for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
    sprintf(histoName, "average_fft_antenna_%i", antenna);
    fftWindow[antenna] = new TH1D(histoName, histoName, numBinsF, minF, maxF);
    sprintf(histoName, "average_vTime_fft_antenna_%i", antenna);
    fftAvTime[antenna] = new TH2D(histoName, histoName, numBinsT, minT, maxT, numBinsF, minF, maxF);

  }


}

void saveHistos(){

  for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
    fftAvTime[antenna]->Write();

  }


}


int openRootFile(char *fileDir, int minRun, int maxRun){

  feventTreeIn = new TChain("eventTree");
  char name[260];
  for(int run=minRun;run<maxRun+1;run++){
    sprintf(name, "%s/run0%i/event0%i.root",fileDir, run,run);
    //    printf("Opening %s\n", name);
    if(feventTreeIn->Add(name)==0)
      printf("No entries in the tree\n");
  }

  feventTreeIn->SetBranchAddress("event",&evPtr);
  feventTreeIn->GetEntry(0);
  unixTimeFirstEvent=evPtr->head.unixTime;

  for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
    lastUnixTimeHours[antenna]=unixTimeFirstEvent;
  }


  return 0;

}


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];


}

 

---

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