ARA ROOT Test BEd Software

analysis/makeL2Files_depricated.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"
//Include FFTtools.h if you want to ask the correlation, etc. tools

//ROOT Includes
#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
TFile *fpIn=0;
TTree *feventTreeIn=0;
RawAraEvent *evPtr=0;
RFWindow *window=0;

//Some more basic global variables
TH1D *fftWindow[RFCHANS_PER_STATION];
TH2D *fftAvTime[RFCHANS_PER_STATION];



void makeHistos();

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

int openRootFile(char *fileName);

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

  if(openRootFile(argv[1])==0){
    std::cout << "Opened input file\n";
  }
  else
    return -1;
  
   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 lastEventRate=0;
   Int_t eventNum=0;
   Int_t eventNumLastHour=0;
   Int_t lastEventNum=0;
   Int_t runNumber=0;//FIXME
   Int_t unixTime=0; //FIXME  could use full unix time in us as well
   Double_t lastUnixTime=0;
   Double_t lastUnixTimeHours=0;

   RawAraHeader *head=0;
   AraHkData *hk=0;
   UsefulAraEvent *realEvent=0;
   TGraph *rfChan;
   TGraph *rfChanFft;
   // TGraph *rfChan[RFCHANS_PER_STATION]={0};
   // TGraph *rfChanFft[RFCHANS_PER_STATION]={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("lastUnixTime", &lastUnixTime, "lastUnixTime/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);
   
   window = new RFWindow();
   L1Tree->Branch("window","RFWindow",&window);

   //Now loop through the 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);
     hk=&(evPtr->hk);

     eventRate=calcEventRate(lastUnixTime, unixTime, lastEventNum, eventNum, 60, lastEventRate);
   
     for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
       
       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
      
       //decide whether to add to the rf window
       TH1D* tempFftHisto = realEvent->getFFTHistForRFChan(antenna);
       if(unixTime>lastUnixTimeHours+3600){
         window->addToAvHisto(antenna, tempFftHisto);
       }
       else
         window->clearAvHisto(antenna);
       
       //Fill our TH1D and TH2D for the average FFt
       TH1D *tempTempFftHisto = realEvent->getFFTHistForRFChan(antenna);
       if(unixTime<lastUnixTimeHours+3600){
         fftWindow->Add(tempTempFftHisto);
         eventNumLastHour++;
       }
       else{
         fftAvTime->Fill


       }
         



       delete rfChan;
       delete rfChanFft;
       delete tempFftHisto;
       

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

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

     delete realEvent;
   
   }
   
   fpIn->Close();

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

   std::cerr << "\n";

}

void makeHistos(){

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

   Int_t numBinsX=11;
   Double_t minX=0.0;
   Double_t maxX=10.0;


  for(int antenna=0;antenna<RFCHANS_PER_STATION;antenna++){
    sprintf(histoName, "average_fft_antenna_%i", antenna);
    fftWindow[antenna] = new TH1D(histoName, histoName, numBinsX, minX, maxX);
    fftAvTime[antenna] = new TH1D(histoName, histoName, numBinsX, minX, maxX, numBinsY, minY, maxY);

  }


}


int openRootFile(char *fileName){

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

   feventTreeIn->SetBranchAddress("event",&evPtr);

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


}

   // //Create TH2 to store the average FFT with the right number of bins
   // feventTreeIn->GetEntry(0);
   // Double_t minX=evPtr->head.unixTime/60;
   // feventTreeIn->GetEntry(feventTreeIn->GetEntries()-1);
   // Double_t maxX=evPtr->head.unixTime/60;
   // Int_t xBins=(maxX-minX)/60;
   // Double_t maxY=1000.0;
   // Double_t minY = 0.0;
   // Int_t yBins = 256;
   // minY = minY - ( (maxY-minY)/yBins/2.0 ); // adjust histogram edges so that the bin centers
   // maxY = maxY + ( (maxY-minY)/yBins/2.0 ); // of the histograms are aligned with graph [add bdf]
   // yBins++;

   // printf("First time %i, lastTime %i number of bins is %i\n", firstTime, endTime, nTimeBins);

   // avFFT = new TH2D("average_FFT_histo", "average_FFT_histo", xBins, minX, maxX, yBins, minY, maxY);









//Things to Do

//Average FFT over an hour

// --> Need 16 TH2D's to store that are x - time y - frequency
// --> Need 16 TH1D's to fill this
// ----> Each TH1D is added to over an hour then averaged at the end


// Some functions like:

// addHisto for antenna
// getAvHisto for antenna
// clearAvHisto

---

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