ARA ROOT Test BEd Software

analysis/pulserResponse.cxx

//pulserResponse.cxx - Make histos of peak cross corr time between pairs of channels for pulser events

//Includes
#include <iostream>
#include <fstream>
#include <cmath>

//AraRoot Includes
#include "UsefulAraEvent.h"
#include "RawAraEvent.h"
#include "AraEventCalibrator.h"
#include "FFTtools.h"


//ROOT Includes
#include "TTree.h"
#include "TFile.h"
#include "TGraph.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TChain.h"


using namespace std;


int pulserResponse(char rootDir[260], int minRun, int maxRun, char outDir[260], char nameOptions[260], int calibType)
{

  char outName[260];
  if(calibType==0)
    sprintf(outName, "%s/pulserResponse_%s_kLatestCalib.root", outDir, nameOptions);
  if(calibType==1)
    sprintf(outName, "%s/pulserResponse_%s_kFirstCalibPlusCables.root", outDir, nameOptions);
  if(calibType==2)
    sprintf(outName, "%s/pulserResponse_%s_kFirstCalib.root", outDir, nameOptions);
  if(calibType==3)
    sprintf(outName, "%s/pulserResponse_%s_kVoltageTime.root", outDir, nameOptions);
  
  TFile *fOut = TFile::Open(outName, "RECREATE");

  TH1D *histCorrelation = new TH1D("pulserCrossCorr", "pulserCrossCorr", 1000, -100, 0);

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

  Double_t deltaT=0;
  correlationTree->Branch("deltaT", &deltaT, "deltaT/D");

  TChain *eventTree = new TChain("eventTree");
  
  for(int i=minRun; i<maxRun+1; i++){
    char inName[260];
    
    sprintf(inName, "%s/run%06d/event%06d.root", rootDir, i, i);
    eventTree->Add(inName);
  }
  
  RawAraEvent *evPtr=0;
   
  eventTree->SetBranchAddress("event",&evPtr); 
    
  Long64_t numEntries=eventTree->GetEntries();
  Long64_t starEvery=numEntries/80;
  if(starEvery==0) starEvery++;
  UsefulAraEvent *realEvent;
  
  cout << numEntries << endl;
  for(Long64_t entry=0; entry<numEntries; entry++){
   
    if(entry%starEvery==0){
      cerr << "*";
    }
    eventTree->GetEntry(entry);


    
    if(calibType==0)
      realEvent = new UsefulAraEvent(evPtr, AraCalType::kLatestCalib);  
    else if(calibType==1)
      realEvent = new UsefulAraEvent(evPtr, AraCalType::kFirstCalibPlusCables);  
    else if(calibType==2)
      realEvent = new UsefulAraEvent(evPtr, AraCalType::kFirstCalib);  
    else if(calibType==3)
      realEvent = new UsefulAraEvent(evPtr, AraCalType::kVoltageTime);
    else
      continue;

 if(realEvent->isCalPulserEvent()==0)
      {
        delete realEvent;
        continue;
      }

      TGraph *chanA = realEvent->getGraphFromRFChan(2);
      TGraph *chanB = realEvent->getGraphFromRFChan(4);

      double *fvoltsA=chanA->GetY();
      double *fvoltsB=chanB->GetY();    
  
    //zero mean the waveform
    //check that the waveforms are sizeable i.e. peak > 300mV

      double maxA=0, maxB=0;

      int nSamples = chanA->GetN();
      double offset=0;
      for(int i=0; i<nSamples; i++){
        offset+=fvoltsA[i]/nSamples;
      }
      for(int i=0; i<nSamples; i++){
        fvoltsA[i]-=offset;
        if(fvoltsA[i]>maxA){
          maxA=fvoltsA[i];
        }
      }

      nSamples=chanB->GetN();
      offset=0;
      for(int i=0; i<nSamples;i++){
        offset+=fvoltsB[i]/nSamples;
      }
      for(int i=0; i<nSamples; i++){
        fvoltsB[i]-=offset;
        if(fvoltsB[i]>maxB){
          maxB=fvoltsB[i];
        }
      }
      if(maxB<300||maxA<300){
        delete realEvent;
        continue;
      }

      //zero mean waveform

      //interpolate the graphs
      //pad to length

      TGraph *chanAint = FFTtools::getInterpolatedGraph(chanA, 0.05);
      TGraph *chanBint = FFTtools::getInterpolatedGraph(chanB, 0.05);

      //correlate the graphs

      TGraph *correlTemp = FFTtools::getCorrelationGraph(chanAint, chanBint);

      double *xCorrel = correlTemp->GetX();
      double *yCorrel = correlTemp->GetY();

      double xMax = 0, yMax = 0;

      //find max correlation value and add to array
      //made some changes to what is stored in histogram
      //now stores entire correlated waveform

      nSamples = correlTemp->GetN();
      for(int i=0; i<nSamples; i++){

        //      histCorrelation->Fill(xCorrel[i], yCorrel[i]);

        if(yCorrel[i]>yMax){
          yMax=yCorrel[i];
          xMax=xCorrel[i];
        }
      }

      deltaT=xMax;
      histCorrelation->Fill(xMax);
      correlationTree->Fill();

      delete chanA;
      delete chanB;
      delete chanAint;
      delete chanBint;
      delete correlTemp;     
      delete realEvent;
      
  }
  //fOut->cd();
  //  histCorrelation->Write();
  fOut->WriteTObject(histCorrelation);
  fOut->WriteTObject(correlationTree);
  fOut->Close();

  
  std::cerr << "\n";
  return 0;
  
}

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

  if(argc<7){
    cout << "Incorrect usage\n";
    cout << argv[0] << " <rootDir> <minRun> <maxRun> <outDir> <nameOptions> <calibType 0=kLatestCalib 1=kFirstCalibPlusCables 2=kFirstCalib 3=kVoltageTime>\n";
    return -1;
  }

  if(atoi(argv[6])<0||atoi(argv[6])>3)
    {
      cout << "Incorrect usage\n";
      cout << argv[0] << " <rootDir> <minRun> <maxRun> <outDir> <nameOptions> <calibType 0=kLatestCalib 1=kFirstCalibPlusCables 2=kFirstCalib 3=kVoltageTime>\n";      
      return -1;
    }

  

  pulserResponse(argv[1], atoi(argv[2]), atoi(argv[3]), argv[4], argv[5], atoi(argv[6]));

  return 0;
}

---

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