ARA ROOT Test BEd Software

analysis/makeFFTHistos.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 "TObject.h"
#include "TObjArray.h"
#include "TList.h"
#include "TChain.h"
#include "TTree.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH1D.h"
#include "TH2.h"
#include "TH2D.h"
#include "TMath.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};
Int_t eventNumLastHour=0;
Long64_t lastUnixTime=0;
Long64_t lastUnixTimeHours[RFCHANS_PER_STATION]={0};
Long64_t unixTimeFirstEvent=0;
Long64_t unixTimeOfFirstHour=0;
Int_t eventNumFirstEvent=0;

void makeHistos();

double calcEventRate(Long64_t lastTime, Long64_t 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++;
  
  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);
    
  }


  // Create variables to store in our output tree
  Int_t eventNum=0;
  Int_t lastEventNum=0;
  Int_t unixTime=0; //FIXME  could use full unix time in us as well
  Long64_t histoTime=0;
  Double_t eventRate=0;
  Double_t lastEventRate=0;
  Double_t rms[RFCHANS_PER_STATION]={0};
  Int_t eventsInHisto[RFCHANS_PER_STATION]={0};

  TTree *mapTree = new TTree("mapTimes", "mapTimes");
  int fillTree=0;
  mapTree->Branch("histoTime", &histoTime, "histoTime/L");
  mapTree->Branch("eventRate", &eventRate, "eventRate/D");
  mapTree->Branch("rms", rms, "rms[16]/D");//FIXME
  mapTree->Branch("eventsInHisto", &eventsInHisto, "eventsInHisto[16]/I");

  UsefulAraEvent *realEvent=0;
  TH1D* tempFft;
  TGraph *rfChan;


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




     for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
       rfChan = realEvent->getGraphFromRFChan(antenna);
       rms[antenna]=rfChan->GetRMS(2);
       tempFft = realEvent->getFFTHistForRFChan(antenna);       

       //Fill our histograms for the average FFT
       if(unixTime<lastUnixTimeHours[antenna]+60*60){
         for(int bin=1; bin<(tempFft->GetNbinsX()+1);bin++){
           double expBin = TMath::Power(10,(tempFft->GetBinContent(bin))/10);
           fftWindow[antenna]->AddBinContent(bin, expBin);
           // fftWindow[antenna]->Add(tempFft);
         }
       }

       else{
         if(numEventsInWindow[antenna]>0) {

           fftWindow[antenna]->Scale(1./numEventsInWindow[antenna]);
           //re-scale to the dB scale

           for(int bin=1; bin<(fftWindow[antenna]->GetNbinsX()+1);bin++){

             double logBin=0;
             if(fftWindow[antenna]>0)
               logBin = 10*TMath::Log10(fftWindow[antenna]->GetBinContent(bin));
               
             fftWindow[antenna]->SetBinContent(bin, logBin);
             // fftWindow[antenna]->Add(tempFft);
           }

         }

         sprintf(histoName, "time_%llu", lastUnixTimeHours[antenna]+60*60);
         if(fpOut->GetDirectory(histoName)==0)
           fpOut->mkdir(histoName);
         fpOut->cd(histoName);
         fftWindow[antenna]->Write();
         fpOut->cd();

         //Now clear the TH1
         fftWindow[antenna]->Reset();

         eventsInHisto[antenna]=numEventsInWindow[antenna];
         numEventsInWindow[antenna]=0;
         histoTime=lastUnixTimeHours[antenna]+60*60;
         lastUnixTimeHours[antenna]=histoTime;

         fillTree=1;

       }
       
       //std::cerr << "deleting graphs and TH1s\n";

       delete rfChan;
       delete tempFft;

     }
     
     if(fillTree){
         eventRate=calcEventRate(lastUnixTime, unixTime, lastEventNum, eventNum, 3600, lastEventRate);

         //      printf("eventRate %e lastUnixTime %llu unixTime %llu lastEventNum %d eventNum %d lastEventRate %e\n",eventRate, lastUnixTime, unixTime, lastEventNum, eventNum, lastEventRate);
         lastUnixTime=histoTime;
         lastEventNum=eventNum;
         lastEventRate=eventRate;

       mapTree->Fill();
     }
     fillTree=0;

     delete realEvent;

   }

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

   //make sure that we save the tree in the right directory
   mapTree->Write();
   
   fpOut->cd();

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

  }


}

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/run%06i/event%06i.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;
  unixTimeOfFirstHour=((unixTimeFirstEvent/3600)*3600);

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

  return 0;

}


double calcEventRate(Long64_t lastTime, Long64_t 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